/[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 2748 by gross, Tue Nov 17 07:32:59 2009 UTC
# Line 1  Line 1 
 /* $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.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
1    
2    /*******************************************************
3    *
4    * Copyright (c) 2003-2009 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
14    
15    #ifdef PASO_MPI
16    #include <mpi.h>
17    #endif
18    #ifdef USE_NETCDF
19    #include <netcdfcpp.h>
20    #endif
21    #include "MeshAdapterFactory.h"
22    #include "FinleyError.h"
23  extern "C" {  extern "C" {
24  #include "finley/finleyC/Finley.h"  #include "esysUtils/blocktimer.h"
 #include "finley/finleyC/Mesh.h"  
 #include "finley/finleyC/RectangularMesh.h"  
25  }  }
 #include "finley/CPPAdapter/FinleyError.h"  
 #include "finley/CPPAdapter/MeshAdapterFactory.h"  
26    
27    #include <boost/python/extract.hpp>
28    
   
 #include <iostream>  
29  #include <sstream>  #include <sstream>
30    
31  using namespace std;  using namespace std;
# Line 31  using namespace escript; Line 33  using namespace escript;
33    
34  namespace finley {  namespace finley {
35    
36    AbstractContinuousDomain* readMesh(const std::string& fileName,  #ifdef USE_NETCDF
37                       int integrationOrder)    // 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        Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
56        AbstractContinuousDomain* temp;
57        Finley_Mesh *mesh_p=NULL;
58        char error_msg[LenErrorMsg_MAX];
59    
60        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        /* 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              TMPMEMFREE(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              TMPMEMFREE(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              TMPMEMFREE(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              TMPMEMFREE(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              TMPMEMFREE(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              TMPMEMFREE(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              TMPMEMFREE(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              TMPMEMFREE(mesh_p->Nodes->Coordinates);
175              throw DataException("Error - load:: unable to recover Nodes_Coordinates from netCDF file: " + *fName);
176            }
177        // Nodes_DofDistribution
178        first_DofComponent = 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_DofComponent[0], mpi_size+1) ) {
182              throw DataException("Error - loadMesh:: unable to recover Nodes_DofDistribution from NetCDF file: " + *fName);
183            }
184    
185        // Nodes_NodeDistribution
186        first_NodeComponent = TMPMEMALLOC(mpi_size+1,index_t);
187            if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) )
188              throw DataException("Error - loadMesh:: unable to read Nodes_NodeDistribution from netCDF file: " + *fName);
189            if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
190              throw DataException("Error - loadMesh:: unable to recover Nodes_NodeDistribution from NetCDF file: " + *fName);
191            }
192    
193            /* read elements */
194            if (Finley_noError()) {
195              Finley_ReferenceElementSet  *refElements= Finley_ReferenceElementSet_alloc((ElementTypeId)Elements_TypeId,mesh_p->order, mesh_p->reduced_order);
196              if (Finley_noError())  {
197                  mesh_p->Elements=Finley_ElementFile_alloc(refElements, mpi_info);
198              }
199              Finley_ReferenceElementSet_dealloc(refElements);
200              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
201              if (Finley_noError()) {
202                  mesh_p->Elements->minColor=0;
203                  mesh_p->Elements->maxColor=num_Elements-1;
204                  if (num_Elements>0) {
205                   // Elements_Id
206                       if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
207                         throw DataException("Error - loadMesh:: unable to read Elements_Id from netCDF file: " + *fName);
208                       if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) ) {
209                         TMPMEMFREE(mesh_p->Elements->Id);
210                         throw DataException("Error - loadMesh:: unable to recover Elements_Id from NetCDF file: " + *fName);
211                       }
212                   // Elements_Tag
213                       if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
214                         throw DataException("Error - loadMesh:: unable to read Elements_Tag from netCDF file: " + *fName);
215                       if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) ) {
216                         TMPMEMFREE(mesh_p->Elements->Tag);
217                         throw DataException("Error - loadMesh:: unable to recover Elements_Tag from NetCDF file: " + *fName);
218                       }
219                   // Elements_Owner
220                       if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
221                         throw DataException("Error - loadMesh:: unable to read Elements_Owner from netCDF file: " + *fName);
222                       if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) ) {
223                         TMPMEMFREE(mesh_p->Elements->Owner);
224                         throw DataException("Error - loadMesh:: unable to recover Elements_Owner from NetCDF file: " + *fName);
225                       }
226                   // Elements_Color
227                       if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
228                         throw DataException("Error - loadMesh:: unable to read Elements_Color from netCDF file: " + *fName);
229                       if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) ) {
230                         TMPMEMFREE(mesh_p->Elements->Color);
231                         throw DataException("Error - loadMesh:: unable to recover Elements_Color from NetCDF file: " + *fName);
232                       }
233                   // Elements_Nodes
234                        int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
235                       if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
236                         TMPMEMFREE(mesh_p->Elements->Nodes);
237                         throw DataException("Error - loadMesh:: unable to read Elements_Nodes from netCDF file: " + *fName);
238                       }
239                       if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
240                         TMPMEMFREE(Elements_Nodes);
241                         throw DataException("Error - load:: unable to recover Elements_Nodes from netCDF file: " + *fName);
242                       }
243                       // Copy temp array into mesh_p->Elements->Nodes
244                       for (int i=0; i<num_Elements; i++) {
245                           for (int j=0; j<num_Elements_numNodes; j++) {
246                               mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
247                               = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
248                           }
249                       }
250                       TMPMEMFREE(Elements_Nodes);
251                    
252                  } /* num_Elements>0 */
253              }
254        }
255    
256            /* get the face elements */
257            if (Finley_noError()) {
258              Finley_ReferenceElementSet *refFaceElements=  Finley_ReferenceElementSet_alloc((ElementTypeId)FaceElements_TypeId ,mesh_p->order, mesh_p->reduced_order);
259              if (Finley_noError())  {
260                  mesh_p->FaceElements=Finley_ElementFile_alloc(refFaceElements, mpi_info);
261              }
262              Finley_ReferenceElementSet_dealloc(refFaceElements);  
263              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
264              if (Finley_noError()) {
265                  mesh_p->FaceElements->minColor=0;
266                  mesh_p->FaceElements->maxColor=num_FaceElements-1;
267                  if (num_FaceElements>0) {
268                   // FaceElements_Id
269                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
270                         throw DataException("Error - loadMesh:: unable to read FaceElements_Id from netCDF file: " + *fName);
271                       if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) ) {
272                         TMPMEMFREE(mesh_p->FaceElements->Id);
273                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Id from NetCDF file: " + *fName);
274                       }
275                   // FaceElements_Tag
276                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
277                         throw DataException("Error - loadMesh:: unable to read FaceElements_Tag from netCDF file: " + *fName);
278                       if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) ) {
279                         TMPMEMFREE(mesh_p->FaceElements->Tag);
280                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Tag from NetCDF file: " + *fName);
281                       }
282                   // FaceElements_Owner
283                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
284                         throw DataException("Error - loadMesh:: unable to read FaceElements_Owner from netCDF file: " + *fName);
285                       if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) ) {
286                         TMPMEMFREE(mesh_p->FaceElements->Owner);
287                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Owner from NetCDF file: " + *fName);
288                       }
289                   // FaceElements_Color
290                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
291                         throw DataException("Error - loadMesh:: unable to read FaceElements_Color from netCDF file: " + *fName);
292                       if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) ) {
293                         TMPMEMFREE(mesh_p->FaceElements->Color);
294                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Color from NetCDF file: " + *fName);
295                       }
296                   // FaceElements_Nodes
297                        int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
298                       if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
299                         TMPMEMFREE(mesh_p->FaceElements->Nodes);
300                         throw DataException("Error - loadMesh:: unable to read FaceElements_Nodes from netCDF file: " + *fName);
301                       }
302                       if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
303                         TMPMEMFREE(FaceElements_Nodes);
304                         throw DataException("Error - load:: unable to recover FaceElements_Nodes from netCDF file: " + *fName);
305                       }
306                     // Copy temp array into mesh_p->FaceElements->Nodes
307                        for (int i=0; i<num_FaceElements; i++) {
308                            for (int j=0; j<num_FaceElements_numNodes; j++) {
309                                mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
310                            }
311                        }
312                        TMPMEMFREE(FaceElements_Nodes);
313                  } /* num_FaceElements>0 */
314              }
315        }
316    
317            /* get the Contact elements */
318            if (Finley_noError()) {
319              Finley_ReferenceElementSet *refContactElements=   Finley_ReferenceElementSet_alloc((ElementTypeId)ContactElements_TypeId,mesh_p->order, mesh_p->reduced_order);
320              if (Finley_noError())  {
321                  mesh_p->ContactElements=Finley_ElementFile_alloc(refContactElements, mpi_info);
322              }
323              Finley_ReferenceElementSet_dealloc(refContactElements);  
324              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->ContactElements, num_ContactElements);
325              if (Finley_noError()) {
326                  mesh_p->ContactElements->minColor=0;
327                  mesh_p->ContactElements->maxColor=num_ContactElements-1;
328                  if (num_ContactElements>0) {
329                   // ContactElements_Id
330                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Id")) )
331                         throw DataException("Error - loadMesh:: unable to read ContactElements_Id from netCDF file: " + *fName);
332                       if (! nc_var_temp->get(&mesh_p->ContactElements->Id[0], num_ContactElements) ) {
333                         TMPMEMFREE(mesh_p->ContactElements->Id);
334                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Id from NetCDF file: " + *fName);
335                       }
336                   // ContactElements_Tag
337                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Tag")) )
338                         throw DataException("Error - loadMesh:: unable to read ContactElements_Tag from netCDF file: " + *fName);
339                       if (! nc_var_temp->get(&mesh_p->ContactElements->Tag[0], num_ContactElements) ) {
340                         TMPMEMFREE(mesh_p->ContactElements->Tag);
341                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Tag from NetCDF file: " + *fName);
342                       }
343                   // ContactElements_Owner
344                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Owner")) )
345                         throw DataException("Error - loadMesh:: unable to read ContactElements_Owner from netCDF file: " + *fName);
346                       if (! nc_var_temp->get(&mesh_p->ContactElements->Owner[0], num_ContactElements) ) {
347                         TMPMEMFREE(mesh_p->ContactElements->Owner);
348                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Owner from NetCDF file: " + *fName);
349                       }
350                   // ContactElements_Color
351                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Color")) )
352                         throw DataException("Error - loadMesh:: unable to read ContactElements_Color from netCDF file: " + *fName);
353                       if (! nc_var_temp->get(&mesh_p->ContactElements->Color[0], num_ContactElements) ) {
354                         TMPMEMFREE(mesh_p->ContactElements->Color);
355                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Color from NetCDF file: " + *fName);
356                       }
357                   // ContactElements_Nodes
358                        int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
359                       if (!(nc_var_temp = dataFile.get_var("ContactElements_Nodes"))) {
360                         TMPMEMFREE(mesh_p->ContactElements->Nodes);
361                         throw DataException("Error - loadMesh:: unable to read ContactElements_Nodes from netCDF file: " + *fName);
362                       }
363                       if (! nc_var_temp->get(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes) ) {
364                         TMPMEMFREE(ContactElements_Nodes);
365                         throw DataException("Error - load:: unable to recover ContactElements_Nodes from netCDF file: " + *fName);
366                       }
367                       // Copy temp array into mesh_p->ContactElements->Nodes
368                        for (int i=0; i<num_ContactElements; i++) {
369                            for (int j=0; j<num_ContactElements_numNodes; j++) {
370                                mesh_p->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]= ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
371                            }
372                        }
373                        TMPMEMFREE(ContactElements_Nodes);
374                     } /* num_ContactElements>0 */
375                  }
376              
377        }
378    
379            /* get the Points (nodal elements) */
380            if (Finley_noError()) {
381              Finley_ReferenceElementSet *refPoints=    Finley_ReferenceElementSet_alloc((ElementTypeId)Points_TypeId,mesh_p->order, mesh_p->reduced_order);
382              if (Finley_noError())  {
383                  mesh_p->Points=Finley_ElementFile_alloc(refPoints, mpi_info);
384              }
385              Finley_ReferenceElementSet_dealloc(refPoints);
386              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Points, num_Points);
387              if (Finley_noError()) {
388                  mesh_p->Points->minColor=0;
389                  mesh_p->Points->maxColor=num_Points-1;
390                  if (num_Points>0) {
391                   // Points_Id
392                       if (! ( nc_var_temp = dataFile.get_var("Points_Id")) )
393                         throw DataException("Error - loadMesh:: unable to read Points_Id from netCDF file: " + *fName);
394                       if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points) ) {
395                         TMPMEMFREE(mesh_p->Points->Id);
396                         throw DataException("Error - loadMesh:: unable to recover Points_Id from NetCDF file: " + *fName);
397                       }
398                   // Points_Tag
399                       if (! ( nc_var_temp = dataFile.get_var("Points_Tag")) )
400                         throw DataException("Error - loadMesh:: unable to read Points_Tag from netCDF file: " + *fName);
401                       if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points) ) {
402                         TMPMEMFREE(mesh_p->Points->Tag);
403                         throw DataException("Error - loadMesh:: unable to recover Points_Tag from NetCDF file: " + *fName);
404                       }
405                   // Points_Owner
406                       if (! ( nc_var_temp = dataFile.get_var("Points_Owner")) )
407                         throw DataException("Error - loadMesh:: unable to read Points_Owner from netCDF file: " + *fName);
408                       if (! nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points) ) {
409                         TMPMEMFREE(mesh_p->Points->Owner);
410                         throw DataException("Error - loadMesh:: unable to recover Points_Owner from NetCDF file: " + *fName);
411                       }
412                   // Points_Color
413                       if (! ( nc_var_temp = dataFile.get_var("Points_Color")) )
414                         throw DataException("Error - loadMesh:: unable to read Points_Color from netCDF file: " + *fName);
415                       if (! nc_var_temp->get(&mesh_p->Points->Color[0], num_Points) ) {
416                         TMPMEMFREE(mesh_p->Points->Color);
417                         throw DataException("Error - loadMesh:: unable to recover Points_Color from NetCDF file: " + *fName);
418                       }
419                   // Points_Nodes
420               int *Points_Nodes = TMPMEMALLOC(num_Points,int);
421                       if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
422                         TMPMEMFREE(mesh_p->Points->Nodes);
423                         throw DataException("Error - loadMesh:: unable to read Points_Nodes from netCDF file: " + *fName);
424                       }
425                       if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
426                         TMPMEMFREE(Points_Nodes);
427                         throw DataException("Error - load:: unable to recover Points_Nodes from netCDF file: " + *fName);
428                       }
429               // Copy temp array into mesh_p->Points->Nodes
430               for (int i=0; i<num_Points; i++) {
431                 mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
432               }
433               TMPMEMFREE(Points_Nodes);
434            
435                  } /* num_Points>0 */
436              }
437            }
438    
439            /* get the tags */
440            if (Finley_noError()) {
441              if (num_Tags>0) {
442                // Temp storage to gather node IDs
443                int *Tags_keys = TMPMEMALLOC(num_Tags, int);
444                char name_temp[4096];
445            int i;
446    
447            // Tags_keys
448                if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) )
449                  throw DataException("Error - loadMesh:: unable to read Tags_keys from netCDF file: " + *fName);
450                if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
451                  TMPMEMFREE(Tags_keys);
452                  throw DataException("Error - loadMesh:: unable to recover Tags_keys from NetCDF file: " + *fName);
453                }
454            for (i=0; i<num_Tags; i++) {
455                  // Retrieve tag name
456                  sprintf(name_temp, "Tags_name_%d", i);
457                  if (! (attr=dataFile.get_att(name_temp)) ) {
458                    sprintf(error_msg,"Error retrieving tag name from NetCDF file '%s'", fName);
459                    throw DataException(error_msg);
460                  }
461                  char *name = attr->as_string(0);
462                  delete attr;
463                  Finley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
464            }
465          }
466        }
467      
468            if (Finley_noError()) Finley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
469            TMPMEMFREE(first_DofComponent);
470            TMPMEMFREE(first_NodeComponent);
471    
472        } /* Finley_noError() after Finley_Mesh_alloc() */
473    
474        checkFinleyError();
475        temp=new MeshAdapter(mesh_p);
476    
477        if (! Finley_noError()) {
478          Finley_Mesh_free(mesh_p);
479        }
480    
481        /* win32 refactor */
482        TMPMEMFREE(fName);
483    
484        blocktimer_increment("LoadMesh()", blocktimer_start);
485        return temp->getPtr();
486    #else
487        throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
488    #endif /* USE_NETCDF */
489      }
490    
491      Domain_ptr readMesh(const std::string& fileName,
492                         int integrationOrder,
493                                         int reducedIntegrationOrder,
494                                         int optimize)
495      {
496        //
497        // create a copy of the filename to overcome the non-constness of call
498        // to Finley_Mesh_read
499        Finley_Mesh* fMesh=0;
500        // Win32 refactor
501        if( fileName.size() == 0 )
502        {
503           throw DataException("Null file name!");
504        }
505    
506        char *fName = TMPMEMALLOC(fileName.size()+1,char);
507        
508        strcpy(fName,fileName.c_str());
509        double blocktimer_start = blocktimer_time();
510    
511        fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
512        checkFinleyError();
513        AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
514        
515        /* win32 refactor */
516        TMPMEMFREE(fName);
517        
518        blocktimer_increment("ReadMesh()", blocktimer_start);
519        return temp->getPtr();
520      }
521    
522      Domain_ptr readGmsh(const std::string& fileName,
523                                         int numDim,
524                                         int integrationOrder,
525                                         int reducedIntegrationOrder,
526                                         int optimize,
527                                 int useMacroElements)
528    {    {
529      //      //
530      // 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
531      // to Finley_Mesh_read      // to Finley_Mesh_read
532      char fName[fileName.size()+1];      Finley_Mesh* fMesh=0;
533        // Win32 refactor
534        if( fileName.size() == 0 )
535        {
536           throw DataException("Null file name!");
537        }
538    
539        char *fName = TMPMEMALLOC(fileName.size()+1,char);
540        
541      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
542      Finley_Mesh* fMesh=Finley_Mesh_read(fName,integrationOrder);      double blocktimer_start = blocktimer_time();
543    
544        fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE), (useMacroElements ? TRUE : FALSE));
545      checkFinleyError();      checkFinleyError();
546      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
547      return temp;      
548        /* win32 refactor */
549        TMPMEMFREE(fName);
550        
551        blocktimer_increment("ReadGmsh()", blocktimer_start);
552        return temp->getPtr();
553    }    }
554    
555    AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,  /*  AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,*/
556      Domain_ptr brick(int n0,int n1,int n2,int order,
557              double l0,double l1,double l2,              double l0,double l1,double l2,
558              int periodic0,int periodic1,              int periodic0,int periodic1,
559              int periodic2,              int periodic2,
560              int integrationOrder,              int integrationOrder,
561              int useElementsOnFace)                      int reducedIntegrationOrder,
562                int useElementsOnFace,
563                int useFullElementOrder,
564                        int optimize)
565    {    {
 //     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;  
         
566      int numElements[]={n0,n1,n2};      int numElements[]={n0,n1,n2};
567      double length[]={l0,l1,l2};      double length[]={l0,l1,l2};
568      int periodic[]={periodic0, periodic1, periodic2};      int periodic[]={periodic0, periodic1, periodic2};
569    
570      //      //
571      // linearInterpolation      // linearInterpolation
572      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=NULL;
573    
574      if (order==1) {      if (order==1) {
575        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,             fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
576                      useElementsOnFace) ;                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
577      } else if (order==2) {      } else if (order==2) {
578        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,             fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
579                       useElementsOnFace) ;                            useElementsOnFace,useFullElementOrder,FALSE, (optimize ? TRUE : FALSE)) ;
580        } else if (order==-1) {
581               fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
582                              useElementsOnFace,useFullElementOrder,TRUE,(optimize ? TRUE : FALSE)) ;
583      } else {      } else {
584        stringstream temp;        stringstream temp;
585        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
# Line 83  namespace finley { Line 589  namespace finley {
589      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
590      checkFinleyError();      checkFinleyError();
591      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
592      return temp;      return temp->getPtr();
593    }    }
594    AbstractContinuousDomain*  rectangle(int n0,int n1,int order,  
595    /*  AbstractContinuousDomain*  rectangle(int n0,int n1,int order,*/
596      Domain_ptr  rectangle(int n0,int n1,int order,
597              double l0, double l1,              double l0, double l1,
598              int periodic0,int periodic1,              int periodic0,int periodic1,
599              int integrationOrder,              int integrationOrder,
600              int useElementsOnFace)                          int reducedIntegrationOrder,
601                int useElementsOnFace,
602                    int useFullElementOrder,
603                            int optimize)
604    {    {
605      int numElements[]={n0,n1};      int numElements[]={n0,n1};
606      double length[]={l0,l1};      double length[]={l0,l1};
607      int periodic[]={periodic0, periodic1};      int periodic[]={periodic0, periodic1};
608    
609      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=0;
610      if (order==1) {      if (order==1) {
611        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,              fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
612                      useElementsOnFace);                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
613      } else if (order==2) {      } else if (order==2) {
614        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,              fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
615                      useElementsOnFace);                            useElementsOnFace,useFullElementOrder,FALSE,(optimize ? TRUE : FALSE));
616        } else if (order==-1) {
617                fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
618                              useElementsOnFace,useFullElementOrder,TRUE,(optimize ? TRUE : FALSE));
619      } else {      } else {
620        stringstream temp;        stringstream temp;
621        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
# Line 111  namespace finley { Line 625  namespace finley {
625      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
626      checkFinleyError();      checkFinleyError();
627      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
628      return temp;      return temp->getPtr();
629    }    }
630    AbstractContinuousDomain*  interval(int n0,int order,double l0,int periodic0,  
631                 int integrationOrder,    Domain_ptr 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;  
652    }      return temp->getPtr();
   AbstractContinuousDomain*  meshMerge(const boost::python::list& meshList)  
   {  
     AbstractContinuousDomain* temp=new MeshAdapter(0);  
     return temp;  
653    }    }
654    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,  
655              double safetyFactor,    Domain_ptr  glueFaces(const boost::python::list& meshList,
656              double tolerance)                                     double safety_factor,
657                               double tolerance,
658                                           int optimize)
659    {    {
660      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
661      return temp;      //
662        // merge the meshes:
663        Domain_ptr merged_meshes=meshMerge(meshList);
664    
665        //
666        // glue the faces:
667        const MeshAdapter* merged_finley_meshes=dynamic_cast<const MeshAdapter*>(merged_meshes.get());
668        fMesh=merged_finley_meshes->getFinley_Mesh();
669        Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
670    
671        //
672        // Convert any finley errors into a C++ exception
673        checkFinleyError();
674        return merged_meshes->getPtr();
675    }    }
676    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,    Domain_ptr  joinFaces(const boost::python::list& meshList,
677              double safety_factor,              double safety_factor,
678              double tolerance)              double tolerance,
679                            int optimize)
680    {    {
681      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
682      return temp;      //
683        // merge the meshes:
684        Domain_ptr merged_meshes=meshMerge(meshList);
685        //
686        // join the faces:
687        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes.get());
688        fMesh=merged_finley_meshes->getFinley_Mesh();
689        Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
690        //
691        // Convert any finley errors into a C++ exception
692        checkFinleyError();
693        return merged_meshes->getPtr();
694    }    }
695    
696  }  // end of namespace    // end of namespace
697    
698    }

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

  ViewVC Help
Powered by ViewVC 1.1.26