/[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 1755 by ksteube, Mon Sep 8 01:34:40 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    { int i; printf("ksteube read  Nodes_DofDistribution:"); for(i=0; i<mpi_size+1; i++) printf(" %d", first_DofComponent[i]); printf("\n"); }
186    
187        // Nodes_NodeDistribution
188        first_NodeComponent = TMPMEMALLOC(mpi_size+1,index_t);
189            if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) )
190              throw DataException("Error - loadMesh:: unable to read Nodes_NodeDistribution from netCDF file: " + *fName);
191            if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
192              throw DataException("Error - loadMesh:: unable to recover Nodes_NodeDistribution from NetCDF file: " + *fName);
193            }
194    { int i; printf("ksteube read  Nodes_NodeDistribution:"); for(i=0; i<mpi_size+1; i++) printf(" %d", first_NodeComponent[i]); printf("\n"); }
195    
196            /* read elements */
197            if (Finley_noError()) {
198              mesh_p->Elements=Finley_ElementFile_alloc((ElementTypeId)Elements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
199              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
200              mesh_p->Elements->minColor=0;
201              mesh_p->Elements->maxColor=num_Elements-1;
202              if (num_Elements>0) {
203                     if (Finley_noError()) {
204                   // Elements_Id
205                       if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
206                         throw DataException("Error - loadMesh:: unable to read Elements_Id from netCDF file: " + *fName);
207                       if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) ) {
208                         TMPMEMFREE(mesh_p->Elements->Id);
209                         throw DataException("Error - loadMesh:: unable to recover Elements_Id from NetCDF file: " + *fName);
210                       }
211                   // Elements_Tag
212                       if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
213                         throw DataException("Error - loadMesh:: unable to read Elements_Tag from netCDF file: " + *fName);
214                       if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) ) {
215                         TMPMEMFREE(mesh_p->Elements->Tag);
216                         throw DataException("Error - loadMesh:: unable to recover Elements_Tag from NetCDF file: " + *fName);
217                       }
218                   // Elements_Owner
219                       if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
220                         throw DataException("Error - loadMesh:: unable to read Elements_Owner from netCDF file: " + *fName);
221                       if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) ) {
222                         TMPMEMFREE(mesh_p->Elements->Owner);
223                         throw DataException("Error - loadMesh:: unable to recover Elements_Owner from NetCDF file: " + *fName);
224                       }
225                   // Elements_Color
226                       if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
227                         throw DataException("Error - loadMesh:: unable to read Elements_Color from netCDF file: " + *fName);
228                       if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) ) {
229                         TMPMEMFREE(mesh_p->Elements->Color);
230                         throw DataException("Error - loadMesh:: unable to recover Elements_Color from NetCDF file: " + *fName);
231                       }
232                   // Elements_Nodes
233               int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
234                       if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
235                         TMPMEMFREE(mesh_p->Elements->Nodes);
236                         throw DataException("Error - loadMesh:: unable to read Elements_Nodes from netCDF file: " + *fName);
237                       }
238                       if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
239                         TMPMEMFREE(Elements_Nodes);
240                         throw DataException("Error - load:: unable to recover Elements_Nodes from netCDF file: " + *fName);
241                       }
242               // Copy temp array into mesh_p->Elements->Nodes
243               for (int i=0; i<num_Elements; i++) {
244                 for (int j=0; j<num_Elements_numNodes; j++) {
245                   mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
246                     = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
247                 }
248               }
249               TMPMEMFREE(Elements_Nodes);
250             }
251          } /* num_Elements>0 */
252        }
253    
254            /* get the face elements */
255            if (Finley_noError()) {
256              mesh_p->FaceElements=Finley_ElementFile_alloc((ElementTypeId)FaceElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
257              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
258              mesh_p->FaceElements->minColor=0;
259              mesh_p->FaceElements->maxColor=num_FaceElements-1;
260              if (num_FaceElements>0) {
261                     if (Finley_noError()) {
262                   // FaceElements_Id
263                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
264                         throw DataException("Error - loadMesh:: unable to read FaceElements_Id from netCDF file: " + *fName);
265                       if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) ) {
266                         TMPMEMFREE(mesh_p->FaceElements->Id);
267                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Id from NetCDF file: " + *fName);
268                       }
269                   // FaceElements_Tag
270                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
271                         throw DataException("Error - loadMesh:: unable to read FaceElements_Tag from netCDF file: " + *fName);
272                       if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) ) {
273                         TMPMEMFREE(mesh_p->FaceElements->Tag);
274                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Tag from NetCDF file: " + *fName);
275                       }
276                   // FaceElements_Owner
277                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
278                         throw DataException("Error - loadMesh:: unable to read FaceElements_Owner from netCDF file: " + *fName);
279                       if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) ) {
280                         TMPMEMFREE(mesh_p->FaceElements->Owner);
281                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Owner from NetCDF file: " + *fName);
282                       }
283                   // FaceElements_Color
284                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
285                         throw DataException("Error - loadMesh:: unable to read FaceElements_Color from netCDF file: " + *fName);
286                       if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) ) {
287                         TMPMEMFREE(mesh_p->FaceElements->Color);
288                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Color from NetCDF file: " + *fName);
289                       }
290                   // FaceElements_Nodes
291               int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
292                       if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
293                         TMPMEMFREE(mesh_p->FaceElements->Nodes);
294                         throw DataException("Error - loadMesh:: unable to read FaceElements_Nodes from netCDF file: " + *fName);
295                       }
296                       if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
297                         TMPMEMFREE(FaceElements_Nodes);
298                         throw DataException("Error - load:: unable to recover FaceElements_Nodes from netCDF file: " + *fName);
299                       }
300               // Copy temp array into mesh_p->FaceElements->Nodes
301               for (int i=0; i<num_FaceElements; i++) {
302                 for (int j=0; j<num_FaceElements_numNodes; j++) {
303                   mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]
304                     = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
305                 }
306               }
307               TMPMEMFREE(FaceElements_Nodes);
308             }
309          } /* num_FaceElements>0 */
310        }
311    
312            /* get the Contact elements */
313            if (Finley_noError()) {
314              mesh_p->ContactElements=Finley_ElementFile_alloc((ElementTypeId)ContactElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
315              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->ContactElements, num_ContactElements);
316              mesh_p->ContactElements->minColor=0;
317              mesh_p->ContactElements->maxColor=num_ContactElements-1;
318              if (num_ContactElements>0) {
319                     if (Finley_noError()) {
320                   // ContactElements_Id
321                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Id")) )
322                         throw DataException("Error - loadMesh:: unable to read ContactElements_Id from netCDF file: " + *fName);
323                       if (! nc_var_temp->get(&mesh_p->ContactElements->Id[0], num_ContactElements) ) {
324                         TMPMEMFREE(mesh_p->ContactElements->Id);
325                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Id from NetCDF file: " + *fName);
326                       }
327                   // ContactElements_Tag
328                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Tag")) )
329                         throw DataException("Error - loadMesh:: unable to read ContactElements_Tag from netCDF file: " + *fName);
330                       if (! nc_var_temp->get(&mesh_p->ContactElements->Tag[0], num_ContactElements) ) {
331                         TMPMEMFREE(mesh_p->ContactElements->Tag);
332                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Tag from NetCDF file: " + *fName);
333                       }
334                   // ContactElements_Owner
335                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Owner")) )
336                         throw DataException("Error - loadMesh:: unable to read ContactElements_Owner from netCDF file: " + *fName);
337                       if (! nc_var_temp->get(&mesh_p->ContactElements->Owner[0], num_ContactElements) ) {
338                         TMPMEMFREE(mesh_p->ContactElements->Owner);
339                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Owner from NetCDF file: " + *fName);
340                       }
341                   // ContactElements_Color
342                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Color")) )
343                         throw DataException("Error - loadMesh:: unable to read ContactElements_Color from netCDF file: " + *fName);
344                       if (! nc_var_temp->get(&mesh_p->ContactElements->Color[0], num_ContactElements) ) {
345                         TMPMEMFREE(mesh_p->ContactElements->Color);
346                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Color from NetCDF file: " + *fName);
347                       }
348                   // ContactElements_Nodes
349               int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
350                       if (!(nc_var_temp = dataFile.get_var("ContactElements_Nodes"))) {
351                         TMPMEMFREE(mesh_p->ContactElements->Nodes);
352                         throw DataException("Error - loadMesh:: unable to read ContactElements_Nodes from netCDF file: " + *fName);
353                       }
354                       if (! nc_var_temp->get(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes) ) {
355                         TMPMEMFREE(ContactElements_Nodes);
356                         throw DataException("Error - load:: unable to recover ContactElements_Nodes from netCDF file: " + *fName);
357                       }
358               // Copy temp array into mesh_p->ContactElements->Nodes
359               for (int i=0; i<num_ContactElements; i++) {
360                 for (int j=0; j<num_ContactElements_numNodes; j++) {
361                   mesh_p->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]
362                     = ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
363                 }
364               }
365               TMPMEMFREE(ContactElements_Nodes);
366             }
367          } /* num_ContactElements>0 */
368        }
369    
370            /* get the Points (nodal elements) */
371            if (Finley_noError()) {
372              mesh_p->Points=Finley_ElementFile_alloc((ElementTypeId)Points_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
373              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Points, num_Points);
374              mesh_p->Points->minColor=0;
375              mesh_p->Points->maxColor=num_Points-1;
376              if (num_Points>0) {
377                 if (Finley_noError()) {
378                   // Points_Id
379                       if (! ( nc_var_temp = dataFile.get_var("Points_Id")) )
380                         throw DataException("Error - loadMesh:: unable to read Points_Id from netCDF file: " + *fName);
381                       if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points) ) {
382                         TMPMEMFREE(mesh_p->Points->Id);
383                         throw DataException("Error - loadMesh:: unable to recover Points_Id from NetCDF file: " + *fName);
384                       }
385                   // Points_Tag
386                       if (! ( nc_var_temp = dataFile.get_var("Points_Tag")) )
387                         throw DataException("Error - loadMesh:: unable to read Points_Tag from netCDF file: " + *fName);
388                       if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points) ) {
389                         TMPMEMFREE(mesh_p->Points->Tag);
390                         throw DataException("Error - loadMesh:: unable to recover Points_Tag from NetCDF file: " + *fName);
391                       }
392                   // Points_Owner
393                       if (! ( nc_var_temp = dataFile.get_var("Points_Owner")) )
394                         throw DataException("Error - loadMesh:: unable to read Points_Owner from netCDF file: " + *fName);
395                       if (! nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points) ) {
396                         TMPMEMFREE(mesh_p->Points->Owner);
397                         throw DataException("Error - loadMesh:: unable to recover Points_Owner from NetCDF file: " + *fName);
398                       }
399                   // Points_Color
400                       if (! ( nc_var_temp = dataFile.get_var("Points_Color")) )
401                         throw DataException("Error - loadMesh:: unable to read Points_Color from netCDF file: " + *fName);
402                       if (! nc_var_temp->get(&mesh_p->Points->Color[0], num_Points) ) {
403                         TMPMEMFREE(mesh_p->Points->Color);
404                         throw DataException("Error - loadMesh:: unable to recover Points_Color from NetCDF file: " + *fName);
405                       }
406                   // Points_Nodes
407               int *Points_Nodes = TMPMEMALLOC(num_Points,int);
408                       if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
409                         TMPMEMFREE(mesh_p->Points->Nodes);
410                         throw DataException("Error - loadMesh:: unable to read Points_Nodes from netCDF file: " + *fName);
411                       }
412                       if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
413                         TMPMEMFREE(Points_Nodes);
414                         throw DataException("Error - load:: unable to recover Points_Nodes from netCDF file: " + *fName);
415                       }
416               // Copy temp array into mesh_p->Points->Nodes
417               for (int i=0; i<num_Points; i++) {
418                 mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
419               }
420               TMPMEMFREE(Points_Nodes);
421             }
422          } /* num_Points>0 */
423        }
424    
425            /* get the tags */
426            if (Finley_noError()) {
427              if (num_Tags>0) {
428                // Temp storage to gather node IDs
429                int *Tags_keys = TMPMEMALLOC(num_Tags, int);
430                char name_temp[4096];
431            int i;
432    
433            // Tags_keys
434                if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) )
435                  throw DataException("Error - loadMesh:: unable to read Tags_keys from netCDF file: " + *fName);
436                if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
437                  TMPMEMFREE(Tags_keys);
438                  throw DataException("Error - loadMesh:: unable to recover Tags_keys from NetCDF file: " + *fName);
439                }
440            for (i=0; i<num_Tags; i++) {
441                  // Retrieve tag name
442                  sprintf(name_temp, "Tags_name_%d", i);
443                  if (! (attr=dataFile.get_att(name_temp)) ) {
444                    sprintf(error_msg,"Error retrieving tag name from NetCDF file '%s'", fName);
445                    throw DataException(error_msg);
446                  }
447                  char *name = attr->as_string(0);
448                  delete attr;
449                  Finley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
450            }
451          }
452        }
453    
454        } /* Finley_noError() after Finley_Mesh_alloc() */
455      
456        if (Finley_noError()) Finley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
457        TMPMEMFREE(first_DofComponent);
458        TMPMEMFREE(first_NodeComponent);
459    
460        checkFinleyError();
461        temp=new MeshAdapter(mesh_p);
462    
463        if (! Finley_noError()) {
464          Finley_Mesh_free(mesh_p);
465        }
466    
467        /* win32 refactor */
468        TMPMEMFREE(fName);
469    
470        blocktimer_increment("LoadMesh()", blocktimer_start);
471        return temp;
472    #else
473        throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
474    #endif /* USE_NETCDF */
475      }
476    
477    AbstractContinuousDomain* readMesh(const std::string& fileName,    AbstractContinuousDomain* readMesh(const std::string& fileName,
478                       int integrationOrder)                       int integrationOrder,
479                                         int reducedIntegrationOrder,
480                                         int optimize)
481    {    {
482      //      //
483      // 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
484      // to Finley_Mesh_read      // to Finley_Mesh_read
485      char fName[fileName.size()+1];      Finley_Mesh* fMesh=0;
486        // Win32 refactor
487        if( fileName.size() == 0 )
488        {
489           throw DataException("Null file name!");
490        }
491    
492        char *fName = TMPMEMALLOC(fileName.size()+1,char);
493        
494      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
495      Finley_Mesh* fMesh=Finley_Mesh_read(fName,integrationOrder);      double blocktimer_start = blocktimer_time();
496    
497        fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
498      checkFinleyError();      checkFinleyError();
499      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
500        
501        /* win32 refactor */
502        TMPMEMFREE(fName);
503        
504        blocktimer_increment("ReadMesh()", blocktimer_start);
505        return temp;
506      }
507    
508      AbstractContinuousDomain* readMeshMPI(const std::string& fileName,
509                         int integrationOrder,
510                                         int reducedIntegrationOrder,
511                                         int optimize)
512      {
513        //
514        // create a copy of the filename to overcome the non-constness of call
515        // to Finley_Mesh_read
516        Finley_Mesh* fMesh=0;
517        // Win32 refactor
518        if( fileName.size() == 0 )
519        {
520           throw DataException("Null file name!");
521        }
522    
523        char *fName = TMPMEMALLOC(fileName.size()+1,char);
524        
525        strcpy(fName,fileName.c_str());
526        double blocktimer_start = blocktimer_time();
527    
528        fMesh=Finley_Mesh_read_MPI(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
529        checkFinleyError();
530        AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
531        
532        /* win32 refactor */
533        TMPMEMFREE(fName);
534        
535        blocktimer_increment("ReadMesh()", blocktimer_start);
536        return temp;
537      }
538    
539      AbstractContinuousDomain* readGmsh(const std::string& fileName,
540                                         int numDim,
541                                         int integrationOrder,
542                                         int reducedIntegrationOrder,
543                                         int optimize)
544      {
545        //
546        // create a copy of the filename to overcome the non-constness of call
547        // to Finley_Mesh_read
548        Finley_Mesh* fMesh=0;
549        // Win32 refactor
550        if( fileName.size() == 0 )
551        {
552           throw DataException("Null file name!");
553        }
554    
555        char *fName = TMPMEMALLOC(fileName.size()+1,char);
556        
557        strcpy(fName,fileName.c_str());
558        double blocktimer_start = blocktimer_time();
559    
560        fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
561        checkFinleyError();
562        AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
563        
564        /* win32 refactor */
565        TMPMEMFREE(fName);
566        
567        blocktimer_increment("ReadGmsh()", blocktimer_start);
568      return temp;      return temp;
569    }    }
570    
# Line 50  namespace finley { Line 573  namespace finley {
573              int periodic0,int periodic1,              int periodic0,int periodic1,
574              int periodic2,              int periodic2,
575              int integrationOrder,              int integrationOrder,
576              int useElementsOnFace)                      int reducedIntegrationOrder,
577                int useElementsOnFace,
578                int useFullElementOrder,
579                        int optimize)
580    {    {
 //     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;  
         
581      int numElements[]={n0,n1,n2};      int numElements[]={n0,n1,n2};
582      double length[]={l0,l1,l2};      double length[]={l0,l1,l2};
583      int periodic[]={periodic0, periodic1, periodic2};      int periodic[]={periodic0, periodic1, periodic2};
584    
585      //      //
586      // linearInterpolation      // linearInterpolation
587      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=NULL;
588    
589      if (order==1) {      if (order==1) {
590        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
591                      useElementsOnFace) ;                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
592      } else if (order==2) {      }
593        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,          else if (order==2) {
594                       useElementsOnFace) ;        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
595                         useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
596      } else {      } else {
597        stringstream temp;        stringstream temp;
598        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
# Line 89  namespace finley { Line 608  namespace finley {
608              double l0, double l1,              double l0, double l1,
609              int periodic0,int periodic1,              int periodic0,int periodic1,
610              int integrationOrder,              int integrationOrder,
611              int useElementsOnFace)                          int reducedIntegrationOrder,
612                int useElementsOnFace,
613                    int useFullElementOrder,
614                            int optimize)
615    {    {
616      int numElements[]={n0,n1};      int numElements[]={n0,n1};
617      double length[]={l0,l1};      double length[]={l0,l1};
618      int periodic[]={periodic0, periodic1};      int periodic[]={periodic0, periodic1};
619    
620      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=0;
621      if (order==1) {      if (order==1) {
622        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
623                      useElementsOnFace);                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
624      } else if (order==2) {      }
625        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,      else if (order==2) {
626                      useElementsOnFace);        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
627      } else {                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
628        }
629        else {
630        stringstream temp;        stringstream temp;
631        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
632        setFinleyError(VALUE_ERROR,temp.str().c_str());        setFinleyError(VALUE_ERROR,temp.str().c_str());
# Line 113  namespace finley { Line 637  namespace finley {
637      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
638      return temp;      return temp;
639    }    }
640    AbstractContinuousDomain*  interval(int n0,int order,double l0,int periodic0,  
641                 int integrationOrder,    AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
642                 int useElementsOnFace)    {
643    {      Finley_Mesh* fMesh=0;
644      int numElements[]={n0};      //
645      double length[]={l0};      // extract the meshes from meshList
646      int periodic[]={periodic0};      int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
647      Finley_Mesh* fMesh;      Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
648      if (order==1) {      for (int i=0;i<numMsh;++i) {
649        fMesh=Finley_RectangularMesh_Line2(numElements, length,periodic,integrationOrder,           AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
650                       useElementsOnFace);           const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
651      } 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());  
652      }      }
653      //      //
654        // merge the meshes:
655        fMesh=Finley_Mesh_merge(numMsh,mshes);
656          TMPMEMFREE(mshes);
657        //
658      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
659      checkFinleyError();      checkFinleyError();
660      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
661      return temp;  
   }  
   AbstractContinuousDomain*  meshMerge(const boost::python::list& meshList)  
   {  
     AbstractContinuousDomain* temp=new MeshAdapter(0);  
662      return temp;      return temp;
663    }    }
664    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,
665              double safetyFactor,                                     double safety_factor,
666              double tolerance)                             double tolerance,
667                                           int optimize)
668    {    {
669      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
670      return temp;      //
671        // merge the meshes:
672        AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
673        //
674        // glue the faces:
675        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
676        fMesh=merged_finley_meshes->getFinley_Mesh();
677        Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
678    
679        //
680        // Convert any finley errors into a C++ exception
681        checkFinleyError();
682        return merged_meshes;
683    }    }
684    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,
685              double safety_factor,              double safety_factor,
686              double tolerance)              double tolerance,
687                            int optimize)
688    {    {
689      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
690      return temp;      //
691        // merge the meshes:
692        AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
693        //
694        // join the faces:
695        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
696        fMesh=merged_finley_meshes->getFinley_Mesh();
697        Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
698        //
699        // Convert any finley errors into a C++ exception
700        checkFinleyError();
701        return merged_meshes;
702    }    }
703    
704  }  // end of namespace    // end of namespace
705    
706    }

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

  ViewVC Help
Powered by ViewVC 1.1.26