/[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

revision 1345 by ksteube, Wed Nov 14 07:53:34 2007 UTC revision 2779 by caltinay, Thu Nov 26 03:51:15 2009 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2009 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  #ifdef PASO_MPI  #ifdef PASO_MPI
16  #include <mpi.h>  #include <mpi.h>
# Line 22  Line 21 
21  #include "MeshAdapterFactory.h"  #include "MeshAdapterFactory.h"
22  #include "FinleyError.h"  #include "FinleyError.h"
23  extern "C" {  extern "C" {
24  #include "escript/blocktimer.h"  #include "esysUtils/blocktimer.h"
25  }  }
26    
27  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
# Line 49  namespace finley { Line 48  namespace finley {
48    }    }
49  #endif  #endif
50    
51    AbstractContinuousDomain* loadMesh(const std::string& fileName)  //   AbstractContinuousDomain* loadMesh(const std::string& fileName)
52      Domain_ptr loadMesh(const std::string& fileName)
53    {    {
54  #ifdef USE_NETCDF  #ifdef USE_NETCDF
55      Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );      Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
56      AbstractContinuousDomain* temp;      AbstractContinuousDomain* temp;
57      Finley_Mesh *mesh_p=NULL;      Finley_Mesh *mesh_p=NULL;
58      char error_msg[LenErrorMsg_MAX];      char error_msg[LenErrorMsg_MAX];
     // create a copy of the filename to overcome the non-constness of call  
     // to Finley_Mesh_read  
     // Win32 refactor  
     char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;  
     strcpy(fName,fileName.c_str());  
59    
60      printf("ksteube finley::loadMesh %s\n", fName);      char *fName = Paso_MPI_appendRankToFileName(fileName.c_str(),
61                                                    mpi_info->size,
62                                                    mpi_info->rank);
63    
64      double blocktimer_start = blocktimer_time();      double blocktimer_start = blocktimer_time();
65      Finley_resetError();      Finley_resetError();
66        int *first_DofComponent, *first_NodeComponent;
67    
68      // Open NetCDF file for reading      // Open NetCDF file for reading
69      NcAtt *attr;      NcAtt *attr;
# Line 78  namespace finley { Line 76  namespace finley {
76        sprintf(error_msg,"loadMesh: Opening file NetCDF %s for reading failed.", fName);        sprintf(error_msg,"loadMesh: Opening file NetCDF %s for reading failed.", fName);
77        Finley_setError(IO_ERROR,error_msg);        Finley_setError(IO_ERROR,error_msg);
78        Paso_MPIInfo_free( mpi_info );        Paso_MPIInfo_free( mpi_info );
79        throw DataException("Error - loadMesh:: Could not read NetCDF file.");        throw DataException(error_msg);
80      }      }
81    
82      // Read NetCDF integer attributes      // Read NetCDF integer attributes
83      int mpi_size        = NetCDF_Get_Int_Attribute(&dataFile, fName, "mpi_size");      int mpi_size            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_size");
84      int mpi_rank        = NetCDF_Get_Int_Attribute(&dataFile, fName, "mpi_rank");      int mpi_rank            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_rank");
85      int numDim          = NetCDF_Get_Int_Attribute(&dataFile, fName, "numDim");      int numDim              = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numDim");
86      int order           = NetCDF_Get_Int_Attribute(&dataFile, fName, "order");      int order               = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"order");
87      int reduced_order       = NetCDF_Get_Int_Attribute(&dataFile, fName, "reduced_order");      int reduced_order           = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"reduced_order");
88      int numNodes        = NetCDF_Get_Int_Attribute(&dataFile, fName, "numNodes");      int numNodes            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numNodes");
89      int num_Elements        = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_Elements");      int num_Elements            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements");
90      int num_FaceElements    = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_FaceElements");      int num_FaceElements        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements");
91      int num_ContactElements = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_ContactElements");      int num_ContactElements     = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_ContactElements");
92      int num_Points      = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_Points");      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      // Read mesh name
113      if (! (attr=dataFile.get_att("Name")) ) {      if (! (attr=dataFile.get_att("Name")) ) {
# Line 101  namespace finley { Line 117  namespace finley {
117      char *name = attr->as_string(0);      char *name = attr->as_string(0);
118      delete attr;      delete attr;
119    
120        string msgPrefix("Error in loadMesh: NetCDF operation failed - ");
121    
122      /* allocate mesh */      /* allocate mesh */
123      mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);      mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
124      if (Finley_noError()) {      if (Finley_noError()) {
125    
126          /* read nodes */          /* read nodes */
127          Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);          Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
128      // Nodes_Id          // Nodes_Id
129          if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )          if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
130            throw DataException("Error - loadMesh:: unable to read Nodes_Id from netCDF file: " + *fName);            throw DataException(msgPrefix+"get_var(Nodes_Id)");
131          if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {          if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {
132            free(&mesh_p->Nodes->Id);            TMPMEMFREE(mesh_p->Nodes->Id);
133            throw DataException("Error - loadMesh:: unable to recover Nodes_Id from NetCDF file: " + *fName);            throw DataException("get(Nodes_Id)");
134          }          }
135  // printf("ksteube Nodes_Id: "); for (int i=0; i<numNodes; i++) { printf(" %d", mesh_p->Nodes->Id[i]); } printf("\n");          // Nodes_Tag
     // Nodes_Tag  
136          if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )          if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
137            throw DataException("Error - loadMesh:: unable to read Nodes_Tag from netCDF file: " + *fName);            throw DataException("get_var(Nodes_Tag)");
138          if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) ) {          if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) ) {
139            free(&mesh_p->Nodes->Tag);            TMPMEMFREE(mesh_p->Nodes->Tag);
140            throw DataException("Error - loadMesh:: unable to recover Nodes_Tag from NetCDF file: " + *fName);            throw DataException("get(Nodes_Tag)");
141          }          }
142      // Nodes_gDOF          // Nodes_gDOF
143          if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )          if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
144            throw DataException("Error - loadMesh:: unable to read Nodes_gDOF from netCDF file: " + *fName);            throw DataException("get_var(Nodes_gDOF)");
145          if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) ) {          if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) ) {
146            free(&mesh_p->Nodes->globalDegreesOfFreedom);            TMPMEMFREE(mesh_p->Nodes->globalDegreesOfFreedom);
147            throw DataException("Error - loadMesh:: unable to recover Nodes_gDOF from NetCDF file: " + *fName);            throw DataException("get(Nodes_gDOF)");
148          }          }
149      // Nodes_gNI          // Nodes_gNI
150          if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )          if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
151            throw DataException("Error - loadMesh:: unable to read Nodes_gNI from netCDF file: " + *fName);            throw DataException("get_var(Nodes_gNI)");
152          if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) ) {          if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) ) {
153            free(&mesh_p->Nodes->globalNodesIndex);            TMPMEMFREE(mesh_p->Nodes->globalNodesIndex);
154            throw DataException("Error - loadMesh:: unable to recover Nodes_gNI from NetCDF file: " + *fName);            throw DataException("get(Nodes_gNI)");
155          }          }
156      // Nodes_grDfI          // Nodes_grDfI
157          if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )          if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
158            throw DataException("Error - loadMesh:: unable to read Nodes_grDfI from netCDF file: " + *fName);            throw DataException("get_var(Nodes_grDfI)");
159          if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) ) {          if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) ) {
160            free(&mesh_p->Nodes->globalReducedDOFIndex);            TMPMEMFREE(mesh_p->Nodes->globalReducedDOFIndex);
161            throw DataException("Error - loadMesh:: unable to recover Nodes_grDfI from NetCDF file: " + *fName);            throw DataException("get(Nodes_grDfI)");
162          }          }
163      // Nodes_grNI          // Nodes_grNI
164          if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )          if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
165            throw DataException("Error - loadMesh:: unable to read Nodes_grNI from netCDF file: " + *fName);            throw DataException("get_var(Nodes_grNI)");
166          if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) ) {          if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) ) {
167            free(&mesh_p->Nodes->globalReducedNodesIndex);            TMPMEMFREE(mesh_p->Nodes->globalReducedNodesIndex);
168            throw DataException("Error - loadMesh:: unable to recover Nodes_grNI from NetCDF file: " + *fName);            throw DataException("get(Nodes_grNI)");
169          }          }
170      // Nodes_Coordinates          // Nodes_Coordinates
171          if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates"))) {          if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates"))) {
172            free(&mesh_p->Nodes->Coordinates);            TMPMEMFREE(mesh_p->Nodes->Coordinates);
173            throw DataException("Error - loadMesh:: unable to read Nodes_Coordinates from netCDF file: " + *fName);            throw DataException("get_var(Nodes_Coordinates)");
174          }          }
175          if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) ) {          if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) ) {
176            free(&mesh_p->Nodes->Coordinates);            TMPMEMFREE(mesh_p->Nodes->Coordinates);
177            throw DataException("Error - load:: unable to recover Nodes_Coordinates from netCDF file: " + *fName);            throw DataException("get(Nodes_Coordinates)");
178            }
179            // Nodes_DofDistribution
180            first_DofComponent = TMPMEMALLOC(mpi_size+1,index_t);
181            if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) )
182              throw DataException("get_var(Nodes_DofDistribution)");
183            if (! nc_var_temp->get(&first_DofComponent[0], mpi_size+1) ) {
184              throw DataException("get(Nodes_DofDistribution)");
185            }
186    
187            // Nodes_NodeDistribution
188            first_NodeComponent = TMPMEMALLOC(mpi_size+1,index_t);
189            if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) )
190              throw DataException("get_var(Nodes_NodeDistribution)");
191            if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
192              throw DataException("get(Nodes_NodeDistribution)");
193          }          }
194    
 #if 0 /* Not yet...finish the above first */  
195          /* read elements */          /* read elements */
196          if (Finley_noError()) {          if (Finley_noError()) {
197               mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);            Finley_ReferenceElementSet  *refElements= Finley_ReferenceElementSet_alloc((ElementTypeId)Elements_TypeId,mesh_p->order, mesh_p->reduced_order);
198               if (Finley_noError()) {            if (Finley_noError())  {
199                   Finley_ElementFile_allocTable(mesh_p->Elements, numEle);                mesh_p->Elements=Finley_ElementFile_alloc(refElements, mpi_info);
200                   mesh_p->Elements->minColor=0;            }
201                   mesh_p->Elements->maxColor=numEle-1;            Finley_ReferenceElementSet_dealloc(refElements);
202                   if (Finley_noError()) {            if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
203           }            if (Finley_noError()) {
204           }                mesh_p->Elements->minColor=0;
205                  mesh_p->Elements->maxColor=num_Elements-1;
206                  if (num_Elements>0) {
207                       // Elements_Id
208                       if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
209                         throw DataException("get_var(Elements_Id)");
210                       if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) ) {
211                         TMPMEMFREE(mesh_p->Elements->Id);
212                         throw DataException("get(Elements_Id)");
213                       }
214                       // Elements_Tag
215                       if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
216                         throw DataException("get_var(Elements_Tag)");
217                       if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) ) {
218                         TMPMEMFREE(mesh_p->Elements->Tag);
219                         throw DataException("get(Elements_Tag)");
220                       }
221                       // Elements_Owner
222                       if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
223                         throw DataException("get_var(Elements_Owner)");
224                       if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) ) {
225                         TMPMEMFREE(mesh_p->Elements->Owner);
226                         throw DataException("get(Elements_Owner)");
227                       }
228                       // Elements_Color
229                       if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
230                         throw DataException("get_var(Elements_Color)");
231                       if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) ) {
232                         TMPMEMFREE(mesh_p->Elements->Color);
233                         throw DataException("get(Elements_Color)");
234                       }
235                       // Elements_Nodes
236                       int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
237                       if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
238                         TMPMEMFREE(mesh_p->Elements->Nodes);
239                         throw DataException("get_var(Elements_Nodes)");
240                       }
241                       if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
242                         TMPMEMFREE(Elements_Nodes);
243                         throw DataException("get(Elements_Nodes)");
244                       }
245                       // Copy temp array into mesh_p->Elements->Nodes
246                       for (int i=0; i<num_Elements; i++) {
247                           for (int j=0; j<num_Elements_numNodes; j++) {
248                               mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
249                               = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
250                           }
251                       }
252                       TMPMEMFREE(Elements_Nodes);
253                    
254                  } /* num_Elements>0 */
255              }
256      }      }
 #endif  
257    
258          /* get the face elements */          /* get the face elements */
259            if (Finley_noError()) {
260              Finley_ReferenceElementSet *refFaceElements=  Finley_ReferenceElementSet_alloc((ElementTypeId)FaceElements_TypeId ,mesh_p->order, mesh_p->reduced_order);
261              if (Finley_noError())  {
262                  mesh_p->FaceElements=Finley_ElementFile_alloc(refFaceElements, mpi_info);
263              }
264              Finley_ReferenceElementSet_dealloc(refFaceElements);  
265              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
266              if (Finley_noError()) {
267                  mesh_p->FaceElements->minColor=0;
268                  mesh_p->FaceElements->maxColor=num_FaceElements-1;
269                  if (num_FaceElements>0) {
270                       // FaceElements_Id
271                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
272                         throw DataException("get_var(FaceElements_Id)");
273                       if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) ) {
274                         TMPMEMFREE(mesh_p->FaceElements->Id);
275                         throw DataException("get(FaceElements_Id)");
276                       }
277                       // FaceElements_Tag
278                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
279                         throw DataException("get_var(FaceElements_Tag)");
280                       if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) ) {
281                         TMPMEMFREE(mesh_p->FaceElements->Tag);
282                         throw DataException("get(FaceElements_Tag)");
283                       }
284                       // FaceElements_Owner
285                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
286                         throw DataException("get_var(FaceElements_Owner)");
287                       if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) ) {
288                         TMPMEMFREE(mesh_p->FaceElements->Owner);
289                         throw DataException("get(FaceElements_Owner)");
290                       }
291                       // FaceElements_Color
292                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
293                         throw DataException("get_var(FaceElements_Color)");
294                       if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) ) {
295                         TMPMEMFREE(mesh_p->FaceElements->Color);
296                         throw DataException("get(FaceElements_Color)");
297                       }
298                       // FaceElements_Nodes
299                       int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
300                       if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
301                         TMPMEMFREE(mesh_p->FaceElements->Nodes);
302                         throw DataException("get_var(FaceElements_Nodes)");
303                       }
304                       if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
305                         TMPMEMFREE(FaceElements_Nodes);
306                         throw DataException("get(FaceElements_Nodes)");
307                       }
308                       // Copy temp array into mesh_p->FaceElements->Nodes
309                       for (int i=0; i<num_FaceElements; i++) {
310                            for (int j=0; j<num_FaceElements_numNodes; j++) {
311                                mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
312                            }
313                        }
314                        TMPMEMFREE(FaceElements_Nodes);
315                  } /* num_FaceElements>0 */
316              }
317        }
318    
319          /* get the Contact face element */          /* get the Contact elements */
320            if (Finley_noError()) {
321              Finley_ReferenceElementSet *refContactElements=   Finley_ReferenceElementSet_alloc((ElementTypeId)ContactElements_TypeId,mesh_p->order, mesh_p->reduced_order);
322              if (Finley_noError())  {
323                  mesh_p->ContactElements=Finley_ElementFile_alloc(refContactElements, mpi_info);
324              }
325              Finley_ReferenceElementSet_dealloc(refContactElements);  
326              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->ContactElements, num_ContactElements);
327              if (Finley_noError()) {
328                  mesh_p->ContactElements->minColor=0;
329                  mesh_p->ContactElements->maxColor=num_ContactElements-1;
330                  if (num_ContactElements>0) {
331                   // ContactElements_Id
332                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Id")) )
333                         throw DataException("get_var(ContactElements_Id)");
334                       if (! nc_var_temp->get(&mesh_p->ContactElements->Id[0], num_ContactElements) ) {
335                         TMPMEMFREE(mesh_p->ContactElements->Id);
336                         throw DataException("get(ContactElements_Id)");
337                       }
338                   // ContactElements_Tag
339                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Tag")) )
340                         throw DataException("get_var(ContactElements_Tag)");
341                       if (! nc_var_temp->get(&mesh_p->ContactElements->Tag[0], num_ContactElements) ) {
342                         TMPMEMFREE(mesh_p->ContactElements->Tag);
343                         throw DataException("get(ContactElements_Tag)");
344                       }
345                   // ContactElements_Owner
346                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Owner")) )
347                         throw DataException("get_var(ContactElements_Owner)");
348                       if (! nc_var_temp->get(&mesh_p->ContactElements->Owner[0], num_ContactElements) ) {
349                         TMPMEMFREE(mesh_p->ContactElements->Owner);
350                         throw DataException("get(ContactElements_Owner)");
351                       }
352                   // ContactElements_Color
353                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Color")) )
354                         throw DataException("get_var(ContactElements_Color)");
355                       if (! nc_var_temp->get(&mesh_p->ContactElements->Color[0], num_ContactElements) ) {
356                         TMPMEMFREE(mesh_p->ContactElements->Color);
357                         throw DataException("get(ContactElements_Color)");
358                       }
359                   // ContactElements_Nodes
360                        int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
361                       if (!(nc_var_temp = dataFile.get_var("ContactElements_Nodes"))) {
362                         TMPMEMFREE(mesh_p->ContactElements->Nodes);
363                         throw DataException("get_var(ContactElements_Nodes)");
364                       }
365                       if (! nc_var_temp->get(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes) ) {
366                         TMPMEMFREE(ContactElements_Nodes);
367                         throw DataException("get(ContactElements_Nodes)");
368                       }
369                       // Copy temp array into mesh_p->ContactElements->Nodes
370                        for (int i=0; i<num_ContactElements; i++) {
371                            for (int j=0; j<num_ContactElements_numNodes; j++) {
372                                mesh_p->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]= ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
373                            }
374                        }
375                        TMPMEMFREE(ContactElements_Nodes);
376                     } /* num_ContactElements>0 */
377                  }
378              
379        }
380    
381          /* get the nodal element */          /* get the Points (nodal elements) */
382            if (Finley_noError()) {
383              Finley_ReferenceElementSet *refPoints=    Finley_ReferenceElementSet_alloc((ElementTypeId)Points_TypeId,mesh_p->order, mesh_p->reduced_order);
384              if (Finley_noError())  {
385                  mesh_p->Points=Finley_ElementFile_alloc(refPoints, mpi_info);
386              }
387              Finley_ReferenceElementSet_dealloc(refPoints);
388              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Points, num_Points);
389              if (Finley_noError()) {
390                  mesh_p->Points->minColor=0;
391                  mesh_p->Points->maxColor=num_Points-1;
392                  if (num_Points>0) {
393                   // Points_Id
394                       if (! ( nc_var_temp = dataFile.get_var("Points_Id")) )
395                         throw DataException("get_var(Points_Id)");
396                       if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points) ) {
397                         TMPMEMFREE(mesh_p->Points->Id);
398                         throw DataException("get(Points_Id)");
399                       }
400                   // Points_Tag
401                       if (! ( nc_var_temp = dataFile.get_var("Points_Tag")) )
402                         throw DataException("get_var(Points_Tag)");
403                       if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points) ) {
404                         TMPMEMFREE(mesh_p->Points->Tag);
405                         throw DataException("get(Points_Tag)");
406                       }
407                   // Points_Owner
408                       if (! ( nc_var_temp = dataFile.get_var("Points_Owner")) )
409                         throw DataException("get_var(Points_Owner)");
410                       if (! nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points) ) {
411                         TMPMEMFREE(mesh_p->Points->Owner);
412                         throw DataException("get(Points_Owner)");
413                       }
414                   // Points_Color
415                       if (! ( nc_var_temp = dataFile.get_var("Points_Color")) )
416                         throw DataException("get_var(Points_Color)");
417                       if (! nc_var_temp->get(&mesh_p->Points->Color[0], num_Points) ) {
418                         TMPMEMFREE(mesh_p->Points->Color);
419                         throw DataException("get(Points_Color)");
420                       }
421                   // Points_Nodes
422               int *Points_Nodes = TMPMEMALLOC(num_Points,int);
423                       if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
424                         TMPMEMFREE(mesh_p->Points->Nodes);
425                         throw DataException("get_var(Points_Nodes)");
426                       }
427                       if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
428                         TMPMEMFREE(Points_Nodes);
429                         throw DataException("get(Points_Nodes)");
430                       }
431               // Copy temp array into mesh_p->Points->Nodes
432               for (int i=0; i<num_Points; i++) {
433                 mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
434               }
435               TMPMEMFREE(Points_Nodes);
436            
437                  } /* num_Points>0 */
438              }
439            }
440    
441          /* get the name tags */          /* get the tags */
442            if (Finley_noError()) {
443              if (num_Tags>0) {
444                // Temp storage to gather node IDs
445                int *Tags_keys = TMPMEMALLOC(num_Tags, int);
446                char name_temp[4096];
447            int i;
448    
449            // Tags_keys
450                if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) )
451                  throw DataException("get_var(Tags_keys)");
452                if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
453                  TMPMEMFREE(Tags_keys);
454                  throw DataException("get(Tags_keys)");
455                }
456            for (i=0; i<num_Tags; i++) {
457                  // Retrieve tag name
458                  sprintf(name_temp, "Tags_name_%d", i);
459                  if (! (attr=dataFile.get_att(name_temp)) ) {
460                    sprintf(error_msg,"Error retrieving tag name from NetCDF file '%s'", fName);
461                    throw DataException(error_msg);
462                  }
463                  char *name = attr->as_string(0);
464                  delete attr;
465                  Finley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
466            }
467          }
468        }
469      
470            if (Finley_noError()) Finley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
471            TMPMEMFREE(first_DofComponent);
472            TMPMEMFREE(first_NodeComponent);
473    
474      } /* Finley_noError() after Finley_Mesh_alloc() */      } /* Finley_noError() after Finley_Mesh_alloc() */
     
 #if 0 /* Not yet...finish the above first */  
     if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);  
     if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);  
 #endif  
475    
476      checkFinleyError();      checkFinleyError();
477      temp=new MeshAdapter(mesh_p);      temp=new MeshAdapter(mesh_p);
# Line 200  namespace finley { Line 484  namespace finley {
484      TMPMEMFREE(fName);      TMPMEMFREE(fName);
485    
486      blocktimer_increment("LoadMesh()", blocktimer_start);      blocktimer_increment("LoadMesh()", blocktimer_start);
487      return temp;      return temp->getPtr();
488  #else  #else
489      throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");      throw DataException("Error - loadMesh: is not compiled with NetCDF. Please contact your installation manager.");
490  #endif /* USE_NETCDF */  #endif /* USE_NETCDF */
491    }    }
492    
493    AbstractContinuousDomain* readMesh(const std::string& fileName,    Domain_ptr readMesh(const std::string& fileName,
494                       int integrationOrder,                       int integrationOrder,
495                                       int reducedIntegrationOrder,                                       int reducedIntegrationOrder,
496                                       int optimize)                                       int optimize)
# Line 216  namespace finley { Line 500  namespace finley {
500      // to Finley_Mesh_read      // to Finley_Mesh_read
501      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
502      // Win32 refactor      // Win32 refactor
503      char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;      if( fileName.size() == 0 )
504        {
505           throw DataException("Null file name!");
506        }
507    
508        char *fName = TMPMEMALLOC(fileName.size()+1,char);
509        
510      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
511      double blocktimer_start = blocktimer_time();      double blocktimer_start = blocktimer_time();
512    
# Line 228  namespace finley { Line 518  namespace finley {
518      TMPMEMFREE(fName);      TMPMEMFREE(fName);
519            
520      blocktimer_increment("ReadMesh()", blocktimer_start);      blocktimer_increment("ReadMesh()", blocktimer_start);
521      return temp;      return temp->getPtr();
522    }    }
523    
524    AbstractContinuousDomain* readGmsh(const std::string& fileName,    Domain_ptr readGmsh(const std::string& fileName,
525                                       int numDim,                                       int numDim,
526                                       int integrationOrder,                                       int integrationOrder,
527                                       int reducedIntegrationOrder,                                       int reducedIntegrationOrder,
528                                       int optimize)                                       int optimize,
529                                 int useMacroElements)
530    {    {
531      //      //
532      // 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
533      // to Finley_Mesh_read      // to Finley_Mesh_read
534      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
535      // Win32 refactor      // Win32 refactor
536      char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;      if( fileName.size() == 0 )
537        {
538           throw DataException("Null file name!");
539        }
540    
541        char *fName = TMPMEMALLOC(fileName.size()+1,char);
542        
543      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
544      double blocktimer_start = blocktimer_time();      double blocktimer_start = blocktimer_time();
545    
546      fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));      fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE), (useMacroElements ? TRUE : FALSE));
547      checkFinleyError();      checkFinleyError();
548      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
549            
# Line 254  namespace finley { Line 551  namespace finley {
551      TMPMEMFREE(fName);      TMPMEMFREE(fName);
552            
553      blocktimer_increment("ReadGmsh()", blocktimer_start);      blocktimer_increment("ReadGmsh()", blocktimer_start);
554      return temp;      return temp->getPtr();
555    }    }
556    
557    AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,  /*  AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,*/
558      Domain_ptr brick(int n0,int n1,int n2,int order,
559              double l0,double l1,double l2,              double l0,double l1,double l2,
560              int periodic0,int periodic1,              int periodic0,int periodic1,
561              int periodic2,              int periodic2,
# Line 276  namespace finley { Line 574  namespace finley {
574      Finley_Mesh* fMesh=NULL;      Finley_Mesh* fMesh=NULL;
575    
576      if (order==1) {      if (order==1) {
577        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,             fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
578                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
579      }      } else if (order==2) {
580          else if (order==2) {             fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
581        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,                            useElementsOnFace,useFullElementOrder,FALSE, (optimize ? TRUE : FALSE)) ;
582                       useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;      } else if (order==-1) {
583               fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
584                              useElementsOnFace,useFullElementOrder,TRUE,(optimize ? TRUE : FALSE)) ;
585      } else {      } else {
586        stringstream temp;        stringstream temp;
587        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
# Line 291  namespace finley { Line 591  namespace finley {
591      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
592      checkFinleyError();      checkFinleyError();
593      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
594      return temp;      return temp->getPtr();
595    }    }
596    AbstractContinuousDomain*  rectangle(int n0,int n1,int order,  
597    /*  AbstractContinuousDomain*  rectangle(int n0,int n1,int order,*/
598      Domain_ptr  rectangle(int n0,int n1,int order,
599              double l0, double l1,              double l0, double l1,
600              int periodic0,int periodic1,              int periodic0,int periodic1,
601              int integrationOrder,              int integrationOrder,
# Line 308  namespace finley { Line 610  namespace finley {
610    
611      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
612      if (order==1) {      if (order==1) {
613        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,              fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
                     useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));  
     }  
     else if (order==2) {  
       fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,  
614                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
615      }      } else if (order==2) {
616      else {              fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
617                              useElementsOnFace,useFullElementOrder,FALSE,(optimize ? TRUE : FALSE));
618        } else if (order==-1) {
619                fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
620                              useElementsOnFace,useFullElementOrder,TRUE,(optimize ? TRUE : FALSE));
621        } else {
622        stringstream temp;        stringstream temp;
623        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
624        setFinleyError(VALUE_ERROR,temp.str().c_str());        setFinleyError(VALUE_ERROR,temp.str().c_str());
# Line 324  namespace finley { Line 627  namespace finley {
627      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
628      checkFinleyError();      checkFinleyError();
629      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
630      return temp;      return temp->getPtr();
631    }    }
632    
633    AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)    Domain_ptr meshMerge(const boost::python::list& meshList)
634    {    {
635      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
636      //      //
# Line 348  namespace finley { Line 651  namespace finley {
651      checkFinleyError();      checkFinleyError();
652      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
653    
654      return temp;      return temp->getPtr();
655    }    }
656    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,  
657      Domain_ptr  glueFaces(const boost::python::list& meshList,
658                                     double safety_factor,                                     double safety_factor,
659                             double tolerance,                             double tolerance,
660                                         int optimize)                                         int optimize)
# Line 358  namespace finley { Line 662  namespace finley {
662      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
663      //      //
664      // merge the meshes:      // merge the meshes:
665      AbstractContinuousDomain* merged_meshes=meshMerge(meshList);      Domain_ptr merged_meshes=meshMerge(meshList);
666    
667      //      //
668      // glue the faces:      // glue the faces:
669      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);      const MeshAdapter* merged_finley_meshes=dynamic_cast<const MeshAdapter*>(merged_meshes.get());
670      fMesh=merged_finley_meshes->getFinley_Mesh();      fMesh=merged_finley_meshes->getFinley_Mesh();
671      Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));      Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
672    
673      //      //
674      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
675      checkFinleyError();      checkFinleyError();
676      return merged_meshes;      return merged_meshes->getPtr();
677    }    }
678    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,    Domain_ptr  joinFaces(const boost::python::list& meshList,
679              double safety_factor,              double safety_factor,
680              double tolerance,              double tolerance,
681                          int optimize)                          int optimize)
# Line 378  namespace finley { Line 683  namespace finley {
683      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
684      //      //
685      // merge the meshes:      // merge the meshes:
686      AbstractContinuousDomain* merged_meshes=meshMerge(meshList);      Domain_ptr merged_meshes=meshMerge(meshList);
687      //      //
688      // join the faces:      // join the faces:
689      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes.get());
690      fMesh=merged_finley_meshes->getFinley_Mesh();      fMesh=merged_finley_meshes->getFinley_Mesh();
691      Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));      Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
692      //      //
693      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
694      checkFinleyError();      checkFinleyError();
695      return merged_meshes;      return merged_meshes->getPtr();
696    }    }
697    
698    // end of namespace    // end of namespace

Legend:
Removed from v.1345  
changed lines
  Added in v.2779

  ViewVC Help
Powered by ViewVC 1.1.26