/[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 1990 by ksteube, Fri Nov 7 04:19:07 2008 UTC
# Line 1  Line 1 
 /* $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.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
1    
2    /*******************************************************
3    *
4    * Copyright (c) 2003-2008 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
14    
15    #ifdef PASO_MPI
16    #include <mpi.h>
17    #endif
18    #ifdef USE_NETCDF
19    #include <netcdfcpp.h>
20    #endif
21    #include "MeshAdapterFactory.h"
22    #include "FinleyError.h"
23  extern "C" {  extern "C" {
24  #include "finley/finleyC/Finley.h"  #include "escript/blocktimer.h"
 #include "finley/finleyC/Mesh.h"  
 #include "finley/finleyC/RectangularMesh.h"  
25  }  }
 #include "finley/CPPAdapter/FinleyError.h"  
 #include "finley/CPPAdapter/MeshAdapterFactory.h"  
26    
27    #include <boost/python/extract.hpp>
28    
   
 #include <iostream>  
29  #include <sstream>  #include <sstream>
30    
31  using namespace std;  using namespace std;
# Line 31  using namespace escript; Line 33  using namespace escript;
33    
34  namespace finley {  namespace finley {
35    
36    AbstractContinuousDomain* readMesh(const std::string& fileName,  #ifdef USE_NETCDF
37                       int integrationOrder)    // A convenience method to retrieve an integer attribute from a NetCDF file
38      int NetCDF_Get_Int_Attribute(NcFile *dataFile, char *fName, char *attr_name) {
39        NcAtt *attr;
40        char error_msg[LenErrorMsg_MAX];
41        if (! (attr=dataFile->get_att(attr_name)) ) {
42          sprintf(error_msg,"Error retrieving integer attribute '%s' from NetCDF file '%s'", attr_name, fName);
43          throw DataException(error_msg);
44        }
45        int temp = attr->as_int(0);
46        delete attr;
47        return(temp);
48      }
49    #endif
50    
51    //   AbstractContinuousDomain* loadMesh(const std::string& fileName)
52      Domain_ptr loadMesh(const std::string& fileName)
53      {
54    #ifdef USE_NETCDF
55        Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
56        AbstractContinuousDomain* temp;
57        Finley_Mesh *mesh_p=NULL;
58        char error_msg[LenErrorMsg_MAX];
59    
60        char *fName = Paso_MPI_appendRankToFileName(fileName.c_str(),
61                                                    mpi_info->size,
62                                                    mpi_info->rank);
63    
64        double blocktimer_start = blocktimer_time();
65        Finley_resetError();
66        int *first_DofComponent, *first_NodeComponent;
67    
68        // Open NetCDF file for reading
69        NcAtt *attr;
70        NcVar *nc_var_temp;
71        // netCDF error handler
72        NcError err(NcError::silent_nonfatal);
73        // Create the NetCDF file.
74        NcFile dataFile(fName, NcFile::ReadOnly);
75        if (!dataFile.is_valid()) {
76          sprintf(error_msg,"loadMesh: Opening file NetCDF %s for reading failed.", fName);
77          Finley_setError(IO_ERROR,error_msg);
78          Paso_MPIInfo_free( mpi_info );
79          throw DataException(error_msg);
80        }
81    
82        // Read NetCDF integer attributes
83        int mpi_size            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_size");
84        int mpi_rank            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_rank");
85        int numDim              = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numDim");
86        int order               = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"order");
87        int reduced_order           = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"reduced_order");
88        int numNodes            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numNodes");
89        int num_Elements            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements");
90        int num_FaceElements        = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements");
91        int num_ContactElements     = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_ContactElements");
92        int num_Points          = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Points");
93        int num_Elements_numNodes       = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements_numNodes");
94        int Elements_TypeId         = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Elements_TypeId");
95        int num_FaceElements_numNodes   = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements_numNodes");
96        int FaceElements_TypeId     = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"FaceElements_TypeId");
97        int num_ContactElements_numNodes    = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_ContactElements_numNodes");
98        int ContactElements_TypeId      = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"ContactElements_TypeId");
99        int Points_TypeId           = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Points_TypeId");
100        int num_Tags            = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Tags");
101    
102        // Verify size and rank
103        if (mpi_info->size != mpi_size) {
104          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);
105          throw DataException(error_msg);
106        }
107        if (mpi_info->rank != mpi_rank) {
108          sprintf(error_msg, "Error loadMesh - The NetCDF file '%s' should be read on CPU #%d instead of %d", fName, mpi_rank, mpi_info->rank);
109          throw DataException(error_msg);
110        }
111    
112        // Read mesh name
113        if (! (attr=dataFile.get_att("Name")) ) {
114          sprintf(error_msg,"Error retrieving mesh name from NetCDF file '%s'", fName);
115          throw DataException(error_msg);
116        }
117        char *name = attr->as_string(0);
118        delete attr;
119    
120        /* allocate mesh */
121        mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
122        if (Finley_noError()) {
123    
124            /* read nodes */
125            Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
126        // Nodes_Id
127            if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
128              throw DataException("Error - loadMesh:: unable to read Nodes_Id from netCDF file: " + *fName);
129            if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {
130              TMPMEMFREE(mesh_p->Nodes->Id);
131              throw DataException("Error - loadMesh:: unable to recover Nodes_Id from NetCDF file: " + *fName);
132            }
133        // Nodes_Tag
134            if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
135              throw DataException("Error - loadMesh:: unable to read Nodes_Tag from netCDF file: " + *fName);
136            if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) ) {
137              TMPMEMFREE(mesh_p->Nodes->Tag);
138              throw DataException("Error - loadMesh:: unable to recover Nodes_Tag from NetCDF file: " + *fName);
139            }
140        // Nodes_gDOF
141            if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
142              throw DataException("Error - loadMesh:: unable to read Nodes_gDOF from netCDF file: " + *fName);
143            if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) ) {
144              TMPMEMFREE(mesh_p->Nodes->globalDegreesOfFreedom);
145              throw DataException("Error - loadMesh:: unable to recover Nodes_gDOF from NetCDF file: " + *fName);
146            }
147        // Nodes_gNI
148            if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
149              throw DataException("Error - loadMesh:: unable to read Nodes_gNI from netCDF file: " + *fName);
150            if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) ) {
151              TMPMEMFREE(mesh_p->Nodes->globalNodesIndex);
152              throw DataException("Error - loadMesh:: unable to recover Nodes_gNI from NetCDF file: " + *fName);
153            }
154        // Nodes_grDfI
155            if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
156              throw DataException("Error - loadMesh:: unable to read Nodes_grDfI from netCDF file: " + *fName);
157            if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) ) {
158              TMPMEMFREE(mesh_p->Nodes->globalReducedDOFIndex);
159              throw DataException("Error - loadMesh:: unable to recover Nodes_grDfI from NetCDF file: " + *fName);
160            }
161        // Nodes_grNI
162            if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
163              throw DataException("Error - loadMesh:: unable to read Nodes_grNI from netCDF file: " + *fName);
164            if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) ) {
165              TMPMEMFREE(mesh_p->Nodes->globalReducedNodesIndex);
166              throw DataException("Error - loadMesh:: unable to recover Nodes_grNI from NetCDF file: " + *fName);
167            }
168        // Nodes_Coordinates
169            if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates"))) {
170              TMPMEMFREE(mesh_p->Nodes->Coordinates);
171              throw DataException("Error - loadMesh:: unable to read Nodes_Coordinates from netCDF file: " + *fName);
172            }
173            if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) ) {
174              TMPMEMFREE(mesh_p->Nodes->Coordinates);
175              throw DataException("Error - load:: unable to recover Nodes_Coordinates from netCDF file: " + *fName);
176            }
177        // Nodes_DofDistribution
178        first_DofComponent = TMPMEMALLOC(mpi_size+1,index_t);
179            if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) )
180              throw DataException("Error - loadMesh:: unable to read Nodes_DofDistribution from netCDF file: " + *fName);
181            if (! nc_var_temp->get(&first_DofComponent[0], mpi_size+1) ) {
182              throw DataException("Error - loadMesh:: unable to recover Nodes_DofDistribution from NetCDF file: " + *fName);
183            }
184    
185        // Nodes_NodeDistribution
186        first_NodeComponent = TMPMEMALLOC(mpi_size+1,index_t);
187            if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) )
188              throw DataException("Error - loadMesh:: unable to read Nodes_NodeDistribution from netCDF file: " + *fName);
189            if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
190              throw DataException("Error - loadMesh:: unable to recover Nodes_NodeDistribution from NetCDF file: " + *fName);
191            }
192    
193            /* read elements */
194            if (Finley_noError()) {
195              mesh_p->Elements=Finley_ElementFile_alloc((ElementTypeId)Elements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
196              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
197              mesh_p->Elements->minColor=0;
198              mesh_p->Elements->maxColor=num_Elements-1;
199              if (num_Elements>0) {
200                     if (Finley_noError()) {
201                   // Elements_Id
202                       if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
203                         throw DataException("Error - loadMesh:: unable to read Elements_Id from netCDF file: " + *fName);
204                       if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) ) {
205                         TMPMEMFREE(mesh_p->Elements->Id);
206                         throw DataException("Error - loadMesh:: unable to recover Elements_Id from NetCDF file: " + *fName);
207                       }
208                   // Elements_Tag
209                       if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
210                         throw DataException("Error - loadMesh:: unable to read Elements_Tag from netCDF file: " + *fName);
211                       if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) ) {
212                         TMPMEMFREE(mesh_p->Elements->Tag);
213                         throw DataException("Error - loadMesh:: unable to recover Elements_Tag from NetCDF file: " + *fName);
214                       }
215                   // Elements_Owner
216                       if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
217                         throw DataException("Error - loadMesh:: unable to read Elements_Owner from netCDF file: " + *fName);
218                       if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) ) {
219                         TMPMEMFREE(mesh_p->Elements->Owner);
220                         throw DataException("Error - loadMesh:: unable to recover Elements_Owner from NetCDF file: " + *fName);
221                       }
222                   // Elements_Color
223                       if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
224                         throw DataException("Error - loadMesh:: unable to read Elements_Color from netCDF file: " + *fName);
225                       if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) ) {
226                         TMPMEMFREE(mesh_p->Elements->Color);
227                         throw DataException("Error - loadMesh:: unable to recover Elements_Color from NetCDF file: " + *fName);
228                       }
229                   // Elements_Nodes
230               int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
231                       if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
232                         TMPMEMFREE(mesh_p->Elements->Nodes);
233                         throw DataException("Error - loadMesh:: unable to read Elements_Nodes from netCDF file: " + *fName);
234                       }
235                       if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
236                         TMPMEMFREE(Elements_Nodes);
237                         throw DataException("Error - load:: unable to recover Elements_Nodes from netCDF file: " + *fName);
238                       }
239               // Copy temp array into mesh_p->Elements->Nodes
240               for (int i=0; i<num_Elements; i++) {
241                 for (int j=0; j<num_Elements_numNodes; j++) {
242                   mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
243                     = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
244                 }
245               }
246               TMPMEMFREE(Elements_Nodes);
247             }
248          } /* num_Elements>0 */
249        }
250    
251            /* get the face elements */
252            if (Finley_noError()) {
253              mesh_p->FaceElements=Finley_ElementFile_alloc((ElementTypeId)FaceElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
254              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
255              mesh_p->FaceElements->minColor=0;
256              mesh_p->FaceElements->maxColor=num_FaceElements-1;
257              if (num_FaceElements>0) {
258                     if (Finley_noError()) {
259                   // FaceElements_Id
260                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
261                         throw DataException("Error - loadMesh:: unable to read FaceElements_Id from netCDF file: " + *fName);
262                       if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) ) {
263                         TMPMEMFREE(mesh_p->FaceElements->Id);
264                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Id from NetCDF file: " + *fName);
265                       }
266                   // FaceElements_Tag
267                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
268                         throw DataException("Error - loadMesh:: unable to read FaceElements_Tag from netCDF file: " + *fName);
269                       if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) ) {
270                         TMPMEMFREE(mesh_p->FaceElements->Tag);
271                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Tag from NetCDF file: " + *fName);
272                       }
273                   // FaceElements_Owner
274                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
275                         throw DataException("Error - loadMesh:: unable to read FaceElements_Owner from netCDF file: " + *fName);
276                       if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) ) {
277                         TMPMEMFREE(mesh_p->FaceElements->Owner);
278                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Owner from NetCDF file: " + *fName);
279                       }
280                   // FaceElements_Color
281                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
282                         throw DataException("Error - loadMesh:: unable to read FaceElements_Color from netCDF file: " + *fName);
283                       if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) ) {
284                         TMPMEMFREE(mesh_p->FaceElements->Color);
285                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Color from NetCDF file: " + *fName);
286                       }
287                   // FaceElements_Nodes
288               int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
289                       if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
290                         TMPMEMFREE(mesh_p->FaceElements->Nodes);
291                         throw DataException("Error - loadMesh:: unable to read FaceElements_Nodes from netCDF file: " + *fName);
292                       }
293                       if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
294                         TMPMEMFREE(FaceElements_Nodes);
295                         throw DataException("Error - load:: unable to recover FaceElements_Nodes from netCDF file: " + *fName);
296                       }
297               // Copy temp array into mesh_p->FaceElements->Nodes
298               for (int i=0; i<num_FaceElements; i++) {
299                 for (int j=0; j<num_FaceElements_numNodes; j++) {
300                   mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]
301                     = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
302                 }
303               }
304               TMPMEMFREE(FaceElements_Nodes);
305             }
306          } /* num_FaceElements>0 */
307        }
308    
309            /* get the Contact elements */
310            if (Finley_noError()) {
311              mesh_p->ContactElements=Finley_ElementFile_alloc((ElementTypeId)ContactElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
312              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->ContactElements, num_ContactElements);
313              mesh_p->ContactElements->minColor=0;
314              mesh_p->ContactElements->maxColor=num_ContactElements-1;
315              if (num_ContactElements>0) {
316                     if (Finley_noError()) {
317                   // ContactElements_Id
318                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Id")) )
319                         throw DataException("Error - loadMesh:: unable to read ContactElements_Id from netCDF file: " + *fName);
320                       if (! nc_var_temp->get(&mesh_p->ContactElements->Id[0], num_ContactElements) ) {
321                         TMPMEMFREE(mesh_p->ContactElements->Id);
322                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Id from NetCDF file: " + *fName);
323                       }
324                   // ContactElements_Tag
325                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Tag")) )
326                         throw DataException("Error - loadMesh:: unable to read ContactElements_Tag from netCDF file: " + *fName);
327                       if (! nc_var_temp->get(&mesh_p->ContactElements->Tag[0], num_ContactElements) ) {
328                         TMPMEMFREE(mesh_p->ContactElements->Tag);
329                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Tag from NetCDF file: " + *fName);
330                       }
331                   // ContactElements_Owner
332                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Owner")) )
333                         throw DataException("Error - loadMesh:: unable to read ContactElements_Owner from netCDF file: " + *fName);
334                       if (! nc_var_temp->get(&mesh_p->ContactElements->Owner[0], num_ContactElements) ) {
335                         TMPMEMFREE(mesh_p->ContactElements->Owner);
336                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Owner from NetCDF file: " + *fName);
337                       }
338                   // ContactElements_Color
339                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Color")) )
340                         throw DataException("Error - loadMesh:: unable to read ContactElements_Color from netCDF file: " + *fName);
341                       if (! nc_var_temp->get(&mesh_p->ContactElements->Color[0], num_ContactElements) ) {
342                         TMPMEMFREE(mesh_p->ContactElements->Color);
343                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Color from NetCDF file: " + *fName);
344                       }
345                   // ContactElements_Nodes
346               int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
347                       if (!(nc_var_temp = dataFile.get_var("ContactElements_Nodes"))) {
348                         TMPMEMFREE(mesh_p->ContactElements->Nodes);
349                         throw DataException("Error - loadMesh:: unable to read ContactElements_Nodes from netCDF file: " + *fName);
350                       }
351                       if (! nc_var_temp->get(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes) ) {
352                         TMPMEMFREE(ContactElements_Nodes);
353                         throw DataException("Error - load:: unable to recover ContactElements_Nodes from netCDF file: " + *fName);
354                       }
355               // Copy temp array into mesh_p->ContactElements->Nodes
356               for (int i=0; i<num_ContactElements; i++) {
357                 for (int j=0; j<num_ContactElements_numNodes; j++) {
358                   mesh_p->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]
359                     = ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
360                 }
361               }
362               TMPMEMFREE(ContactElements_Nodes);
363             }
364          } /* num_ContactElements>0 */
365        }
366    
367            /* get the Points (nodal elements) */
368            if (Finley_noError()) {
369              mesh_p->Points=Finley_ElementFile_alloc((ElementTypeId)Points_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
370              if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Points, num_Points);
371              mesh_p->Points->minColor=0;
372              mesh_p->Points->maxColor=num_Points-1;
373              if (num_Points>0) {
374                 if (Finley_noError()) {
375                   // Points_Id
376                       if (! ( nc_var_temp = dataFile.get_var("Points_Id")) )
377                         throw DataException("Error - loadMesh:: unable to read Points_Id from netCDF file: " + *fName);
378                       if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points) ) {
379                         TMPMEMFREE(mesh_p->Points->Id);
380                         throw DataException("Error - loadMesh:: unable to recover Points_Id from NetCDF file: " + *fName);
381                       }
382                   // Points_Tag
383                       if (! ( nc_var_temp = dataFile.get_var("Points_Tag")) )
384                         throw DataException("Error - loadMesh:: unable to read Points_Tag from netCDF file: " + *fName);
385                       if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points) ) {
386                         TMPMEMFREE(mesh_p->Points->Tag);
387                         throw DataException("Error - loadMesh:: unable to recover Points_Tag from NetCDF file: " + *fName);
388                       }
389                   // Points_Owner
390                       if (! ( nc_var_temp = dataFile.get_var("Points_Owner")) )
391                         throw DataException("Error - loadMesh:: unable to read Points_Owner from netCDF file: " + *fName);
392                       if (! nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points) ) {
393                         TMPMEMFREE(mesh_p->Points->Owner);
394                         throw DataException("Error - loadMesh:: unable to recover Points_Owner from NetCDF file: " + *fName);
395                       }
396                   // Points_Color
397                       if (! ( nc_var_temp = dataFile.get_var("Points_Color")) )
398                         throw DataException("Error - loadMesh:: unable to read Points_Color from netCDF file: " + *fName);
399                       if (! nc_var_temp->get(&mesh_p->Points->Color[0], num_Points) ) {
400                         TMPMEMFREE(mesh_p->Points->Color);
401                         throw DataException("Error - loadMesh:: unable to recover Points_Color from NetCDF file: " + *fName);
402                       }
403                   // Points_Nodes
404               int *Points_Nodes = TMPMEMALLOC(num_Points,int);
405                       if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
406                         TMPMEMFREE(mesh_p->Points->Nodes);
407                         throw DataException("Error - loadMesh:: unable to read Points_Nodes from netCDF file: " + *fName);
408                       }
409                       if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
410                         TMPMEMFREE(Points_Nodes);
411                         throw DataException("Error - load:: unable to recover Points_Nodes from netCDF file: " + *fName);
412                       }
413               // Copy temp array into mesh_p->Points->Nodes
414               for (int i=0; i<num_Points; i++) {
415                 mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
416               }
417               TMPMEMFREE(Points_Nodes);
418             }
419          } /* num_Points>0 */
420        }
421    
422            /* get the tags */
423            if (Finley_noError()) {
424              if (num_Tags>0) {
425                // Temp storage to gather node IDs
426                int *Tags_keys = TMPMEMALLOC(num_Tags, int);
427                char name_temp[4096];
428            int i;
429    
430            // Tags_keys
431                if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) )
432                  throw DataException("Error - loadMesh:: unable to read Tags_keys from netCDF file: " + *fName);
433                if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
434                  TMPMEMFREE(Tags_keys);
435                  throw DataException("Error - loadMesh:: unable to recover Tags_keys from NetCDF file: " + *fName);
436                }
437            for (i=0; i<num_Tags; i++) {
438                  // Retrieve tag name
439                  sprintf(name_temp, "Tags_name_%d", i);
440                  if (! (attr=dataFile.get_att(name_temp)) ) {
441                    sprintf(error_msg,"Error retrieving tag name from NetCDF file '%s'", fName);
442                    throw DataException(error_msg);
443                  }
444                  char *name = attr->as_string(0);
445                  delete attr;
446                  Finley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
447            }
448          }
449        }
450      
451            if (Finley_noError()) Finley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
452            TMPMEMFREE(first_DofComponent);
453            TMPMEMFREE(first_NodeComponent);
454    
455        } /* Finley_noError() after Finley_Mesh_alloc() */
456    
457        checkFinleyError();
458        temp=new MeshAdapter(mesh_p);
459    
460        if (! Finley_noError()) {
461          Finley_Mesh_free(mesh_p);
462        }
463    
464        /* win32 refactor */
465        TMPMEMFREE(fName);
466    
467        blocktimer_increment("LoadMesh()", blocktimer_start);
468        return temp->getPtr();
469    #else
470        throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
471    #endif /* USE_NETCDF */
472      }
473    
474      Domain_ptr readMesh(const std::string& fileName,
475                         int integrationOrder,
476                                         int reducedIntegrationOrder,
477                                         int optimize)
478    {    {
479      //      //
480      // 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
481      // to Finley_Mesh_read      // to Finley_Mesh_read
482      char fName[fileName.size()+1];      Finley_Mesh* fMesh=0;
483        // Win32 refactor
484        if( fileName.size() == 0 )
485        {
486           throw DataException("Null file name!");
487        }
488    
489        char *fName = TMPMEMALLOC(fileName.size()+1,char);
490        
491        strcpy(fName,fileName.c_str());
492        double blocktimer_start = blocktimer_time();
493    
494        fMesh=Finley_Mesh_read_MPI(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
495        checkFinleyError();
496        AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
497        
498        /* win32 refactor */
499        TMPMEMFREE(fName);
500        
501        blocktimer_increment("ReadMesh()", blocktimer_start);
502        return temp->getPtr();
503      }
504    
505      Domain_ptr readGmsh(const std::string& fileName,
506                                         int numDim,
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());      strcpy(fName,fileName.c_str());
524      Finley_Mesh* fMesh=Finley_Mesh_read(fName,integrationOrder);      double blocktimer_start = blocktimer_time();
525    
526        fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
527      checkFinleyError();      checkFinleyError();
528      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
529      return temp;      
530        /* win32 refactor */
531        TMPMEMFREE(fName);
532        
533        blocktimer_increment("ReadGmsh()", blocktimer_start);
534        return temp->getPtr();
535    }    }
536    
537    AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,  /*  AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,*/
538      Domain_ptr brick(int n0,int n1,int n2,int order,
539              double l0,double l1,double l2,              double l0,double l1,double l2,
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 83  namespace finley { Line 569  namespace finley {
569      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
570      checkFinleyError();      checkFinleyError();
571      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
572      return temp;      return temp->getPtr();
573    }    }
574    AbstractContinuousDomain*  rectangle(int n0,int n1,int order,  
575    /*  AbstractContinuousDomain*  rectangle(int n0,int n1,int order,*/
576      Domain_ptr  rectangle(int n0,int n1,int order,
577              double l0, double l1,              double l0, double l1,
578              int periodic0,int periodic1,              int periodic0,int periodic1,
579              int integrationOrder,              int integrationOrder,
580              int useElementsOnFace)                          int reducedIntegrationOrder,
581                int useElementsOnFace,
582                    int useFullElementOrder,
583                            int optimize)
584    {    {
585      int numElements[]={n0,n1};      int numElements[]={n0,n1};
586      double length[]={l0,l1};      double length[]={l0,l1};
587      int periodic[]={periodic0, periodic1};      int periodic[]={periodic0, periodic1};
588    
589      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=0;
590      if (order==1) {      if (order==1) {
591        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
592                      useElementsOnFace);                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
593      } else if (order==2) {      }
594        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,      else if (order==2) {
595                      useElementsOnFace);        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
596      } else {                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
597        }
598        else {
599        stringstream temp;        stringstream temp;
600        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
601        setFinleyError(VALUE_ERROR,temp.str().c_str());        setFinleyError(VALUE_ERROR,temp.str().c_str());
# Line 111  namespace finley { Line 604  namespace finley {
604      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
605      checkFinleyError();      checkFinleyError();
606      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
607      return temp;      return temp->getPtr();
608    }    }
609    AbstractContinuousDomain*  interval(int n0,int order,double l0,int periodic0,  
610                 int integrationOrder,    Domain_ptr meshMerge(const boost::python::list& meshList)
611                 int useElementsOnFace)    {
612    {      Finley_Mesh* fMesh=0;
613      int numElements[]={n0};      //
614      double length[]={l0};      // extract the meshes from meshList
615      int periodic[]={periodic0};      int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
616      Finley_Mesh* fMesh;      Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
617      if (order==1) {      for (int i=0;i<numMsh;++i) {
618        fMesh=Finley_RectangularMesh_Line2(numElements, length,periodic,integrationOrder,           AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
619                       useElementsOnFace);           const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
620      } else if (order==2) {           mshes[i]=finley_meshListMember->getFinley_Mesh();
       fMesh=Finley_RectangularMesh_Line3(numElements,length,periodic,integrationOrder,  
                      useElementsOnFace);  
     } else {  
       stringstream temp;  
       temp << "Illegal interpolation order: " << order;  
       setFinleyError(VALUE_ERROR,temp.str().c_str());  
621      }      }
622      //      //
623        // merge the meshes:
624        fMesh=Finley_Mesh_merge(numMsh,mshes);
625          TMPMEMFREE(mshes);
626        //
627      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
628      checkFinleyError();      checkFinleyError();
629      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
630      return temp;  
631    }      return temp->getPtr();
   AbstractContinuousDomain*  meshMerge(const boost::python::list& meshList)  
   {  
     AbstractContinuousDomain* temp=new MeshAdapter(0);  
     return temp;  
632    }    }
633    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,  
634              double safetyFactor,    Domain_ptr  glueFaces(const boost::python::list& meshList,
635              double tolerance)                                     double safety_factor,
636                               double tolerance,
637                                           int optimize)
638    {    {
639      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
640      return temp;      //
641        // merge the meshes:
642        Domain_ptr merged_meshes=meshMerge(meshList);
643    
644        //
645        // glue the faces:
646        const MeshAdapter* merged_finley_meshes=dynamic_cast<const MeshAdapter*>(merged_meshes.get());
647        fMesh=merged_finley_meshes->getFinley_Mesh();
648        Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
649    
650        //
651        // Convert any finley errors into a C++ exception
652        checkFinleyError();
653        return merged_meshes->getPtr();
654    }    }
655    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,    Domain_ptr  joinFaces(const boost::python::list& meshList,
656              double safety_factor,              double safety_factor,
657              double tolerance)              double tolerance,
658                            int optimize)
659    {    {
660      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
661      return temp;      //
662        // merge the meshes:
663        Domain_ptr merged_meshes=meshMerge(meshList);
664        //
665        // join the faces:
666        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes.get());
667        fMesh=merged_finley_meshes->getFinley_Mesh();
668        Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
669        //
670        // Convert any finley errors into a C++ exception
671        checkFinleyError();
672        return merged_meshes->getPtr();
673    }    }
674    
675  }  // end of namespace    // end of namespace
676    
677    }

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

  ViewVC Help
Powered by ViewVC 1.1.26