/[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 1347 by ksteube, Fri Nov 16 05:37:07 2007 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));
     fMesh=Finley_Mesh_read(fName,integrationOrder);  
 #else  
     {  
       stringstream temp;  
       temp << "Unable to read meshes from file under MPI yet...";  
       setFinleyError(VALUE_ERROR,temp.str().c_str());  
     }  
 #endif  
479      checkFinleyError();      checkFinleyError();
480      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
481            
482      /* win32 refactor */      /* win32 refactor */
483      TMPMEMFREE(fName);      TMPMEMFREE(fName);
484            
485        blocktimer_increment("ReadMesh()", blocktimer_start);
486        return temp;
487      }
488    
489      AbstractContinuousDomain* readGmsh(const std::string& fileName,
490                                         int numDim,
491                                         int integrationOrder,
492                                         int reducedIntegrationOrder,
493                                         int optimize)
494      {
495        //
496        // create a copy of the filename to overcome the non-constness of call
497        // to Finley_Mesh_read
498        Finley_Mesh* fMesh=0;
499        // Win32 refactor
500        char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;
501        strcpy(fName,fileName.c_str());
502        double blocktimer_start = blocktimer_time();
503    
504        fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
505        checkFinleyError();
506        AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
507        
508        /* win32 refactor */
509        TMPMEMFREE(fName);
510        
511        blocktimer_increment("ReadGmsh()", blocktimer_start);
512      return temp;      return temp;
513    }    }
514    
# Line 59  namespace finley { Line 517  namespace finley {
517              int periodic0,int periodic1,              int periodic0,int periodic1,
518              int periodic2,              int periodic2,
519              int integrationOrder,              int integrationOrder,
520              int useElementsOnFace)                      int reducedIntegrationOrder,
521                int useElementsOnFace,
522                int useFullElementOrder,
523                        int optimize)
524    {    {
 //     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;  
         
525      int numElements[]={n0,n1,n2};      int numElements[]={n0,n1,n2};
526      double length[]={l0,l1,l2};      double length[]={l0,l1,l2};
527      int periodic[]={periodic0, periodic1, periodic2};      int periodic[]={periodic0, periodic1, periodic2};
# Line 79  namespace finley { Line 531  namespace finley {
531      Finley_Mesh* fMesh=NULL;      Finley_Mesh* fMesh=NULL;
532    
533      if (order==1) {      if (order==1) {
534        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
535                      useElementsOnFace) ;                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
536      }      }
 #ifndef PASO_MPI  
537          else if (order==2) {          else if (order==2) {
538        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
539                       useElementsOnFace) ;                       useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
540      } else {      } else {
541        stringstream temp;        stringstream temp;
542        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
543        setFinleyError(VALUE_ERROR,temp.str().c_str());        setFinleyError(VALUE_ERROR,temp.str().c_str());
544      }      }
 #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  
545      //      //
546      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
547      checkFinleyError();      checkFinleyError();
# Line 108  namespace finley { Line 552  namespace finley {
552              double l0, double l1,              double l0, double l1,
553              int periodic0,int periodic1,              int periodic0,int periodic1,
554              int integrationOrder,              int integrationOrder,
555              int useElementsOnFace)                          int reducedIntegrationOrder,
556                int useElementsOnFace,
557                    int useFullElementOrder,
558                            int optimize)
559    {    {
560      int numElements[]={n0,n1};      int numElements[]={n0,n1};
561      double length[]={l0,l1};      double length[]={l0,l1};
# Line 116  namespace finley { Line 563  namespace finley {
563    
564      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
565      if (order==1) {      if (order==1) {
566        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
567                      useElementsOnFace);                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
568      }      }
 #ifndef PASO_MPI  
569      else if (order==2) {      else if (order==2) {
570        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
571                      useElementsOnFace);                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
572      }      }
 #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*  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  
573      else {      else {
574        stringstream temp;        stringstream temp;
575        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
# Line 165  namespace finley { Line 581  namespace finley {
581      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
582      return temp;      return temp;
583    }    }
584    
585    AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)    AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
586    {    {
587      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
588      //      //
589      // extract the meshes from meshList      // extract the meshes from meshList
 #ifndef PASO_MPI  
590      int numMsh=boost::python::extract<int>(meshList.attr("__len__")());      int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
591      Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;      Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
592      for (int i=0;i<numMsh;++i) {      for (int i=0;i<numMsh;++i) {
# Line 181  namespace finley { Line 597  namespace finley {
597      //      //
598      // merge the meshes:      // merge the meshes:
599      fMesh=Finley_Mesh_merge(numMsh,mshes);      fMesh=Finley_Mesh_merge(numMsh,mshes);
600  #else        TMPMEMFREE(mshes);
     {  
       stringstream temp;  
       temp << "meshMerge() not available in MPI yet...";  
       setFinleyError(VALUE_ERROR,temp.str().c_str());  
     }  
 #endif  
601      //      //
602      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
603      checkFinleyError();      checkFinleyError();
604      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
     TMPMEMFREE(mshes);  
605    
606      return temp;      return temp;
607    }    }
608    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,
609              double safety_factor,                                     double safety_factor,
610              double tolerance)                             double tolerance,
611                                           int optimize)
612    {    {
613      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
 #ifndef PASO_MPI  
614      //      //
615      // merge the meshes:      // merge the meshes:
616      AbstractContinuousDomain* merged_meshes=meshMerge(meshList);      AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
# Line 209  namespace finley { Line 618  namespace finley {
618      // glue the faces:      // glue the faces:
619      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
620      fMesh=merged_finley_meshes->getFinley_Mesh();      fMesh=merged_finley_meshes->getFinley_Mesh();
621      Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance);      Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
622    
623      //      //
624      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
625      checkFinleyError();      checkFinleyError();
626      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  
   
627    }    }
628    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,
629              double safety_factor,              double safety_factor,
630              double tolerance)              double tolerance,
631                            int optimize)
632    {    {
633      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
634      //      //
635      // merge the meshes:      // merge the meshes:
 #ifndef PASO_MPI  
636      AbstractContinuousDomain* merged_meshes=meshMerge(meshList);      AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
637      //      //
638      // join the faces:      // join the faces:
639      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
640      fMesh=merged_finley_meshes->getFinley_Mesh();      fMesh=merged_finley_meshes->getFinley_Mesh();
641      Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance);      Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
642      //      //
643      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
644      checkFinleyError();      checkFinleyError();
645      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  
646    }    }
647    
648  }  // end of namespace    // end of namespace
649    
650    }

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

  ViewVC Help
Powered by ViewVC 1.1.26