/[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 1312 by ksteube, Mon Sep 24 06:18:44 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>
17  #endif  #endif
18    #ifdef USE_NETCDF
19    #include <netcdfcpp.h>
20    #endif
21  #include "MeshAdapterFactory.h"  #include "MeshAdapterFactory.h"
22  #include "FinleyError.h"  #include "FinleyError.h"
23    extern "C" {
24    #include "esysUtils/blocktimer.h"
25    }
26    
27  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
28    
# Line 28  using namespace escript; Line 33  using namespace escript;
33    
34  namespace finley {  namespace finley {
35    
36    AbstractContinuousDomain* loadMesh(const std::string& fileName)  #ifdef USE_NETCDF
37      // A convenience method to retrieve an integer attribute from a NetCDF file
38      int NetCDF_Get_Int_Attribute(NcFile *dataFile, char *fName, char *attr_name) {
39        NcAtt *attr;
40        char error_msg[LenErrorMsg_MAX];
41        if (! (attr=dataFile->get_att(attr_name)) ) {
42          sprintf(error_msg,"Error retrieving integer attribute '%s' from NetCDF file '%s'", attr_name, fName);
43          throw DataException(error_msg);
44        }
45        int temp = attr->as_int(0);
46        delete attr;
47        return(temp);
48      }
49    #endif
50    
51    //   AbstractContinuousDomain* loadMesh(const std::string& fileName)
52      Domain_ptr loadMesh(const std::string& fileName)
53    {    {
54      //  #ifdef USE_NETCDF
55      // create a copy of the filename to overcome the non-constness of call      Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
56      // to Finley_Mesh_read      AbstractContinuousDomain* temp;
57      Finley_Mesh* fMesh=0;      Finley_Mesh *mesh_p=NULL;
58      // Win32 refactor      char error_msg[LenErrorMsg_MAX];
59      char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;  
60      strcpy(fName,fileName.c_str());      char *fName = Paso_MPI_appendRankToFileName(fileName.c_str(),
61                                                    mpi_info->size,
62                                                    mpi_info->rank);
63    
64        double blocktimer_start = blocktimer_time();
65        Finley_resetError();
66        int *first_DofComponent, *first_NodeComponent;
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        string msgPrefix("Error in loadMesh: NetCDF operation failed - ");
121    
122        /* allocate mesh */
123        mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
124        if (Finley_noError()) {
125    
126            /* read nodes */
127            Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
128            // Nodes_Id
129            if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
130              throw DataException(msgPrefix+"get_var(Nodes_Id)");
131            if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {
132              TMPMEMFREE(mesh_p->Nodes->Id);
133              throw DataException("get(Nodes_Id)");
134            }
135            // Nodes_Tag
136            if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
137              throw DataException("get_var(Nodes_Tag)");
138            if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) ) {
139              TMPMEMFREE(mesh_p->Nodes->Tag);
140              throw DataException("get(Nodes_Tag)");
141            }
142            // Nodes_gDOF
143            if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
144              throw DataException("get_var(Nodes_gDOF)");
145            if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) ) {
146              TMPMEMFREE(mesh_p->Nodes->globalDegreesOfFreedom);
147              throw DataException("get(Nodes_gDOF)");
148            }
149            // Nodes_gNI
150            if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
151              throw DataException("get_var(Nodes_gNI)");
152            if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) ) {
153              TMPMEMFREE(mesh_p->Nodes->globalNodesIndex);
154              throw DataException("get(Nodes_gNI)");
155            }
156            // Nodes_grDfI
157            if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
158              throw DataException("get_var(Nodes_grDfI)");
159            if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) ) {
160              TMPMEMFREE(mesh_p->Nodes->globalReducedDOFIndex);
161              throw DataException("get(Nodes_grDfI)");
162            }
163            // Nodes_grNI
164            if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
165              throw DataException("get_var(Nodes_grNI)");
166            if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) ) {
167              TMPMEMFREE(mesh_p->Nodes->globalReducedNodesIndex);
168              throw DataException("get(Nodes_grNI)");
169            }
170            // Nodes_Coordinates
171            if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates"))) {
172              TMPMEMFREE(mesh_p->Nodes->Coordinates);
173              throw DataException("get_var(Nodes_Coordinates)");
174            }
175            if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) ) {
176              TMPMEMFREE(mesh_p->Nodes->Coordinates);
177              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    
195            /* read elements */
196            if (Finley_noError()) {
197              Finley_ReferenceElementSet  *refElements= Finley_ReferenceElementSet_alloc((ElementTypeId)Elements_TypeId,mesh_p->order, mesh_p->reduced_order);
198              if (Finley_noError())  {
199                  mesh_p->Elements=Finley_ElementFile_alloc(refElements, mpi_info);
200              }
201              Finley_ReferenceElementSet_dealloc(refElements);
202              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        }
257    
258            /* 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 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 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 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() */
475    
     fMesh=Finley_Mesh_load(fName);  
476      checkFinleyError();      checkFinleyError();
477      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      temp=new MeshAdapter(mesh_p);
478        
479        if (! Finley_noError()) {
480          Finley_Mesh_free(mesh_p);
481        }
482    
483      /* win32 refactor */      /* win32 refactor */
484      TMPMEMFREE(fName);      TMPMEMFREE(fName);
485        
486      return temp;      blocktimer_increment("LoadMesh()", blocktimer_start);
487        return temp->getPtr();
488    #else
489        throw DataException("Error - loadMesh: is not compiled with NetCDF. Please contact your installation manager.");
490    #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 58  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();
512    
513      fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));      fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
514      checkFinleyError();      checkFinleyError();
# Line 68  namespace finley { Line 517  namespace finley {
517      /* win32 refactor */      /* win32 refactor */
518      TMPMEMFREE(fName);      TMPMEMFREE(fName);
519            
520      return temp;      blocktimer_increment("ReadMesh()", blocktimer_start);
521        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();
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            
550      /* win32 refactor */      /* win32 refactor */
551      TMPMEMFREE(fName);      TMPMEMFREE(fName);
552            
553      return temp;      blocktimer_increment("ReadGmsh()", blocktimer_start);
554        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 114  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 129  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 146  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 162  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 186  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 196  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 216  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.1312  
changed lines
  Added in v.2779

  ViewVC Help
Powered by ViewVC 1.1.26