/[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 1628 by phornby, Fri Jul 11 13:12:46 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        // I don't think the strdup is needed since Paso_MPI_appendRankToFileName
62        // does it's own allocation.
63        // char *fName = Paso_MPI_appendRankToFileName(strdup(fileName.c_str()),
64        //                                             mpi_info->size,
65        //                                             mpi_info->rank);
66    
67        char *fName = Paso_MPI_appendRankToFileName(fileName.c_str(),
68                                                    mpi_info->size,
69                                                    mpi_info->rank);
70    
71        double blocktimer_start = blocktimer_time();
72        Finley_resetError();
73    
74        // Open NetCDF file for reading
75        NcAtt *attr;
76        NcVar *nc_var_temp;
77        // netCDF error handler
78        NcError err(NcError::silent_nonfatal);
79        // Create the NetCDF file.
80        NcFile dataFile(fName, NcFile::ReadOnly);
81        if (!dataFile.is_valid()) {
82          sprintf(error_msg,"loadMesh: Opening file NetCDF %s for reading failed.", fName);
83          Finley_setError(IO_ERROR,error_msg);
84          Paso_MPIInfo_free( mpi_info );
85          throw DataException(error_msg);
86        }
87    
88        // Read NetCDF integer attributes
89        int mpi_size            = NetCDF_Get_Int_Attribute(&dataFile, fName, "mpi_size");
90        int mpi_rank            = NetCDF_Get_Int_Attribute(&dataFile, fName, "mpi_rank");
91        int numDim              = NetCDF_Get_Int_Attribute(&dataFile, fName, "numDim");
92        int order               = NetCDF_Get_Int_Attribute(&dataFile, fName, "order");
93        int reduced_order           = NetCDF_Get_Int_Attribute(&dataFile, fName, "reduced_order");
94        int numNodes            = NetCDF_Get_Int_Attribute(&dataFile, fName, "numNodes");
95        int num_Elements            = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_Elements");
96        int num_FaceElements        = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_FaceElements");
97        int num_ContactElements     = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_ContactElements");
98        int num_Points          = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_Points");
99        int num_Elements_numNodes       = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_Elements_numNodes");
100        int Elements_TypeId         = NetCDF_Get_Int_Attribute(&dataFile, fName, "Elements_TypeId");
101        int num_FaceElements_numNodes   = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_FaceElements_numNodes");
102        int FaceElements_TypeId     = NetCDF_Get_Int_Attribute(&dataFile, fName, "FaceElements_TypeId");
103        int num_ContactElements_numNodes    = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_ContactElements_numNodes");
104        int ContactElements_TypeId      = NetCDF_Get_Int_Attribute(&dataFile, fName, "ContactElements_TypeId");
105        int Points_TypeId           = NetCDF_Get_Int_Attribute(&dataFile, fName, "Points_TypeId");
106        int num_Tags            = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_Tags");
107    
108        // Verify size and rank
109        if (mpi_info->size != mpi_size) {
110          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);
111          throw DataException(error_msg);
112        }
113        if (mpi_info->rank != mpi_rank) {
114          sprintf(error_msg, "Error loadMesh - The NetCDF file '%s' should be read on CPU #%d instead of %d", fName, mpi_rank, mpi_info->rank);
115          throw DataException(error_msg);
116        }
117    
118        // Read mesh name
119        if (! (attr=dataFile.get_att("Name")) ) {
120          sprintf(error_msg,"Error retrieving mesh name from NetCDF file '%s'", fName);
121          throw DataException(error_msg);
122        }
123        char *name = attr->as_string(0);
124        delete attr;
125    
126        /* allocate mesh */
127        mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
128        if (Finley_noError()) {
129    
130            /* read nodes */
131            Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
132        // Nodes_Id
133            if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
134              throw DataException("Error - loadMesh:: unable to read Nodes_Id from netCDF file: " + *fName);
135            if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {
136              free(&mesh_p->Nodes->Id);
137              throw DataException("Error - loadMesh:: unable to recover Nodes_Id from NetCDF file: " + *fName);
138            }
139        // Nodes_Tag
140            if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
141              throw DataException("Error - loadMesh:: unable to read Nodes_Tag from netCDF file: " + *fName);
142            if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) ) {
143              free(&mesh_p->Nodes->Tag);
144              throw DataException("Error - loadMesh:: unable to recover Nodes_Tag from NetCDF file: " + *fName);
145            }
146        // Nodes_gDOF
147            if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
148              throw DataException("Error - loadMesh:: unable to read Nodes_gDOF from netCDF file: " + *fName);
149            if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) ) {
150              free(&mesh_p->Nodes->globalDegreesOfFreedom);
151              throw DataException("Error - loadMesh:: unable to recover Nodes_gDOF from NetCDF file: " + *fName);
152            }
153        // Nodes_gNI
154            if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
155              throw DataException("Error - loadMesh:: unable to read Nodes_gNI from netCDF file: " + *fName);
156            if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) ) {
157              free(&mesh_p->Nodes->globalNodesIndex);
158              throw DataException("Error - loadMesh:: unable to recover Nodes_gNI from NetCDF file: " + *fName);
159            }
160        // Nodes_grDfI
161            if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
162              throw DataException("Error - loadMesh:: unable to read Nodes_grDfI from netCDF file: " + *fName);
163            if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) ) {
164              free(&mesh_p->Nodes->globalReducedDOFIndex);
165              throw DataException("Error - loadMesh:: unable to recover Nodes_grDfI from NetCDF file: " + *fName);
166            }
167        // Nodes_grNI
168            if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
169              throw DataException("Error - loadMesh:: unable to read Nodes_grNI from netCDF file: " + *fName);
170            if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) ) {
171              free(&mesh_p->Nodes->globalReducedNodesIndex);
172              throw DataException("Error - loadMesh:: unable to recover Nodes_grNI from NetCDF file: " + *fName);
173            }
174        // Nodes_Coordinates
175            if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates"))) {
176              free(&mesh_p->Nodes->Coordinates);
177              throw DataException("Error - loadMesh:: unable to read Nodes_Coordinates from netCDF file: " + *fName);
178            }
179            if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) ) {
180              free(&mesh_p->Nodes->Coordinates);
181              throw DataException("Error - load:: unable to recover Nodes_Coordinates from netCDF file: " + *fName);
182            }
183        // Nodes_DofDistribution
184        int *first_component = TMPMEMALLOC(mpi_size+1,index_t);
185            if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) )
186              throw DataException("Error - loadMesh:: unable to read Nodes_DofDistribution from netCDF file: " + *fName);
187            if (! nc_var_temp->get(&first_component[0], mpi_size+1) ) {
188              free(&first_component);
189              throw DataException("Error - loadMesh:: unable to recover Nodes_DofDistribution from NetCDF file: " + *fName);
190            }
191        mesh_p->Nodes->degreesOfFreedomDistribution=Paso_Distribution_alloc(mesh_p->Nodes->MPIInfo,first_component,1,0);
192        TMPMEMFREE(first_component);
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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                         free(&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                  free(&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, mesh_p->Nodes->degreesOfFreedomDistribution->first_component);
455    
456        checkFinleyError();
457        temp=new MeshAdapter(mesh_p);
458    
459        if (! Finley_noError()) {
460          Finley_Mesh_free(mesh_p);
461        }
462    
463        /* win32 refactor */
464        TMPMEMFREE(fName);
465    
466        blocktimer_increment("LoadMesh()", blocktimer_start);
467        return temp;
468    #else
469        throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
470    #endif /* USE_NETCDF */
471      }
472    
473    AbstractContinuousDomain* readMesh(const std::string& fileName,    AbstractContinuousDomain* readMesh(const std::string& fileName,
474                       int integrationOrder)                       int integrationOrder,
475                                         int reducedIntegrationOrder,
476                                         int optimize)
477    {    {
478      //      //
479      // 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
480      // to Finley_Mesh_read      // to Finley_Mesh_read
481      char fName[fileName.size()+1];      Finley_Mesh* fMesh=0;
482        // Win32 refactor
483        if( fileName.size() == 0 )
484        {
485           throw DataException("Null file name!");
486        }
487    
488        char *fName = TMPMEMALLOC(fileName.size()+1,char);
489        
490      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
491      Finley_Mesh* fMesh=Finley_Mesh_read(fName,integrationOrder);      double blocktimer_start = blocktimer_time();
492    
493        fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
494      checkFinleyError();      checkFinleyError();
495      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
496        
497        /* win32 refactor */
498        TMPMEMFREE(fName);
499        
500        blocktimer_increment("ReadMesh()", blocktimer_start);
501        return temp;
502      }
503    
504      AbstractContinuousDomain* readMeshMPI(const std::string& fileName,
505                         int integrationOrder,
506                                         int reducedIntegrationOrder,
507                                         int optimize)
508      {
509        //
510        // create a copy of the filename to overcome the non-constness of call
511        // to Finley_Mesh_read
512        Finley_Mesh* fMesh=0;
513        // Win32 refactor
514        if( fileName.size() == 0 )
515        {
516           throw DataException("Null file name!");
517        }
518    
519        char *fName = TMPMEMALLOC(fileName.size()+1,char);
520        
521        strcpy(fName,fileName.c_str());
522        double blocktimer_start = blocktimer_time();
523    
524        fMesh=Finley_Mesh_read_MPI(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
525        checkFinleyError();
526        AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
527        
528        /* win32 refactor */
529        TMPMEMFREE(fName);
530        
531        blocktimer_increment("ReadMesh()", blocktimer_start);
532        return temp;
533      }
534    
535      AbstractContinuousDomain* readGmsh(const std::string& fileName,
536                                         int numDim,
537                                         int integrationOrder,
538                                         int reducedIntegrationOrder,
539                                         int optimize)
540      {
541        //
542        // create a copy of the filename to overcome the non-constness of call
543        // to Finley_Mesh_read
544        Finley_Mesh* fMesh=0;
545        // Win32 refactor
546        if( fileName.size() == 0 )
547        {
548           throw DataException("Null file name!");
549        }
550    
551        char *fName = TMPMEMALLOC(fileName.size()+1,char);
552        
553        strcpy(fName,fileName.c_str());
554        double blocktimer_start = blocktimer_time();
555    
556        fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
557        checkFinleyError();
558        AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
559        
560        /* win32 refactor */
561        TMPMEMFREE(fName);
562        
563        blocktimer_increment("ReadGmsh()", blocktimer_start);
564      return temp;      return temp;
565    }    }
566    
# Line 50  namespace finley { Line 569  namespace finley {
569              int periodic0,int periodic1,              int periodic0,int periodic1,
570              int periodic2,              int periodic2,
571              int integrationOrder,              int integrationOrder,
572              int useElementsOnFace)                      int reducedIntegrationOrder,
573                int useElementsOnFace,
574                int useFullElementOrder,
575                        int optimize)
576    {    {
 //     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;  
         
577      int numElements[]={n0,n1,n2};      int numElements[]={n0,n1,n2};
578      double length[]={l0,l1,l2};      double length[]={l0,l1,l2};
579      int periodic[]={periodic0, periodic1, periodic2};      int periodic[]={periodic0, periodic1, periodic2};
580    
581      //      //
582      // linearInterpolation      // linearInterpolation
583      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=NULL;
584    
585      if (order==1) {      if (order==1) {
586        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
587                      useElementsOnFace) ;                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
588      } else if (order==2) {      }
589        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,          else if (order==2) {
590                       useElementsOnFace) ;        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
591                         useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
592      } else {      } else {
593        stringstream temp;        stringstream temp;
594        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
# Line 89  namespace finley { Line 604  namespace finley {
604              double l0, double l1,              double l0, double l1,
605              int periodic0,int periodic1,              int periodic0,int periodic1,
606              int integrationOrder,              int integrationOrder,
607              int useElementsOnFace)                          int reducedIntegrationOrder,
608                int useElementsOnFace,
609                    int useFullElementOrder,
610                            int optimize)
611    {    {
612      int numElements[]={n0,n1};      int numElements[]={n0,n1};
613      double length[]={l0,l1};      double length[]={l0,l1};
614      int periodic[]={periodic0, periodic1};      int periodic[]={periodic0, periodic1};
615    
616      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=0;
617      if (order==1) {      if (order==1) {
618        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
619                      useElementsOnFace);                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
620      } else if (order==2) {      }
621        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,      else if (order==2) {
622                      useElementsOnFace);        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
623      } else {                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
624        }
625        else {
626        stringstream temp;        stringstream temp;
627        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
628        setFinleyError(VALUE_ERROR,temp.str().c_str());        setFinleyError(VALUE_ERROR,temp.str().c_str());
# Line 113  namespace finley { Line 633  namespace finley {
633      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
634      return temp;      return temp;
635    }    }
636    AbstractContinuousDomain*  interval(int n0,int order,double l0,int periodic0,  
637                 int integrationOrder,    AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
638                 int useElementsOnFace)    {
639    {      Finley_Mesh* fMesh=0;
640      int numElements[]={n0};      //
641      double length[]={l0};      // extract the meshes from meshList
642      int periodic[]={periodic0};      int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
643      Finley_Mesh* fMesh;      Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
644      if (order==1) {      for (int i=0;i<numMsh;++i) {
645        fMesh=Finley_RectangularMesh_Line2(numElements, length,periodic,integrationOrder,           AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
646                       useElementsOnFace);           const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
647      } 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());  
648      }      }
649      //      //
650        // merge the meshes:
651        fMesh=Finley_Mesh_merge(numMsh,mshes);
652          TMPMEMFREE(mshes);
653        //
654      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
655      checkFinleyError();      checkFinleyError();
656      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
657      return temp;  
   }  
   AbstractContinuousDomain*  meshMerge(const boost::python::list& meshList)  
   {  
     AbstractContinuousDomain* temp=new MeshAdapter(0);  
658      return temp;      return temp;
659    }    }
660    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,
661              double safetyFactor,                                     double safety_factor,
662              double tolerance)                             double tolerance,
663                                           int optimize)
664    {    {
665      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
666      return temp;      //
667        // merge the meshes:
668        AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
669        //
670        // glue the faces:
671        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
672        fMesh=merged_finley_meshes->getFinley_Mesh();
673        Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
674    
675        //
676        // Convert any finley errors into a C++ exception
677        checkFinleyError();
678        return merged_meshes;
679    }    }
680    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,
681              double safety_factor,              double safety_factor,
682              double tolerance)              double tolerance,
683                            int optimize)
684    {    {
685      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
686      return temp;      //
687        // merge the meshes:
688        AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
689        //
690        // join the faces:
691        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
692        fMesh=merged_finley_meshes->getFinley_Mesh();
693        Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
694        //
695        // Convert any finley errors into a C++ exception
696        checkFinleyError();
697        return merged_meshes;
698    }    }
699    
700  }  // end of namespace    // end of namespace
701    
702    }

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

  ViewVC Help
Powered by ViewVC 1.1.26