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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26