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

Diff of /branches/more_shared_ptrs_from_1812/finley/src/CPPAdapter/MeshAdapterFactory.cpp

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

trunk/esys2/finley/src/CPPAdapter/MeshAdapterFactory.cpp revision 110 by jgs, Mon Feb 14 04:14:42 2005 UTC trunk/finley/src/CPPAdapter/MeshAdapterFactory.cpp revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 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-2008 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 "escript/blocktimer.h"
 #include "finley/finleyC/Mesh.h"  
 #include "finley/finleyC/RectangularMesh.h"  
25  }  }
 #include "finley/CPPAdapter/FinleyError.h"  
 #include "finley/CPPAdapter/MeshAdapterFactory.h"  
26    
 #include <iostream>  
 #include <sstream>  
27  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
28    
29    #include <sstream>
30    
31  using namespace std;  using namespace std;
32  using namespace escript;  using namespace escript;
33    
34  namespace finley {  namespace finley {
35    
36    #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      {
53    #ifdef USE_NETCDF
54        bool optimize=FALSE; // Don't optimize since this would cause problems with Data().dump()
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              mesh_p->Elements=Finley_ElementFile_alloc((ElementTypeId)Elements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
196              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
197              mesh_p->Elements->minColor=0;
198              mesh_p->Elements->maxColor=num_Elements-1;
199              if (num_Elements>0) {
200                     if (Finley_noError()) {
201                   // Elements_Id
202                       if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
203                         throw DataException("Error - loadMesh:: unable to read Elements_Id from netCDF file: " + *fName);
204                       if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) ) {
205                         TMPMEMFREE(mesh_p->Elements->Id);
206                         throw DataException("Error - loadMesh:: unable to recover Elements_Id from NetCDF file: " + *fName);
207                       }
208                   // Elements_Tag
209                       if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
210                         throw DataException("Error - loadMesh:: unable to read Elements_Tag from netCDF file: " + *fName);
211                       if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) ) {
212                         TMPMEMFREE(mesh_p->Elements->Tag);
213                         throw DataException("Error - loadMesh:: unable to recover Elements_Tag from NetCDF file: " + *fName);
214                       }
215                   // Elements_Owner
216                       if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
217                         throw DataException("Error - loadMesh:: unable to read Elements_Owner from netCDF file: " + *fName);
218                       if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) ) {
219                         TMPMEMFREE(mesh_p->Elements->Owner);
220                         throw DataException("Error - loadMesh:: unable to recover Elements_Owner from NetCDF file: " + *fName);
221                       }
222                   // Elements_Color
223                       if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
224                         throw DataException("Error - loadMesh:: unable to read Elements_Color from netCDF file: " + *fName);
225                       if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) ) {
226                         TMPMEMFREE(mesh_p->Elements->Color);
227                         throw DataException("Error - loadMesh:: unable to recover Elements_Color from NetCDF file: " + *fName);
228                       }
229                   // Elements_Nodes
230               int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
231                       if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
232                         TMPMEMFREE(mesh_p->Elements->Nodes);
233                         throw DataException("Error - loadMesh:: unable to read Elements_Nodes from netCDF file: " + *fName);
234                       }
235                       if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
236                         TMPMEMFREE(Elements_Nodes);
237                         throw DataException("Error - load:: unable to recover Elements_Nodes from netCDF file: " + *fName);
238                       }
239               // Copy temp array into mesh_p->Elements->Nodes
240               for (int i=0; i<num_Elements; i++) {
241                 for (int j=0; j<num_Elements_numNodes; j++) {
242                   mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
243                     = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
244                 }
245               }
246               TMPMEMFREE(Elements_Nodes);
247             }
248          } /* num_Elements>0 */
249        }
250    
251            /* get the face elements */
252            if (Finley_noError()) {
253              mesh_p->FaceElements=Finley_ElementFile_alloc((ElementTypeId)FaceElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
254              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
255              mesh_p->FaceElements->minColor=0;
256              mesh_p->FaceElements->maxColor=num_FaceElements-1;
257              if (num_FaceElements>0) {
258                     if (Finley_noError()) {
259                   // FaceElements_Id
260                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
261                         throw DataException("Error - loadMesh:: unable to read FaceElements_Id from netCDF file: " + *fName);
262                       if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) ) {
263                         TMPMEMFREE(mesh_p->FaceElements->Id);
264                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Id from NetCDF file: " + *fName);
265                       }
266                   // FaceElements_Tag
267                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
268                         throw DataException("Error - loadMesh:: unable to read FaceElements_Tag from netCDF file: " + *fName);
269                       if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) ) {
270                         TMPMEMFREE(mesh_p->FaceElements->Tag);
271                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Tag from NetCDF file: " + *fName);
272                       }
273                   // FaceElements_Owner
274                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
275                         throw DataException("Error - loadMesh:: unable to read FaceElements_Owner from netCDF file: " + *fName);
276                       if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) ) {
277                         TMPMEMFREE(mesh_p->FaceElements->Owner);
278                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Owner from NetCDF file: " + *fName);
279                       }
280                   // FaceElements_Color
281                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
282                         throw DataException("Error - loadMesh:: unable to read FaceElements_Color from netCDF file: " + *fName);
283                       if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) ) {
284                         TMPMEMFREE(mesh_p->FaceElements->Color);
285                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Color from NetCDF file: " + *fName);
286                       }
287                   // FaceElements_Nodes
288               int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
289                       if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
290                         TMPMEMFREE(mesh_p->FaceElements->Nodes);
291                         throw DataException("Error - loadMesh:: unable to read FaceElements_Nodes from netCDF file: " + *fName);
292                       }
293                       if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
294                         TMPMEMFREE(FaceElements_Nodes);
295                         throw DataException("Error - load:: unable to recover FaceElements_Nodes from netCDF file: " + *fName);
296                       }
297               // Copy temp array into mesh_p->FaceElements->Nodes
298               for (int i=0; i<num_FaceElements; i++) {
299                 for (int j=0; j<num_FaceElements_numNodes; j++) {
300                   mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]
301                     = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
302                 }
303               }
304               TMPMEMFREE(FaceElements_Nodes);
305             }
306          } /* num_FaceElements>0 */
307        }
308    
309            /* get the Contact elements */
310            if (Finley_noError()) {
311              mesh_p->ContactElements=Finley_ElementFile_alloc((ElementTypeId)ContactElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
312              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->ContactElements, num_ContactElements);
313              mesh_p->ContactElements->minColor=0;
314              mesh_p->ContactElements->maxColor=num_ContactElements-1;
315              if (num_ContactElements>0) {
316                     if (Finley_noError()) {
317                   // ContactElements_Id
318                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Id")) )
319                         throw DataException("Error - loadMesh:: unable to read ContactElements_Id from netCDF file: " + *fName);
320                       if (! nc_var_temp->get(&mesh_p->ContactElements->Id[0], num_ContactElements) ) {
321                         TMPMEMFREE(mesh_p->ContactElements->Id);
322                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Id from NetCDF file: " + *fName);
323                       }
324                   // ContactElements_Tag
325                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Tag")) )
326                         throw DataException("Error - loadMesh:: unable to read ContactElements_Tag from netCDF file: " + *fName);
327                       if (! nc_var_temp->get(&mesh_p->ContactElements->Tag[0], num_ContactElements) ) {
328                         TMPMEMFREE(mesh_p->ContactElements->Tag);
329                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Tag from NetCDF file: " + *fName);
330                       }
331                   // ContactElements_Owner
332                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Owner")) )
333                         throw DataException("Error - loadMesh:: unable to read ContactElements_Owner from netCDF file: " + *fName);
334                       if (! nc_var_temp->get(&mesh_p->ContactElements->Owner[0], num_ContactElements) ) {
335                         TMPMEMFREE(mesh_p->ContactElements->Owner);
336                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Owner from NetCDF file: " + *fName);
337                       }
338                   // ContactElements_Color
339                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Color")) )
340                         throw DataException("Error - loadMesh:: unable to read ContactElements_Color from netCDF file: " + *fName);
341                       if (! nc_var_temp->get(&mesh_p->ContactElements->Color[0], num_ContactElements) ) {
342                         TMPMEMFREE(mesh_p->ContactElements->Color);
343                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Color from NetCDF file: " + *fName);
344                       }
345                   // ContactElements_Nodes
346               int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
347                       if (!(nc_var_temp = dataFile.get_var("ContactElements_Nodes"))) {
348                         TMPMEMFREE(mesh_p->ContactElements->Nodes);
349                         throw DataException("Error - loadMesh:: unable to read ContactElements_Nodes from netCDF file: " + *fName);
350                       }
351                       if (! nc_var_temp->get(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes) ) {
352                         TMPMEMFREE(ContactElements_Nodes);
353                         throw DataException("Error - load:: unable to recover ContactElements_Nodes from netCDF file: " + *fName);
354                       }
355               // Copy temp array into mesh_p->ContactElements->Nodes
356               for (int i=0; i<num_ContactElements; i++) {
357                 for (int j=0; j<num_ContactElements_numNodes; j++) {
358                   mesh_p->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]
359                     = ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
360                 }
361               }
362               TMPMEMFREE(ContactElements_Nodes);
363             }
364          } /* num_ContactElements>0 */
365        }
366    
367            /* get the Points (nodal elements) */
368            if (Finley_noError()) {
369              mesh_p->Points=Finley_ElementFile_alloc((ElementTypeId)Points_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
370              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Points, num_Points);
371              mesh_p->Points->minColor=0;
372              mesh_p->Points->maxColor=num_Points-1;
373              if (num_Points>0) {
374                 if (Finley_noError()) {
375                   // Points_Id
376                       if (! ( nc_var_temp = dataFile.get_var("Points_Id")) )
377                         throw DataException("Error - loadMesh:: unable to read Points_Id from netCDF file: " + *fName);
378                       if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points) ) {
379                         TMPMEMFREE(mesh_p->Points->Id);
380                         throw DataException("Error - loadMesh:: unable to recover Points_Id from NetCDF file: " + *fName);
381                       }
382                   // Points_Tag
383                       if (! ( nc_var_temp = dataFile.get_var("Points_Tag")) )
384                         throw DataException("Error - loadMesh:: unable to read Points_Tag from netCDF file: " + *fName);
385                       if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points) ) {
386                         TMPMEMFREE(mesh_p->Points->Tag);
387                         throw DataException("Error - loadMesh:: unable to recover Points_Tag from NetCDF file: " + *fName);
388                       }
389                   // Points_Owner
390                       if (! ( nc_var_temp = dataFile.get_var("Points_Owner")) )
391                         throw DataException("Error - loadMesh:: unable to read Points_Owner from netCDF file: " + *fName);
392                       if (! nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points) ) {
393                         TMPMEMFREE(mesh_p->Points->Owner);
394                         throw DataException("Error - loadMesh:: unable to recover Points_Owner from NetCDF file: " + *fName);
395                       }
396                   // Points_Color
397                       if (! ( nc_var_temp = dataFile.get_var("Points_Color")) )
398                         throw DataException("Error - loadMesh:: unable to read Points_Color from netCDF file: " + *fName);
399                       if (! nc_var_temp->get(&mesh_p->Points->Color[0], num_Points) ) {
400                         TMPMEMFREE(mesh_p->Points->Color);
401                         throw DataException("Error - loadMesh:: unable to recover Points_Color from NetCDF file: " + *fName);
402                       }
403                   // Points_Nodes
404               int *Points_Nodes = TMPMEMALLOC(num_Points,int);
405                       if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
406                         TMPMEMFREE(mesh_p->Points->Nodes);
407                         throw DataException("Error - loadMesh:: unable to read Points_Nodes from netCDF file: " + *fName);
408                       }
409                       if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
410                         TMPMEMFREE(Points_Nodes);
411                         throw DataException("Error - load:: unable to recover Points_Nodes from netCDF file: " + *fName);
412                       }
413               // Copy temp array into mesh_p->Points->Nodes
414               for (int i=0; i<num_Points; i++) {
415                 mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
416               }
417               TMPMEMFREE(Points_Nodes);
418             }
419          } /* num_Points>0 */
420        }
421    
422            /* get the tags */
423            if (Finley_noError()) {
424              if (num_Tags>0) {
425                // Temp storage to gather node IDs
426                int *Tags_keys = TMPMEMALLOC(num_Tags, int);
427                char name_temp[4096];
428            int i;
429    
430            // Tags_keys
431                if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) )
432                  throw DataException("Error - loadMesh:: unable to read Tags_keys from netCDF file: " + *fName);
433                if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
434                  TMPMEMFREE(Tags_keys);
435                  throw DataException("Error - loadMesh:: unable to recover Tags_keys from NetCDF file: " + *fName);
436                }
437            for (i=0; i<num_Tags; i++) {
438                  // Retrieve tag name
439                  sprintf(name_temp, "Tags_name_%d", i);
440                  if (! (attr=dataFile.get_att(name_temp)) ) {
441                    sprintf(error_msg,"Error retrieving tag name from NetCDF file '%s'", fName);
442                    throw DataException(error_msg);
443                  }
444                  char *name = attr->as_string(0);
445                  delete attr;
446                  Finley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
447            }
448          }
449        }
450    
451        } /* Finley_noError() after Finley_Mesh_alloc() */
452      
453        if (Finley_noError()) Finley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
454        TMPMEMFREE(first_DofComponent);
455        TMPMEMFREE(first_NodeComponent);
456    
457        checkFinleyError();
458        temp=new MeshAdapter(mesh_p);
459    
460        if (! Finley_noError()) {
461          Finley_Mesh_free(mesh_p);
462        }
463    
464        /* win32 refactor */
465        TMPMEMFREE(fName);
466    
467        blocktimer_increment("LoadMesh()", blocktimer_start);
468        return temp;
469    #else
470        throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
471    #endif /* USE_NETCDF */
472      }
473    
474    AbstractContinuousDomain* readMesh(const std::string& fileName,    AbstractContinuousDomain* readMesh(const std::string& fileName,
475                       int integrationOrder)                       int integrationOrder,
476                                         int reducedIntegrationOrder,
477                                         int optimize)
478    {    {
479      //      //
480      // 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
481      // to Finley_Mesh_read      // to Finley_Mesh_read
482      char fName[fileName.size()+1];      Finley_Mesh* fMesh=0;
483        // Win32 refactor
484        if( fileName.size() == 0 )
485        {
486           throw DataException("Null file name!");
487        }
488    
489        char *fName = TMPMEMALLOC(fileName.size()+1,char);
490        
491      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
492      Finley_Mesh* fMesh=Finley_Mesh_read(fName,integrationOrder);      double blocktimer_start = blocktimer_time();
493    
494        fMesh=Finley_Mesh_read_MPI(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
495      checkFinleyError();      checkFinleyError();
496      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
497        
498        /* win32 refactor */
499        TMPMEMFREE(fName);
500        
501        blocktimer_increment("ReadMesh()", blocktimer_start);
502        return temp;
503      }
504    
505      AbstractContinuousDomain* readGmsh(const std::string& fileName,
506                                         int numDim,
507                                         int integrationOrder,
508                                         int reducedIntegrationOrder,
509                                         int optimize)
510      {
511        //
512        // create a copy of the filename to overcome the non-constness of call
513        // to Finley_Mesh_read
514        Finley_Mesh* fMesh=0;
515        // Win32 refactor
516        if( fileName.size() == 0 )
517        {
518           throw DataException("Null file name!");
519        }
520    
521        char *fName = TMPMEMALLOC(fileName.size()+1,char);
522        
523        strcpy(fName,fileName.c_str());
524        double blocktimer_start = blocktimer_time();
525    
526        fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
527        checkFinleyError();
528        AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
529        
530        /* win32 refactor */
531        TMPMEMFREE(fName);
532        
533        blocktimer_increment("ReadGmsh()", blocktimer_start);
534      return temp;      return temp;
535    }    }
536    
# Line 49  namespace finley { Line 539  namespace finley {
539              int periodic0,int periodic1,              int periodic0,int periodic1,
540              int periodic2,              int periodic2,
541              int integrationOrder,              int integrationOrder,
542              int useElementsOnFace)                      int reducedIntegrationOrder,
543                int useElementsOnFace,
544                int useFullElementOrder,
545                        int optimize)
546    {    {
 //     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;  
         
547      int numElements[]={n0,n1,n2};      int numElements[]={n0,n1,n2};
548      double length[]={l0,l1,l2};      double length[]={l0,l1,l2};
549      int periodic[]={periodic0, periodic1, periodic2};      int periodic[]={periodic0, periodic1, periodic2};
550    
551      //      //
552      // linearInterpolation      // linearInterpolation
553      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=NULL;
554    
555      if (order==1) {      if (order==1) {
556        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
557                      useElementsOnFace) ;                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
558      } else if (order==2) {      }
559        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,          else if (order==2) {
560                       useElementsOnFace) ;        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
561                         useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
562      } else {      } else {
563        stringstream temp;        stringstream temp;
564        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
# Line 88  namespace finley { Line 574  namespace finley {
574              double l0, double l1,              double l0, double l1,
575              int periodic0,int periodic1,              int periodic0,int periodic1,
576              int integrationOrder,              int integrationOrder,
577              int useElementsOnFace)                          int reducedIntegrationOrder,
578                int useElementsOnFace,
579                    int useFullElementOrder,
580                            int optimize)
581    {    {
582      int numElements[]={n0,n1};      int numElements[]={n0,n1};
583      double length[]={l0,l1};      double length[]={l0,l1};
584      int periodic[]={periodic0, periodic1};      int periodic[]={periodic0, periodic1};
585    
586      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=0;
587      if (order==1) {      if (order==1) {
588        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
589                      useElementsOnFace);                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
     } else if (order==2) {  
       fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,  
                     useElementsOnFace);  
     } else {  
       stringstream temp;  
       temp << "Illegal interpolation order: " << order;  
       setFinleyError(VALUE_ERROR,temp.str().c_str());  
590      }      }
591      //      else if (order==2) {
592      // Convert any finley errors into a C++ exception        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
593      checkFinleyError();                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
594      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      }
595      return temp;      else {
   }  
   AbstractContinuousDomain*  interval(int n0,int order,double l0,int periodic0,  
                int integrationOrder,  
                int useElementsOnFace)  
   {  
     int numElements[]={n0};  
     double length[]={l0};  
     int periodic[]={periodic0};  
     Finley_Mesh* fMesh;  
     if (order==1) {  
       fMesh=Finley_RectangularMesh_Line2(numElements, length,periodic,integrationOrder,  
                      useElementsOnFace);  
     } else if (order==2) {  
       fMesh=Finley_RectangularMesh_Line3(numElements,length,periodic,integrationOrder,  
                      useElementsOnFace);  
     } else {  
596        stringstream temp;        stringstream temp;
597        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
598        setFinleyError(VALUE_ERROR,temp.str().c_str());        setFinleyError(VALUE_ERROR,temp.str().c_str());
# Line 137  namespace finley { Line 603  namespace finley {
603      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
604      return temp;      return temp;
605    }    }
606    
607    AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)    AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
608    {    {
609      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
610      //      //
611      // extract the meshes from meshList      // extract the meshes from meshList
612      int numMsh=boost::python::extract<int>(meshList.attr("__len__")());      int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
613      Finley_Mesh* mshes[numMsh];      Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
614      for (int i=0;i<numMsh;++i) {      for (int i=0;i<numMsh;++i) {
615           AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);           AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
616           const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);           const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
# Line 152  namespace finley { Line 619  namespace finley {
619      //      //
620      // merge the meshes:      // merge the meshes:
621      fMesh=Finley_Mesh_merge(numMsh,mshes);      fMesh=Finley_Mesh_merge(numMsh,mshes);
622          TMPMEMFREE(mshes);
623      //      //
624      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
625      checkFinleyError();      checkFinleyError();
626      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
627    
628      return temp;      return temp;
629    }    }
630    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,
631              double safety_factor,                                     double safety_factor,
632              double tolerance)                             double tolerance,
633                                           int optimize)
634    {    {
635      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
636      //      //
# Line 170  namespace finley { Line 640  namespace finley {
640      // glue the faces:      // glue the faces:
641      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
642      fMesh=merged_finley_meshes->getFinley_Mesh();      fMesh=merged_finley_meshes->getFinley_Mesh();
643      Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance);      Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
644    
645      //      //
646      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
647      checkFinleyError();      checkFinleyError();
# Line 178  namespace finley { Line 649  namespace finley {
649    }    }
650    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,
651              double safety_factor,              double safety_factor,
652              double tolerance)              double tolerance,
653                            int optimize)
654    {    {
655      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
656      //      //
# Line 188  namespace finley { Line 660  namespace finley {
660      // join the faces:      // join the faces:
661      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
662      fMesh=merged_finley_meshes->getFinley_Mesh();      fMesh=merged_finley_meshes->getFinley_Mesh();
663      Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance);      Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
664      //      //
665      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
666      checkFinleyError();      checkFinleyError();
667      return merged_meshes;      return merged_meshes;
668    }    }
669    
670  }  // end of namespace    // end of namespace
671    
672    }

Legend:
Removed from v.110  
changed lines
  Added in v.1811

  ViewVC Help
Powered by ViewVC 1.1.26