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

Diff of /trunk/finley/src/CPPAdapter/MeshAdapterFactory.cpp

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

trunk/esys2/finley/src/CPPAdapter/MeshAdapterFactory.cpp revision 102 by jgs, Wed Dec 15 07:08:39 2004 UTC trunk/finley/src/CPPAdapter/MeshAdapterFactory.cpp revision 1360 by ksteube, Thu Dec 13 05:05:00 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"
23    #include "FinleyError.h"
24  extern "C" {  extern "C" {
25  #include "finley/finleyC/Finley.h"  #include "escript/blocktimer.h"
 #include "finley/finleyC/Mesh.h"  
 #include "finley/finleyC/RectangularMesh.h"  
26  }  }
 #include "finley/CPPAdapter/FinleyError.h"  
 #include "finley/CPPAdapter/MeshAdapterFactory.h"  
   
27    
28    #include <boost/python/extract.hpp>
29    
 #include <iostream>  
30  #include <sstream>  #include <sstream>
31    
32  using namespace std;  using namespace std;
# Line 31  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
471      // to Finley_Mesh_read      // to Finley_Mesh_read
472      char fName[fileName.size()+1];      Finley_Mesh* fMesh=0;
473        // Win32 refactor
474        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      Finley_Mesh* fMesh=Finley_Mesh_read(fName,integrationOrder);      double blocktimer_start = blocktimer_time();
477    
478        fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
479      checkFinleyError();      checkFinleyError();
480      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
481        
482        /* win32 refactor */
483        TMPMEMFREE(fName);
484        
485        blocktimer_increment("ReadMesh()", blocktimer_start);
486        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();
505        AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
506        
507        /* win32 refactor */
508        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 50  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};
553    
554      //      //
555      // linearInterpolation      // linearInterpolation
556      Finley_Mesh* fMesh;      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      } else if (order==2) {      }
562        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,          else if (order==2) {
563                       useElementsOnFace) ;        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
564                         useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
565      } else {      } else {
566        stringstream temp;        stringstream temp;
567        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
# Line 89  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};
587      int periodic[]={periodic0, periodic1};      int periodic[]={periodic0, periodic1};
588    
589      Finley_Mesh* fMesh;      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      } else if (order==2) {      }
594        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,      else if (order==2) {
595                      useElementsOnFace);        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
596      } else {                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
597        }
598        else {
599        stringstream temp;        stringstream temp;
600        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
601        setFinleyError(VALUE_ERROR,temp.str().c_str());        setFinleyError(VALUE_ERROR,temp.str().c_str());
# Line 113  namespace finley { Line 606  namespace finley {
606      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
607      return temp;      return temp;
608    }    }
609    AbstractContinuousDomain*  interval(int n0,int order,double l0,int periodic0,  
610                 int integrationOrder,    AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
611                 int useElementsOnFace)    {
612    {      Finley_Mesh* fMesh=0;
613      int numElements[]={n0};      //
614      double length[]={l0};      // extract the meshes from meshList
615      int periodic[]={periodic0};      int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
616      Finley_Mesh* fMesh;      Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
617      if (order==1) {      for (int i=0;i<numMsh;++i) {
618        fMesh=Finley_RectangularMesh_Line2(numElements, length,periodic,integrationOrder,           AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
619                       useElementsOnFace);           const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
620      } else if (order==2) {           mshes[i]=finley_meshListMember->getFinley_Mesh();
       fMesh=Finley_RectangularMesh_Line3(numElements,length,periodic,integrationOrder,  
                      useElementsOnFace);  
     } else {  
       stringstream temp;  
       temp << "Illegal interpolation order: " << order;  
       setFinleyError(VALUE_ERROR,temp.str().c_str());  
621      }      }
622      //      //
623        // merge the meshes:
624        fMesh=Finley_Mesh_merge(numMsh,mshes);
625          TMPMEMFREE(mshes);
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);
630      return temp;  
   }  
   AbstractContinuousDomain*  meshMerge(const boost::python::list& meshList)  
   {  
     AbstractContinuousDomain* temp=new MeshAdapter(0);  
631      return temp;      return temp;
632    }    }
633    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,
634              double safetyFactor,                                     double safety_factor,
635              double tolerance)                             double tolerance,
636                                           int optimize)
637    {    {
638      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
639      return temp;      //
640        // merge the meshes:
641        AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
642        //
643        // glue the faces:
644        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
645        fMesh=merged_finley_meshes->getFinley_Mesh();
646        Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
647    
648        //
649        // Convert any finley errors into a C++ exception
650        checkFinleyError();
651        return merged_meshes;
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      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
659      return temp;      //
660        // merge the meshes:
661        AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
662        //
663        // join the faces:
664        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
665        fMesh=merged_finley_meshes->getFinley_Mesh();
666        Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
667        //
668        // Convert any finley errors into a C++ exception
669        checkFinleyError();
670        return merged_meshes;
671    }    }
672    
673  }  // end of namespace    // end of namespace
674    
675    }

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

  ViewVC Help
Powered by ViewVC 1.1.26