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

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

  ViewVC Help
Powered by ViewVC 1.1.26