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

Diff of /branches/more_shared_ptrs_from_1812/finley/src/CPPAdapter/MeshAdapterFactory.cpp

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

trunk/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 1807 by ksteube, Thu Sep 25 01:04:51 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
482        // to Finley_Mesh_read
483        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());
493        double blocktimer_start = blocktimer_time();
494    
495        fMesh=Finley_Mesh_read_MPI(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
496        checkFinleyError();
497        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* readGmsh(const std::string& fileName,
507                                         int numDim,
508                                         int integrationOrder,
509                                         int reducedIntegrationOrder,
510                                         int optimize)
511    {    {
512      //      //
513      // create a copy of the filename to overcome the non-constness of call      // create a copy of the filename to overcome the non-constness of call
514      // to Finley_Mesh_read      // to Finley_Mesh_read
515      char fName[fileName.size()+1];      Finley_Mesh* fMesh=0;
516        // Win32 refactor
517        if( fileName.size() == 0 )
518        {
519           throw DataException("Null file name!");
520        }
521    
522        char *fName = TMPMEMALLOC(fileName.size()+1,char);
523        
524      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
525      Finley_Mesh* fMesh=Finley_Mesh_read(fName,integrationOrder);      double blocktimer_start = blocktimer_time();
526    
527        fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
528      checkFinleyError();      checkFinleyError();
529      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
530        
531        /* win32 refactor */
532        TMPMEMFREE(fName);
533        
534        blocktimer_increment("ReadGmsh()", blocktimer_start);
535      return temp;      return temp;
536    }    }
537    
# Line 50  namespace finley { Line 540  namespace finley {
540              int periodic0,int periodic1,              int periodic0,int periodic1,
541              int periodic2,              int periodic2,
542              int integrationOrder,              int integrationOrder,
543              int useElementsOnFace)                      int reducedIntegrationOrder,
544                int useElementsOnFace,
545                int useFullElementOrder,
546                        int optimize)
547    {    {
 //     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;  
         
548      int numElements[]={n0,n1,n2};      int numElements[]={n0,n1,n2};
549      double length[]={l0,l1,l2};      double length[]={l0,l1,l2};
550      int periodic[]={periodic0, periodic1, periodic2};      int periodic[]={periodic0, periodic1, periodic2};
551    
552      //      //
553      // linearInterpolation      // linearInterpolation
554      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=NULL;
555    
556      if (order==1) {      if (order==1) {
557        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
558                      useElementsOnFace) ;                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
559      } else if (order==2) {      }
560        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,          else if (order==2) {
561                       useElementsOnFace) ;        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
562                         useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
563      } else {      } else {
564        stringstream temp;        stringstream temp;
565        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
# Line 89  namespace finley { Line 575  namespace finley {
575              double l0, double l1,              double l0, double l1,
576              int periodic0,int periodic1,              int periodic0,int periodic1,
577              int integrationOrder,              int integrationOrder,
578              int useElementsOnFace)                          int reducedIntegrationOrder,
579                int useElementsOnFace,
580                    int useFullElementOrder,
581                            int optimize)
582    {    {
583      int numElements[]={n0,n1};      int numElements[]={n0,n1};
584      double length[]={l0,l1};      double length[]={l0,l1};
585      int periodic[]={periodic0, periodic1};      int periodic[]={periodic0, periodic1};
586    
587      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=0;
588      if (order==1) {      if (order==1) {
589        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
590                      useElementsOnFace);                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
591      } else if (order==2) {      }
592        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,      else if (order==2) {
593                      useElementsOnFace);        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
594      } else {                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
595        }
596        else {
597        stringstream temp;        stringstream temp;
598        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
599        setFinleyError(VALUE_ERROR,temp.str().c_str());        setFinleyError(VALUE_ERROR,temp.str().c_str());
# Line 113  namespace finley { Line 604  namespace finley {
604      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
605      return temp;      return temp;
606    }    }
607    AbstractContinuousDomain*  interval(int n0,int order,double l0,int periodic0,  
608                 int integrationOrder,    AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
609                 int useElementsOnFace)    {
610    {      Finley_Mesh* fMesh=0;
611      int numElements[]={n0};      //
612      double length[]={l0};      // extract the meshes from meshList
613      int periodic[]={periodic0};      int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
614      Finley_Mesh* fMesh;      Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
615      if (order==1) {      for (int i=0;i<numMsh;++i) {
616        fMesh=Finley_RectangularMesh_Line2(numElements, length,periodic,integrationOrder,           AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
617                       useElementsOnFace);           const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
618      } 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());  
619      }      }
620      //      //
621        // merge the meshes:
622        fMesh=Finley_Mesh_merge(numMsh,mshes);
623          TMPMEMFREE(mshes);
624        //
625      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
626      checkFinleyError();      checkFinleyError();
627      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
628      return temp;  
   }  
   AbstractContinuousDomain*  meshMerge(const boost::python::list& meshList)  
   {  
     AbstractContinuousDomain* temp=new MeshAdapter(0);  
629      return temp;      return temp;
630    }    }
631    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,
632              double safetyFactor,                                     double safety_factor,
633              double tolerance)                             double tolerance,
634                                           int optimize)
635    {    {
636      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
637      return temp;      //
638        // merge the meshes:
639        AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
640        //
641        // glue the faces:
642        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
643        fMesh=merged_finley_meshes->getFinley_Mesh();
644        Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
645    
646        //
647        // Convert any finley errors into a C++ exception
648        checkFinleyError();
649        return merged_meshes;
650    }    }
651    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,
652              double safety_factor,              double safety_factor,
653              double tolerance)              double tolerance,
654                            int optimize)
655    {    {
656      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
657      return temp;      //
658        // merge the meshes:
659        AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
660        //
661        // join the faces:
662        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
663        fMesh=merged_finley_meshes->getFinley_Mesh();
664        Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
665        //
666        // Convert any finley errors into a C++ exception
667        checkFinleyError();
668        return merged_meshes;
669    }    }
670    
671  }  // end of namespace    // end of namespace
672    
673    }

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

  ViewVC Help
Powered by ViewVC 1.1.26