/[escript]/trunk/finley/src/CPPAdapter/MeshAdapterFactory.cpp
ViewVC logotype

Diff of /trunk/finley/src/CPPAdapter/MeshAdapterFactory.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/esys2/finley/src/CPPAdapter/MeshAdapterFactory.cpp revision 102 by jgs, Wed Dec 15 07:08:39 2004 UTC trunk/finley/src/CPPAdapter/MeshAdapterFactory.cpp revision 1753 by ksteube, Sun Sep 7 22:01:23 2008 UTC
# Line 1  Line 1 
1    
2  /* $Id$ */  /* $Id$ */
 /*  
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  
  *                                                                            *  
  * This software is the property of ACcESS. No part of this code              *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that person has a software    *  
  * license agreement with ACcESS.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
3    
4    /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16    #ifdef PASO_MPI
17    #include <mpi.h>
18    #endif
19    #ifdef USE_NETCDF
20    #include <netcdfcpp.h>
21    #endif
22    #include "MeshAdapterFactory.h"
23    #include "FinleyError.h"
24  extern "C" {  extern "C" {
25  #include "finley/finleyC/Finley.h"  #include "escript/blocktimer.h"
 #include "finley/finleyC/Mesh.h"  
 #include "finley/finleyC/RectangularMesh.h"  
26  }  }
 #include "finley/CPPAdapter/FinleyError.h"  
 #include "finley/CPPAdapter/MeshAdapterFactory.h"  
   
27    
28    #include <boost/python/extract.hpp>
29    
 #include <iostream>  
30  #include <sstream>  #include <sstream>
31    
32  using namespace std;  using namespace std;
# Line 31  using namespace escript; Line 34  using namespace escript;
34    
35  namespace finley {  namespace finley {
36    
37    #ifdef USE_NETCDF
38      // A convenience method to retrieve an integer attribute from a NetCDF file
39      int NetCDF_Get_Int_Attribute(NcFile *dataFile, char *fName, char *attr_name) {
40        NcAtt *attr;
41        char error_msg[LenErrorMsg_MAX];
42        if (! (attr=dataFile->get_att(attr_name)) ) {
43          sprintf(error_msg,"Error retrieving integer attribute '%s' from NetCDF file '%s'", attr_name, fName);
44          throw DataException(error_msg);
45        }
46        int temp = attr->as_int(0);
47        delete attr;
48        return(temp);
49      }
50    #endif
51    
52      AbstractContinuousDomain* loadMesh(const std::string& fileName)
53      {
54    #ifdef USE_NETCDF
55        bool optimize=FALSE; // Don't optimize since this would cause problems with Data().dump()
56        Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
57        AbstractContinuousDomain* temp;
58        Finley_Mesh *mesh_p=NULL;
59        char error_msg[LenErrorMsg_MAX];
60    
61        char *fName = Paso_MPI_appendRankToFileName(fileName.c_str(),
62                                                    mpi_info->size,
63                                                    mpi_info->rank);
64    
65        double blocktimer_start = blocktimer_time();
66        Finley_resetError();
67    
68        // Open NetCDF file for reading
69        NcAtt *attr;
70        NcVar *nc_var_temp;
71        // netCDF error handler
72        NcError err(NcError::silent_nonfatal);
73        // Create the NetCDF file.
74        NcFile dataFile(fName, NcFile::ReadOnly);
75        if (!dataFile.is_valid()) {
76          sprintf(error_msg,"loadMesh: Opening file NetCDF %s for reading failed.", fName);
77          Finley_setError(IO_ERROR,error_msg);
78          Paso_MPIInfo_free( mpi_info );
79          throw DataException(error_msg);
80        }
81    
82        // Read NetCDF integer attributes
83        int mpi_size            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_size");
84        int mpi_rank            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_rank");
85        int numDim              = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numDim");
86        int order               = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"order");
87        int reduced_order           = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"reduced_order");
88        int numNodes            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numNodes");
89        int num_Elements            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements");
90        int num_FaceElements        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements");
91        int num_ContactElements     = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_ContactElements");
92        int num_Points          = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Points");
93        int num_Elements_numNodes       = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements_numNodes");
94        int Elements_TypeId         = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Elements_TypeId");
95        int num_FaceElements_numNodes   = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements_numNodes");
96        int FaceElements_TypeId     = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"FaceElements_TypeId");
97        int num_ContactElements_numNodes    = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_ContactElements_numNodes");
98        int ContactElements_TypeId      = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"ContactElements_TypeId");
99        int Points_TypeId           = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Points_TypeId");
100        int num_Tags            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Tags");
101    
102        // Verify size and rank
103        if (mpi_info->size != mpi_size) {
104          sprintf(error_msg, "Error loadMesh - The NetCDF file '%s' can only be read on %d CPUs instead of %d", fName, mpi_size, mpi_info->size);
105          throw DataException(error_msg);
106        }
107        if (mpi_info->rank != mpi_rank) {
108          sprintf(error_msg, "Error loadMesh - The NetCDF file '%s' should be read on CPU #%d instead of %d", fName, mpi_rank, mpi_info->rank);
109          throw DataException(error_msg);
110        }
111    
112        // Read mesh name
113        if (! (attr=dataFile.get_att("Name")) ) {
114          sprintf(error_msg,"Error retrieving mesh name from NetCDF file '%s'", fName);
115          throw DataException(error_msg);
116        }
117        char *name = attr->as_string(0);
118        delete attr;
119    
120        /* allocate mesh */
121        mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
122        if (Finley_noError()) {
123    
124            /* read nodes */
125            Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
126        // Nodes_Id
127            if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
128              throw DataException("Error - loadMesh:: unable to read Nodes_Id from netCDF file: " + *fName);
129            if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {
130              free(&mesh_p->Nodes->Id);
131              throw DataException("Error - loadMesh:: unable to recover Nodes_Id from NetCDF file: " + *fName);
132            }
133        // Nodes_Tag
134            if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
135              throw DataException("Error - loadMesh:: unable to read Nodes_Tag from netCDF file: " + *fName);
136            if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) ) {
137              free(&mesh_p->Nodes->Tag);
138              throw DataException("Error - loadMesh:: unable to recover Nodes_Tag from NetCDF file: " + *fName);
139            }
140        // Nodes_gDOF
141            if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
142              throw DataException("Error - loadMesh:: unable to read Nodes_gDOF from netCDF file: " + *fName);
143            if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) ) {
144              free(&mesh_p->Nodes->globalDegreesOfFreedom);
145              throw DataException("Error - loadMesh:: unable to recover Nodes_gDOF from NetCDF file: " + *fName);
146            }
147        // Nodes_gNI
148            if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
149              throw DataException("Error - loadMesh:: unable to read Nodes_gNI from netCDF file: " + *fName);
150            if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) ) {
151              free(&mesh_p->Nodes->globalNodesIndex);
152              throw DataException("Error - loadMesh:: unable to recover Nodes_gNI from NetCDF file: " + *fName);
153            }
154        // Nodes_grDfI
155            if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
156              throw DataException("Error - loadMesh:: unable to read Nodes_grDfI from netCDF file: " + *fName);
157            if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) ) {
158              free(&mesh_p->Nodes->globalReducedDOFIndex);
159              throw DataException("Error - loadMesh:: unable to recover Nodes_grDfI from NetCDF file: " + *fName);
160            }
161        // Nodes_grNI
162            if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
163              throw DataException("Error - loadMesh:: unable to read Nodes_grNI from netCDF file: " + *fName);
164            if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) ) {
165              free(&mesh_p->Nodes->globalReducedNodesIndex);
166              throw DataException("Error - loadMesh:: unable to recover Nodes_grNI from NetCDF file: " + *fName);
167            }
168        // Nodes_Coordinates
169            if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates"))) {
170              free(&mesh_p->Nodes->Coordinates);
171              throw DataException("Error - loadMesh:: unable to read Nodes_Coordinates from netCDF file: " + *fName);
172            }
173            if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) ) {
174              free(&mesh_p->Nodes->Coordinates);
175              throw DataException("Error - load:: unable to recover Nodes_Coordinates from netCDF file: " + *fName);
176            }
177        // Nodes_DofDistribution
178        int *first_component = TMPMEMALLOC(mpi_size+1,index_t);
179            if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) )
180              throw DataException("Error - loadMesh:: unable to read Nodes_DofDistribution from netCDF file: " + *fName);
181            if (! nc_var_temp->get(&first_component[0], mpi_size+1) ) {
182              free(&first_component);
183              throw DataException("Error - loadMesh:: unable to recover Nodes_DofDistribution from NetCDF file: " + *fName);
184            }
185        mesh_p->Nodes->degreesOfFreedomDistribution=Paso_Distribution_alloc(mesh_p->Nodes->MPIInfo,first_component,1,0);
186        TMPMEMFREE(first_component);
187    
188            /* read elements */
189            if (Finley_noError()) {
190              mesh_p->Elements=Finley_ElementFile_alloc((ElementTypeId)Elements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
191              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
192              mesh_p->Elements->minColor=0;
193              mesh_p->Elements->maxColor=num_Elements-1;
194              if (num_Elements>0) {
195                     if (Finley_noError()) {
196                   // Elements_Id
197                       if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
198                         throw DataException("Error - loadMesh:: unable to read Elements_Id from netCDF file: " + *fName);
199                       if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) ) {
200                         free(&mesh_p->Elements->Id);
201                         throw DataException("Error - loadMesh:: unable to recover Elements_Id from NetCDF file: " + *fName);
202                       }
203                   // Elements_Tag
204                       if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
205                         throw DataException("Error - loadMesh:: unable to read Elements_Tag from netCDF file: " + *fName);
206                       if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) ) {
207                         free(&mesh_p->Elements->Tag);
208                         throw DataException("Error - loadMesh:: unable to recover Elements_Tag from NetCDF file: " + *fName);
209                       }
210                   // Elements_Owner
211                       if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
212                         throw DataException("Error - loadMesh:: unable to read Elements_Owner from netCDF file: " + *fName);
213                       if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) ) {
214                         free(&mesh_p->Elements->Owner);
215                         throw DataException("Error - loadMesh:: unable to recover Elements_Owner from NetCDF file: " + *fName);
216                       }
217                   // Elements_Color
218                       if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
219                         throw DataException("Error - loadMesh:: unable to read Elements_Color from netCDF file: " + *fName);
220                       if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) ) {
221                         free(&mesh_p->Elements->Color);
222                         throw DataException("Error - loadMesh:: unable to recover Elements_Color from NetCDF file: " + *fName);
223                       }
224                   // Elements_Nodes
225               int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
226                       if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
227                         free(&mesh_p->Elements->Nodes);
228                         throw DataException("Error - loadMesh:: unable to read Elements_Nodes from netCDF file: " + *fName);
229                       }
230                       if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
231                         free(&Elements_Nodes);
232                         throw DataException("Error - load:: unable to recover Elements_Nodes from netCDF file: " + *fName);
233                       }
234               // Copy temp array into mesh_p->Elements->Nodes
235               for (int i=0; i<num_Elements; i++) {
236                 for (int j=0; j<num_Elements_numNodes; j++) {
237                   mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
238                     = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
239                 }
240               }
241               TMPMEMFREE(Elements_Nodes);
242             }
243          } /* num_Elements>0 */
244        }
245    
246            /* get the face elements */
247            if (Finley_noError()) {
248              mesh_p->FaceElements=Finley_ElementFile_alloc((ElementTypeId)FaceElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
249              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
250              mesh_p->FaceElements->minColor=0;
251              mesh_p->FaceElements->maxColor=num_FaceElements-1;
252              if (num_FaceElements>0) {
253                     if (Finley_noError()) {
254                   // FaceElements_Id
255                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
256                         throw DataException("Error - loadMesh:: unable to read FaceElements_Id from netCDF file: " + *fName);
257                       if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) ) {
258                         free(&mesh_p->FaceElements->Id);
259                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Id from NetCDF file: " + *fName);
260                       }
261                   // FaceElements_Tag
262                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
263                         throw DataException("Error - loadMesh:: unable to read FaceElements_Tag from netCDF file: " + *fName);
264                       if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) ) {
265                         free(&mesh_p->FaceElements->Tag);
266                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Tag from NetCDF file: " + *fName);
267                       }
268                   // FaceElements_Owner
269                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
270                         throw DataException("Error - loadMesh:: unable to read FaceElements_Owner from netCDF file: " + *fName);
271                       if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) ) {
272                         free(&mesh_p->FaceElements->Owner);
273                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Owner from NetCDF file: " + *fName);
274                       }
275                   // FaceElements_Color
276                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
277                         throw DataException("Error - loadMesh:: unable to read FaceElements_Color from netCDF file: " + *fName);
278                       if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) ) {
279                         free(&mesh_p->FaceElements->Color);
280                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Color from NetCDF file: " + *fName);
281                       }
282                   // FaceElements_Nodes
283               int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
284                       if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
285                         free(&mesh_p->FaceElements->Nodes);
286                         throw DataException("Error - loadMesh:: unable to read FaceElements_Nodes from netCDF file: " + *fName);
287                       }
288                       if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
289                         free(&FaceElements_Nodes);
290                         throw DataException("Error - load:: unable to recover FaceElements_Nodes from netCDF file: " + *fName);
291                       }
292               // Copy temp array into mesh_p->FaceElements->Nodes
293               for (int i=0; i<num_FaceElements; i++) {
294                 for (int j=0; j<num_FaceElements_numNodes; j++) {
295                   mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]
296                     = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
297                 }
298               }
299               TMPMEMFREE(FaceElements_Nodes);
300             }
301          } /* num_FaceElements>0 */
302        }
303    
304            /* get the Contact elements */
305            if (Finley_noError()) {
306              mesh_p->ContactElements=Finley_ElementFile_alloc((ElementTypeId)ContactElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
307              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->ContactElements, num_ContactElements);
308              mesh_p->ContactElements->minColor=0;
309              mesh_p->ContactElements->maxColor=num_ContactElements-1;
310              if (num_ContactElements>0) {
311                     if (Finley_noError()) {
312                   // ContactElements_Id
313                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Id")) )
314                         throw DataException("Error - loadMesh:: unable to read ContactElements_Id from netCDF file: " + *fName);
315                       if (! nc_var_temp->get(&mesh_p->ContactElements->Id[0], num_ContactElements) ) {
316                         free(&mesh_p->ContactElements->Id);
317                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Id from NetCDF file: " + *fName);
318                       }
319                   // ContactElements_Tag
320                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Tag")) )
321                         throw DataException("Error - loadMesh:: unable to read ContactElements_Tag from netCDF file: " + *fName);
322                       if (! nc_var_temp->get(&mesh_p->ContactElements->Tag[0], num_ContactElements) ) {
323                         free(&mesh_p->ContactElements->Tag);
324                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Tag from NetCDF file: " + *fName);
325                       }
326                   // ContactElements_Owner
327                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Owner")) )
328                         throw DataException("Error - loadMesh:: unable to read ContactElements_Owner from netCDF file: " + *fName);
329                       if (! nc_var_temp->get(&mesh_p->ContactElements->Owner[0], num_ContactElements) ) {
330                         free(&mesh_p->ContactElements->Owner);
331                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Owner from NetCDF file: " + *fName);
332                       }
333                   // ContactElements_Color
334                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Color")) )
335                         throw DataException("Error - loadMesh:: unable to read ContactElements_Color from netCDF file: " + *fName);
336                       if (! nc_var_temp->get(&mesh_p->ContactElements->Color[0], num_ContactElements) ) {
337                         free(&mesh_p->ContactElements->Color);
338                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Color from NetCDF file: " + *fName);
339                       }
340                   // ContactElements_Nodes
341               int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
342                       if (!(nc_var_temp = dataFile.get_var("ContactElements_Nodes"))) {
343                         free(&mesh_p->ContactElements->Nodes);
344                         throw DataException("Error - loadMesh:: unable to read ContactElements_Nodes from netCDF file: " + *fName);
345                       }
346                       if (! nc_var_temp->get(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes) ) {
347                         free(&ContactElements_Nodes);
348                         throw DataException("Error - load:: unable to recover ContactElements_Nodes from netCDF file: " + *fName);
349                       }
350               // Copy temp array into mesh_p->ContactElements->Nodes
351               for (int i=0; i<num_ContactElements; i++) {
352                 for (int j=0; j<num_ContactElements_numNodes; j++) {
353                   mesh_p->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]
354                     = ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
355                 }
356               }
357               TMPMEMFREE(ContactElements_Nodes);
358             }
359          } /* num_ContactElements>0 */
360        }
361    
362            /* get the Points (nodal elements) */
363            if (Finley_noError()) {
364              mesh_p->Points=Finley_ElementFile_alloc((ElementTypeId)Points_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
365              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Points, num_Points);
366              mesh_p->Points->minColor=0;
367              mesh_p->Points->maxColor=num_Points-1;
368              if (num_Points>0) {
369                 if (Finley_noError()) {
370                   // Points_Id
371                       if (! ( nc_var_temp = dataFile.get_var("Points_Id")) )
372                         throw DataException("Error - loadMesh:: unable to read Points_Id from netCDF file: " + *fName);
373                       if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points) ) {
374                         free(&mesh_p->Points->Id);
375                         throw DataException("Error - loadMesh:: unable to recover Points_Id from NetCDF file: " + *fName);
376                       }
377                   // Points_Tag
378                       if (! ( nc_var_temp = dataFile.get_var("Points_Tag")) )
379                         throw DataException("Error - loadMesh:: unable to read Points_Tag from netCDF file: " + *fName);
380                       if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points) ) {
381                         free(&mesh_p->Points->Tag);
382                         throw DataException("Error - loadMesh:: unable to recover Points_Tag from NetCDF file: " + *fName);
383                       }
384                   // Points_Owner
385                       if (! ( nc_var_temp = dataFile.get_var("Points_Owner")) )
386                         throw DataException("Error - loadMesh:: unable to read Points_Owner from netCDF file: " + *fName);
387                       if (! nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points) ) {
388                         free(&mesh_p->Points->Owner);
389                         throw DataException("Error - loadMesh:: unable to recover Points_Owner from NetCDF file: " + *fName);
390                       }
391                   // Points_Color
392                       if (! ( nc_var_temp = dataFile.get_var("Points_Color")) )
393                         throw DataException("Error - loadMesh:: unable to read Points_Color from netCDF file: " + *fName);
394                       if (! nc_var_temp->get(&mesh_p->Points->Color[0], num_Points) ) {
395                         free(&mesh_p->Points->Color);
396                         throw DataException("Error - loadMesh:: unable to recover Points_Color from NetCDF file: " + *fName);
397                       }
398                   // Points_Nodes
399               int *Points_Nodes = TMPMEMALLOC(num_Points,int);
400                       if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
401                         free(&mesh_p->Points->Nodes);
402                         throw DataException("Error - loadMesh:: unable to read Points_Nodes from netCDF file: " + *fName);
403                       }
404                       if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
405                         free(&Points_Nodes);
406                         throw DataException("Error - load:: unable to recover Points_Nodes from netCDF file: " + *fName);
407                       }
408               // Copy temp array into mesh_p->Points->Nodes
409               for (int i=0; i<num_Points; i++) {
410                 mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
411               }
412               TMPMEMFREE(Points_Nodes);
413             }
414          } /* num_Points>0 */
415        }
416    
417            /* get the tags */
418            if (Finley_noError()) {
419              if (num_Tags>0) {
420                // Temp storage to gather node IDs
421                int *Tags_keys = TMPMEMALLOC(num_Tags, int);
422                char name_temp[4096];
423            int i;
424    
425            // Tags_keys
426                if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) )
427                  throw DataException("Error - loadMesh:: unable to read Tags_keys from netCDF file: " + *fName);
428                if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
429                  free(&Tags_keys);
430                  throw DataException("Error - loadMesh:: unable to recover Tags_keys from NetCDF file: " + *fName);
431                }
432            for (i=0; i<num_Tags; i++) {
433                  // Retrieve tag name
434                  sprintf(name_temp, "Tags_name_%d", i);
435                  if (! (attr=dataFile.get_att(name_temp)) ) {
436                    sprintf(error_msg,"Error retrieving tag name from NetCDF file '%s'", fName);
437                    throw DataException(error_msg);
438                  }
439                  char *name = attr->as_string(0);
440                  delete attr;
441                  Finley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
442            }
443          }
444        }
445    
446        } /* Finley_noError() after Finley_Mesh_alloc() */
447      
448        if (Finley_noError()) Finley_Mesh_createMappings(mesh_p, mesh_p->Nodes->degreesOfFreedomDistribution->first_component, 0);
449    
450        checkFinleyError();
451        temp=new MeshAdapter(mesh_p);
452    
453        if (! Finley_noError()) {
454          Finley_Mesh_free(mesh_p);
455        }
456    
457        /* win32 refactor */
458        TMPMEMFREE(fName);
459    
460        blocktimer_increment("LoadMesh()", blocktimer_start);
461        return temp;
462    #else
463        throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
464    #endif /* USE_NETCDF */
465      }
466    
467    AbstractContinuousDomain* readMesh(const std::string& fileName,    AbstractContinuousDomain* readMesh(const std::string& fileName,
468                       int integrationOrder)                       int integrationOrder,
469                                         int reducedIntegrationOrder,
470                                         int optimize)
471    {    {
472      //      //
473      // create a copy of the filename to overcome the non-constness of call      // create a copy of the filename to overcome the non-constness of call
474      // to Finley_Mesh_read      // to Finley_Mesh_read
475      char fName[fileName.size()+1];      Finley_Mesh* fMesh=0;
476        // Win32 refactor
477        if( fileName.size() == 0 )
478        {
479           throw DataException("Null file name!");
480        }
481    
482        char *fName = TMPMEMALLOC(fileName.size()+1,char);
483        
484      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
485      Finley_Mesh* fMesh=Finley_Mesh_read(fName,integrationOrder);      double blocktimer_start = blocktimer_time();
486    
487        fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
488      checkFinleyError();      checkFinleyError();
489      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
490        
491        /* win32 refactor */
492        TMPMEMFREE(fName);
493        
494        blocktimer_increment("ReadMesh()", blocktimer_start);
495        return temp;
496      }
497    
498      AbstractContinuousDomain* readMeshMPI(const std::string& fileName,
499                         int integrationOrder,
500                                         int reducedIntegrationOrder,
501                                         int optimize)
502      {
503        //
504        // create a copy of the filename to overcome the non-constness of call
505        // to Finley_Mesh_read
506        Finley_Mesh* fMesh=0;
507        // Win32 refactor
508        if( fileName.size() == 0 )
509        {
510           throw DataException("Null file name!");
511        }
512    
513        char *fName = TMPMEMALLOC(fileName.size()+1,char);
514        
515        strcpy(fName,fileName.c_str());
516        double blocktimer_start = blocktimer_time();
517    
518        fMesh=Finley_Mesh_read_MPI(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
519        checkFinleyError();
520        AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
521        
522        /* win32 refactor */
523        TMPMEMFREE(fName);
524        
525        blocktimer_increment("ReadMesh()", blocktimer_start);
526        return temp;
527      }
528    
529      AbstractContinuousDomain* readGmsh(const std::string& fileName,
530                                         int numDim,
531                                         int integrationOrder,
532                                         int reducedIntegrationOrder,
533                                         int optimize)
534      {
535        //
536        // create a copy of the filename to overcome the non-constness of call
537        // to Finley_Mesh_read
538        Finley_Mesh* fMesh=0;
539        // Win32 refactor
540        if( fileName.size() == 0 )
541        {
542           throw DataException("Null file name!");
543        }
544    
545        char *fName = TMPMEMALLOC(fileName.size()+1,char);
546        
547        strcpy(fName,fileName.c_str());
548        double blocktimer_start = blocktimer_time();
549    
550        fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
551        checkFinleyError();
552        AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
553        
554        /* win32 refactor */
555        TMPMEMFREE(fName);
556        
557        blocktimer_increment("ReadGmsh()", blocktimer_start);
558      return temp;      return temp;
559    }    }
560    
# Line 50  namespace finley { Line 563  namespace finley {
563              int periodic0,int periodic1,              int periodic0,int periodic1,
564              int periodic2,              int periodic2,
565              int integrationOrder,              int integrationOrder,
566              int useElementsOnFace)                      int reducedIntegrationOrder,
567                int useElementsOnFace,
568                int useFullElementOrder,
569                        int optimize)
570    {    {
 //     cout << "n0=" << n0 << " n1=" << n1 << " n2=" << n2  
 //   << " order=" << order  
 //   << " l0=" << l0 << " l1=" << l1 << " l2=" << l2  
 //   << " periodic0=" << periodic0  
 //   << " periodic1=" << periodic1  
 //   << " periodic2=" << periodic2  
 //   << " integerationOrder=" << integrationOrder  
 //   << " useElementsOnFace=" << useElementsOnFace << endl;  
         
571      int numElements[]={n0,n1,n2};      int numElements[]={n0,n1,n2};
572      double length[]={l0,l1,l2};      double length[]={l0,l1,l2};
573      int periodic[]={periodic0, periodic1, periodic2};      int periodic[]={periodic0, periodic1, periodic2};
574    
575      //      //
576      // linearInterpolation      // linearInterpolation
577      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=NULL;
578    
579      if (order==1) {      if (order==1) {
580        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
581                      useElementsOnFace) ;                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
582      } else if (order==2) {      }
583        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,          else if (order==2) {
584                       useElementsOnFace) ;        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
585                         useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
586      } else {      } else {
587        stringstream temp;        stringstream temp;
588        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
# Line 89  namespace finley { Line 598  namespace finley {
598              double l0, double l1,              double l0, double l1,
599              int periodic0,int periodic1,              int periodic0,int periodic1,
600              int integrationOrder,              int integrationOrder,
601              int useElementsOnFace)                          int reducedIntegrationOrder,
602                int useElementsOnFace,
603                    int useFullElementOrder,
604                            int optimize)
605    {    {
606      int numElements[]={n0,n1};      int numElements[]={n0,n1};
607      double length[]={l0,l1};      double length[]={l0,l1};
608      int periodic[]={periodic0, periodic1};      int periodic[]={periodic0, periodic1};
609    
610      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=0;
611      if (order==1) {      if (order==1) {
612        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
613                      useElementsOnFace);                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
614      } else if (order==2) {      }
615        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,      else if (order==2) {
616                      useElementsOnFace);        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
617      } else {                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
618        }
619        else {
620        stringstream temp;        stringstream temp;
621        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
622        setFinleyError(VALUE_ERROR,temp.str().c_str());        setFinleyError(VALUE_ERROR,temp.str().c_str());
# Line 113  namespace finley { Line 627  namespace finley {
627      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
628      return temp;      return temp;
629    }    }
630    AbstractContinuousDomain*  interval(int n0,int order,double l0,int periodic0,  
631                 int integrationOrder,    AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
632                 int useElementsOnFace)    {
633    {      Finley_Mesh* fMesh=0;
634      int numElements[]={n0};      //
635      double length[]={l0};      // extract the meshes from meshList
636      int periodic[]={periodic0};      int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
637      Finley_Mesh* fMesh;      Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
638      if (order==1) {      for (int i=0;i<numMsh;++i) {
639        fMesh=Finley_RectangularMesh_Line2(numElements, length,periodic,integrationOrder,           AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
640                       useElementsOnFace);           const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
641      } else if (order==2) {           mshes[i]=finley_meshListMember->getFinley_Mesh();
       fMesh=Finley_RectangularMesh_Line3(numElements,length,periodic,integrationOrder,  
                      useElementsOnFace);  
     } else {  
       stringstream temp;  
       temp << "Illegal interpolation order: " << order;  
       setFinleyError(VALUE_ERROR,temp.str().c_str());  
642      }      }
643      //      //
644        // merge the meshes:
645        fMesh=Finley_Mesh_merge(numMsh,mshes);
646          TMPMEMFREE(mshes);
647        //
648      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
649      checkFinleyError();      checkFinleyError();
650      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
651      return temp;  
   }  
   AbstractContinuousDomain*  meshMerge(const boost::python::list& meshList)  
   {  
     AbstractContinuousDomain* temp=new MeshAdapter(0);  
652      return temp;      return temp;
653    }    }
654    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,
655              double safetyFactor,                                     double safety_factor,
656              double tolerance)                             double tolerance,
657                                           int optimize)
658    {    {
659      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
660      return temp;      //
661        // merge the meshes:
662        AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
663        //
664        // glue the faces:
665        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
666        fMesh=merged_finley_meshes->getFinley_Mesh();
667        Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
668    
669        //
670        // Convert any finley errors into a C++ exception
671        checkFinleyError();
672        return merged_meshes;
673    }    }
674    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,
675              double safety_factor,              double safety_factor,
676              double tolerance)              double tolerance,
677                            int optimize)
678    {    {
679      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
680      return temp;      //
681        // merge the meshes:
682        AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
683        //
684        // join the faces:
685        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
686        fMesh=merged_finley_meshes->getFinley_Mesh();
687        Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
688        //
689        // Convert any finley errors into a C++ exception
690        checkFinleyError();
691        return merged_meshes;
692    }    }
693    
694  }  // end of namespace    // end of namespace
695    
696    }

Legend:
Removed from v.102  
changed lines
  Added in v.1753

  ViewVC Help
Powered by ViewVC 1.1.26