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

Legend:
Removed from v.757  
changed lines
  Added in v.1821

  ViewVC Help
Powered by ViewVC 1.1.26