/[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

revision 757 by woo409, Mon Jun 26 13:12:56 2006 UTC revision 1346 by ksteube, Wed Nov 14 22:48:12 2007 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"  #include "MeshAdapterFactory.h"
23  #include "FinleyError.h"  #include "FinleyError.h"
24    extern "C" {
25    #include "escript/blocktimer.h"
26    }
27    
28  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
29    
# Line 25  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; // ksteube should this be an argument to LoadMesh?
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        // create a copy of the filename to overcome the non-constness of call
61        // to Finley_Mesh_read
62        // Win32 refactor
63        char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;
64        strcpy(fName,fileName.c_str());
65    
66        printf("ksteube finley::loadMesh %s\n", fName);
67    
68        double blocktimer_start = blocktimer_time();
69        Finley_resetError();
70    
71        // Open NetCDF file for reading
72        NcAtt *attr;
73        NcVar *nc_var_temp;
74        // netCDF error handler
75        NcError err(NcError::silent_nonfatal);
76        // Create the NetCDF file.
77        NcFile dataFile(fName, NcFile::ReadOnly);
78        if (!dataFile.is_valid()) {
79          sprintf(error_msg,"loadMesh: Opening file NetCDF %s for reading failed.", fName);
80          Finley_setError(IO_ERROR,error_msg);
81          Paso_MPIInfo_free( mpi_info );
82          throw DataException("Error - loadMesh:: Could not read NetCDF file.");
83        }
84    
85        // Read NetCDF integer attributes
86        int mpi_size        = NetCDF_Get_Int_Attribute(&dataFile, fName, "mpi_size");
87        int mpi_rank        = NetCDF_Get_Int_Attribute(&dataFile, fName, "mpi_rank");
88        int numDim          = NetCDF_Get_Int_Attribute(&dataFile, fName, "numDim");
89        int order           = NetCDF_Get_Int_Attribute(&dataFile, fName, "order");
90        int reduced_order       = NetCDF_Get_Int_Attribute(&dataFile, fName, "reduced_order");
91        int numNodes        = NetCDF_Get_Int_Attribute(&dataFile, fName, "numNodes");
92        int num_Elements        = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_Elements");
93        int num_FaceElements    = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_FaceElements");
94        int num_ContactElements = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_ContactElements");
95        int num_Points      = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_Points");
96    
97        // Read mesh name
98        if (! (attr=dataFile.get_att("Name")) ) {
99          sprintf(error_msg,"Error retrieving mesh name from NetCDF file '%s'", fName);
100          throw DataException(error_msg);
101        }
102        char *name = attr->as_string(0);
103        delete attr;
104    
105        /* allocate mesh */
106        mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
107        if (Finley_noError()) {
108    
109            /* read nodes */
110            Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
111        // Nodes_Id
112            if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
113              throw DataException("Error - loadMesh:: unable to read Nodes_Id from netCDF file: " + *fName);
114            if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {
115              free(&mesh_p->Nodes->Id);
116              throw DataException("Error - loadMesh:: unable to recover Nodes_Id from NetCDF file: " + *fName);
117            }
118        // Nodes_Tag
119            if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
120              throw DataException("Error - loadMesh:: unable to read Nodes_Tag from netCDF file: " + *fName);
121            if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) ) {
122              free(&mesh_p->Nodes->Tag);
123              throw DataException("Error - loadMesh:: unable to recover Nodes_Tag from NetCDF file: " + *fName);
124            }
125        // Nodes_gDOF
126            if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
127              throw DataException("Error - loadMesh:: unable to read Nodes_gDOF from netCDF file: " + *fName);
128            if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) ) {
129              free(&mesh_p->Nodes->globalDegreesOfFreedom);
130              throw DataException("Error - loadMesh:: unable to recover Nodes_gDOF from NetCDF file: " + *fName);
131            }
132        // Nodes_gNI
133            if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
134              throw DataException("Error - loadMesh:: unable to read Nodes_gNI from netCDF file: " + *fName);
135            if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) ) {
136              free(&mesh_p->Nodes->globalNodesIndex);
137              throw DataException("Error - loadMesh:: unable to recover Nodes_gNI from NetCDF file: " + *fName);
138            }
139        // Nodes_grDfI
140            if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
141              throw DataException("Error - loadMesh:: unable to read Nodes_grDfI from netCDF file: " + *fName);
142            if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) ) {
143              free(&mesh_p->Nodes->globalReducedDOFIndex);
144              throw DataException("Error - loadMesh:: unable to recover Nodes_grDfI from NetCDF file: " + *fName);
145            }
146        // Nodes_grNI
147            if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
148              throw DataException("Error - loadMesh:: unable to read Nodes_grNI from netCDF file: " + *fName);
149            if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) ) {
150              free(&mesh_p->Nodes->globalReducedNodesIndex);
151              throw DataException("Error - loadMesh:: unable to recover Nodes_grNI from NetCDF file: " + *fName);
152            }
153        // Nodes_Coordinates
154            if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates"))) {
155              free(&mesh_p->Nodes->Coordinates);
156              throw DataException("Error - loadMesh:: unable to read Nodes_Coordinates from netCDF file: " + *fName);
157            }
158            if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) ) {
159              free(&mesh_p->Nodes->Coordinates);
160              throw DataException("Error - load:: unable to recover Nodes_Coordinates from netCDF file: " + *fName);
161            }
162    
163            /* read elements */
164            if (Finley_noError()) {
165              if (num_Elements>0) {
166             int num_Elements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_Elements_numNodes");
167             int Elements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, "Elements_TypeId");
168                 mesh_p->Elements=Finley_ElementFile_alloc((ElementTypeId)Elements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
169                 if (Finley_noError()) {
170                     Finley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
171                     mesh_p->Elements->minColor=0;
172                     mesh_p->Elements->maxColor=num_Elements-1;
173                     if (Finley_noError()) {
174                   // Elements_Id
175                       if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
176                         throw DataException("Error - loadMesh:: unable to read Elements_Id from netCDF file: " + *fName);
177                       if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) ) {
178                         free(&mesh_p->Elements->Id);
179                         throw DataException("Error - loadMesh:: unable to recover Elements_Id from NetCDF file: " + *fName);
180                       }
181                   // Elements_Tag
182                       if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
183                         throw DataException("Error - loadMesh:: unable to read Elements_Tag from netCDF file: " + *fName);
184                       if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) ) {
185                         free(&mesh_p->Elements->Tag);
186                         throw DataException("Error - loadMesh:: unable to recover Elements_Tag from NetCDF file: " + *fName);
187                       }
188                   // Elements_Owner
189                       if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
190                         throw DataException("Error - loadMesh:: unable to read Elements_Owner from netCDF file: " + *fName);
191                       if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) ) {
192                         free(&mesh_p->Elements->Owner);
193                         throw DataException("Error - loadMesh:: unable to recover Elements_Owner from NetCDF file: " + *fName);
194                       }
195                   // Elements_Color
196                       if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
197                         throw DataException("Error - loadMesh:: unable to read Elements_Color from netCDF file: " + *fName);
198                       if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) ) {
199                         free(&mesh_p->Elements->Color);
200                         throw DataException("Error - loadMesh:: unable to recover Elements_Color from NetCDF file: " + *fName);
201                       }
202                   // Elements_Nodes
203               int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
204                       if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
205                         free(&mesh_p->Elements->Nodes);
206                         throw DataException("Error - loadMesh:: unable to read Elements_Nodes from netCDF file: " + *fName);
207                       }
208                       if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
209                         free(&Elements_Nodes);
210                         throw DataException("Error - load:: unable to recover Elements_Nodes from netCDF file: " + *fName);
211                       }
212               // Copy temp array into mesh_p->Elements->Nodes
213               for (int i=0; i<num_Elements; i++) {
214                 for (int j=0; j<num_Elements_numNodes; j++) {
215                   mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
216                     = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
217                 }
218               }
219               TMPMEMFREE(Elements_Nodes);
220             }
221             }
222          } /* num_Elements>0 */
223        }
224    
225            /* get the face elements */
226            if (Finley_noError()) {
227              if (num_FaceElements>0) {
228             int num_FaceElements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_FaceElements_numNodes");
229             int FaceElements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, "FaceElements_TypeId");
230                 mesh_p->FaceElements=Finley_ElementFile_alloc((ElementTypeId)FaceElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
231                 if (Finley_noError()) {
232                     Finley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
233                     mesh_p->FaceElements->minColor=0;
234                     mesh_p->FaceElements->maxColor=num_FaceElements-1;
235                     if (Finley_noError()) {
236                   // FaceElements_Id
237                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
238                         throw DataException("Error - loadMesh:: unable to read FaceElements_Id from netCDF file: " + *fName);
239                       if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) ) {
240                         free(&mesh_p->FaceElements->Id);
241                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Id from NetCDF file: " + *fName);
242                       }
243                   // FaceElements_Tag
244                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
245                         throw DataException("Error - loadMesh:: unable to read FaceElements_Tag from netCDF file: " + *fName);
246                       if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) ) {
247                         free(&mesh_p->FaceElements->Tag);
248                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Tag from NetCDF file: " + *fName);
249                       }
250                   // FaceElements_Owner
251                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
252                         throw DataException("Error - loadMesh:: unable to read FaceElements_Owner from netCDF file: " + *fName);
253                       if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) ) {
254                         free(&mesh_p->FaceElements->Owner);
255                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Owner from NetCDF file: " + *fName);
256                       }
257                   // FaceElements_Color
258                       if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
259                         throw DataException("Error - loadMesh:: unable to read FaceElements_Color from netCDF file: " + *fName);
260                       if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) ) {
261                         free(&mesh_p->FaceElements->Color);
262                         throw DataException("Error - loadMesh:: unable to recover FaceElements_Color from NetCDF file: " + *fName);
263                       }
264                   // FaceElements_Nodes
265               int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
266                       if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
267                         free(&mesh_p->FaceElements->Nodes);
268                         throw DataException("Error - loadMesh:: unable to read FaceElements_Nodes from netCDF file: " + *fName);
269                       }
270                       if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
271                         free(&FaceElements_Nodes);
272                         throw DataException("Error - load:: unable to recover FaceElements_Nodes from netCDF file: " + *fName);
273                       }
274               // Copy temp array into mesh_p->FaceElements->Nodes
275               for (int i=0; i<num_FaceElements; i++) {
276                 for (int j=0; j<num_FaceElements_numNodes; j++) {
277                   mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]
278                     = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
279                 }
280               }
281               TMPMEMFREE(FaceElements_Nodes);
282             }
283             }
284          } /* num_FaceElements>0 */
285        }
286    
287            /* get the Contact elements */
288            if (Finley_noError()) {
289              if (num_ContactElements>0) {
290             int num_ContactElements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_ContactElements_numNodes");
291             int ContactElements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, "ContactElements_TypeId");
292                 mesh_p->ContactElements=Finley_ElementFile_alloc((ElementTypeId)ContactElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
293                 if (Finley_noError()) {
294                     Finley_ElementFile_allocTable(mesh_p->ContactElements, num_ContactElements);
295                     mesh_p->ContactElements->minColor=0;
296                     mesh_p->ContactElements->maxColor=num_ContactElements-1;
297                     if (Finley_noError()) {
298                   // ContactElements_Id
299                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Id")) )
300                         throw DataException("Error - loadMesh:: unable to read ContactElements_Id from netCDF file: " + *fName);
301                       if (! nc_var_temp->get(&mesh_p->ContactElements->Id[0], num_ContactElements) ) {
302                         free(&mesh_p->ContactElements->Id);
303                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Id from NetCDF file: " + *fName);
304                       }
305                   // ContactElements_Tag
306                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Tag")) )
307                         throw DataException("Error - loadMesh:: unable to read ContactElements_Tag from netCDF file: " + *fName);
308                       if (! nc_var_temp->get(&mesh_p->ContactElements->Tag[0], num_ContactElements) ) {
309                         free(&mesh_p->ContactElements->Tag);
310                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Tag from NetCDF file: " + *fName);
311                       }
312                   // ContactElements_Owner
313                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Owner")) )
314                         throw DataException("Error - loadMesh:: unable to read ContactElements_Owner from netCDF file: " + *fName);
315                       if (! nc_var_temp->get(&mesh_p->ContactElements->Owner[0], num_ContactElements) ) {
316                         free(&mesh_p->ContactElements->Owner);
317                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Owner from NetCDF file: " + *fName);
318                       }
319                   // ContactElements_Color
320                       if (! ( nc_var_temp = dataFile.get_var("ContactElements_Color")) )
321                         throw DataException("Error - loadMesh:: unable to read ContactElements_Color from netCDF file: " + *fName);
322                       if (! nc_var_temp->get(&mesh_p->ContactElements->Color[0], num_ContactElements) ) {
323                         free(&mesh_p->ContactElements->Color);
324                         throw DataException("Error - loadMesh:: unable to recover ContactElements_Color from NetCDF file: " + *fName);
325                       }
326                   // ContactElements_Nodes
327               int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
328                       if (!(nc_var_temp = dataFile.get_var("ContactElements_Nodes"))) {
329                         free(&mesh_p->ContactElements->Nodes);
330                         throw DataException("Error - loadMesh:: unable to read ContactElements_Nodes from netCDF file: " + *fName);
331                       }
332                       if (! nc_var_temp->get(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes) ) {
333                         free(&ContactElements_Nodes);
334                         throw DataException("Error - load:: unable to recover ContactElements_Nodes from netCDF file: " + *fName);
335                       }
336               // Copy temp array into mesh_p->ContactElements->Nodes
337               for (int i=0; i<num_ContactElements; i++) {
338                 for (int j=0; j<num_ContactElements_numNodes; j++) {
339                   mesh_p->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]
340                     = ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
341                 }
342               }
343               TMPMEMFREE(ContactElements_Nodes);
344             }
345             }
346          } /* num_ContactElements>0 */
347        }
348    
349            /* get the Points (nodal elements) */
350            if (Finley_noError()) {
351              if (num_Points>0) {
352             int Points_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, "Points_TypeId");
353                 mesh_p->Points=Finley_ElementFile_alloc((ElementTypeId)Points_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
354                 if (Finley_noError()) {
355                     Finley_ElementFile_allocTable(mesh_p->Points, num_Points);
356                     mesh_p->Points->minColor=0;
357                     mesh_p->Points->maxColor=num_Points-1;
358                     if (Finley_noError()) {
359                   // Points_Id
360                       if (! ( nc_var_temp = dataFile.get_var("Points_Id")) )
361                         throw DataException("Error - loadMesh:: unable to read Points_Id from netCDF file: " + *fName);
362                       if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points) ) {
363                         free(&mesh_p->Points->Id);
364                         throw DataException("Error - loadMesh:: unable to recover Points_Id from NetCDF file: " + *fName);
365                       }
366                   // Points_Tag
367                       if (! ( nc_var_temp = dataFile.get_var("Points_Tag")) )
368                         throw DataException("Error - loadMesh:: unable to read Points_Tag from netCDF file: " + *fName);
369                       if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points) ) {
370                         free(&mesh_p->Points->Tag);
371                         throw DataException("Error - loadMesh:: unable to recover Points_Tag from NetCDF file: " + *fName);
372                       }
373                   // Points_Owner
374                       if (! ( nc_var_temp = dataFile.get_var("Points_Owner")) )
375                         throw DataException("Error - loadMesh:: unable to read Points_Owner from netCDF file: " + *fName);
376                       if (! nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points) ) {
377                         free(&mesh_p->Points->Owner);
378                         throw DataException("Error - loadMesh:: unable to recover Points_Owner from NetCDF file: " + *fName);
379                       }
380                   // Points_Color
381                       if (! ( nc_var_temp = dataFile.get_var("Points_Color")) )
382                         throw DataException("Error - loadMesh:: unable to read Points_Color from netCDF file: " + *fName);
383                       if (! nc_var_temp->get(&mesh_p->Points->Color[0], num_Points) ) {
384                         free(&mesh_p->Points->Color);
385                         throw DataException("Error - loadMesh:: unable to recover Points_Color from NetCDF file: " + *fName);
386                       }
387                   // Points_Nodes
388               int *Points_Nodes = TMPMEMALLOC(num_Points,int);
389                       if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
390                         free(&mesh_p->Points->Nodes);
391                         throw DataException("Error - loadMesh:: unable to read Points_Nodes from netCDF file: " + *fName);
392                       }
393                       if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
394                         free(&Points_Nodes);
395                         throw DataException("Error - load:: unable to recover Points_Nodes from netCDF file: " + *fName);
396                       }
397               // Copy temp array into mesh_p->Points->Nodes
398               for (int i=0; i<num_Points; i++) {
399                 mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
400               }
401               TMPMEMFREE(Points_Nodes);
402             }
403             }
404          } /* num_Points>0 */
405        }
406    
407            /* get the name tags */
408    
409        } /* Finley_noError() after Finley_Mesh_alloc() */
410      
411        if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
412        if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
413    
414        checkFinleyError();
415        temp=new MeshAdapter(mesh_p);
416    
417        if (! Finley_noError()) {
418          Finley_Mesh_free(mesh_p);
419        }
420    
421        /* win32 refactor */
422        TMPMEMFREE(fName);
423    
424        blocktimer_increment("LoadMesh()", blocktimer_start);
425        return temp;
426    #else
427        throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
428    #endif /* USE_NETCDF */
429      }
430    
431    AbstractContinuousDomain* readMesh(const std::string& fileName,    AbstractContinuousDomain* readMesh(const std::string& fileName,
432                       int integrationOrder)                       int integrationOrder,
433                                         int reducedIntegrationOrder,
434                                         int optimize)
435    {    {
436      //      //
437      // 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
# Line 35  namespace finley { Line 440  namespace finley {
440      // Win32 refactor      // Win32 refactor
441      char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;      char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;
442      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
443        double blocktimer_start = blocktimer_time();
444    
445  #ifndef PASO_MPI      fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
446      fMesh=Finley_Mesh_read(fName,integrationOrder);      checkFinleyError();
447  #else      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
448      {      
449        stringstream temp;      /* win32 refactor */
450        temp << "Unable to read meshes from file under MPI yet...";      TMPMEMFREE(fName);
451        setFinleyError(VALUE_ERROR,temp.str().c_str());      
452      }      blocktimer_increment("ReadMesh()", blocktimer_start);
453  #endif      return temp;
454      }
455    
456      AbstractContinuousDomain* readGmsh(const std::string& fileName,
457                                         int numDim,
458                                         int integrationOrder,
459                                         int reducedIntegrationOrder,
460                                         int optimize)
461      {
462        //
463        // create a copy of the filename to overcome the non-constness of call
464        // to Finley_Mesh_read
465        Finley_Mesh* fMesh=0;
466        // Win32 refactor
467        char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;
468        strcpy(fName,fileName.c_str());
469        double blocktimer_start = blocktimer_time();
470    
471        fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
472      checkFinleyError();      checkFinleyError();
473      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
474            
475      /* win32 refactor */      /* win32 refactor */
476      TMPMEMFREE(fName);      TMPMEMFREE(fName);
477            
478        blocktimer_increment("ReadGmsh()", blocktimer_start);
479      return temp;      return temp;
480    }    }
481    
# Line 59  namespace finley { Line 484  namespace finley {
484              int periodic0,int periodic1,              int periodic0,int periodic1,
485              int periodic2,              int periodic2,
486              int integrationOrder,              int integrationOrder,
487              int useElementsOnFace)                      int reducedIntegrationOrder,
488                int useElementsOnFace,
489                int useFullElementOrder,
490                        int optimize)
491    {    {
 //     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;  
         
492      int numElements[]={n0,n1,n2};      int numElements[]={n0,n1,n2};
493      double length[]={l0,l1,l2};      double length[]={l0,l1,l2};
494      int periodic[]={periodic0, periodic1, periodic2};      int periodic[]={periodic0, periodic1, periodic2};
# Line 79  namespace finley { Line 498  namespace finley {
498      Finley_Mesh* fMesh=NULL;      Finley_Mesh* fMesh=NULL;
499    
500      if (order==1) {      if (order==1) {
501        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
502                      useElementsOnFace) ;                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
503      }      }
 #ifndef PASO_MPI  
504          else if (order==2) {          else if (order==2) {
505        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
506                       useElementsOnFace) ;                       useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
507      } else {      } else {
508        stringstream temp;        stringstream temp;
509        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
510        setFinleyError(VALUE_ERROR,temp.str().c_str());        setFinleyError(VALUE_ERROR,temp.str().c_str());
511      }      }
 #else  
         else {  
       stringstream temp;  
       temp << "type of 3D regular mesh requested is unavailable under MPI\nOR\nIllegal interpolation order: " << order;  
       setFinleyError(VALUE_ERROR,temp.str().c_str());  
     }  
 #endif  
512      //      //
513      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
514      checkFinleyError();      checkFinleyError();
# Line 108  namespace finley { Line 519  namespace finley {
519              double l0, double l1,              double l0, double l1,
520              int periodic0,int periodic1,              int periodic0,int periodic1,
521              int integrationOrder,              int integrationOrder,
522              int useElementsOnFace)                          int reducedIntegrationOrder,
523                int useElementsOnFace,
524                    int useFullElementOrder,
525                            int optimize)
526    {    {
527      int numElements[]={n0,n1};      int numElements[]={n0,n1};
528      double length[]={l0,l1};      double length[]={l0,l1};
# Line 116  namespace finley { Line 530  namespace finley {
530    
531      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
532      if (order==1) {      if (order==1) {
533        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
534                      useElementsOnFace);                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
535      }      }
 #ifndef PASO_MPI  
536      else if (order==2) {      else if (order==2) {
537        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
538                      useElementsOnFace);                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
539      }      }
 #endif  
     else {  
       stringstream temp;  
       temp << "Illegal interpolation order: " << order;  
       setFinleyError(VALUE_ERROR,temp.str().c_str());  
     }  
     //  
     // Convert any finley errors into a C++ exception  
     checkFinleyError();  
     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);  
     return temp;  
   }  
   AbstractContinuousDomain*  interval(int n0,int order,double l0,int periodic0,  
                int integrationOrder,  
                int useElementsOnFace)  
   {  
     int numElements[]={n0};  
     double length[]={l0};  
     int periodic[]={periodic0};  
     Finley_Mesh* fMesh;  
     if (order==1) {  
       fMesh=Finley_RectangularMesh_Line2(numElements, length,periodic,integrationOrder,  
                      useElementsOnFace);  
     }  
 #ifndef PASO_MPI  
     else if (order==2) {  
       fMesh=Finley_RectangularMesh_Line3(numElements,length,periodic,integrationOrder,  
                      useElementsOnFace);  
     }  
 #endif  
540      else {      else {
541        stringstream temp;        stringstream temp;
542        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
# Line 165  namespace finley { Line 548  namespace finley {
548      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
549      return temp;      return temp;
550    }    }
551    
552    AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)    AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
553    {    {
554      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
555      //      //
556      // extract the meshes from meshList      // extract the meshes from meshList
 #ifndef PASO_MPI  
557      int numMsh=boost::python::extract<int>(meshList.attr("__len__")());      int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
558      Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;      Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
559      for (int i=0;i<numMsh;++i) {      for (int i=0;i<numMsh;++i) {
# Line 181  namespace finley { Line 564  namespace finley {
564      //      //
565      // merge the meshes:      // merge the meshes:
566      fMesh=Finley_Mesh_merge(numMsh,mshes);      fMesh=Finley_Mesh_merge(numMsh,mshes);
567  #else        TMPMEMFREE(mshes);
     {  
       stringstream temp;  
       temp << "meshMerge() not available in MPI yet...";  
       setFinleyError(VALUE_ERROR,temp.str().c_str());  
     }  
 #endif  
568      //      //
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);
     TMPMEMFREE(mshes);  
572    
573      return temp;      return temp;
574    }    }
575    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,
576              double safety_factor,                                     double safety_factor,
577              double tolerance)                             double tolerance,
578                                           int optimize)
579    {    {
580      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
 #ifndef PASO_MPI  
581      //      //
582      // merge the meshes:      // merge the meshes:
583      AbstractContinuousDomain* merged_meshes=meshMerge(meshList);      AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
# Line 209  namespace finley { Line 585  namespace finley {
585      // glue the faces:      // glue the faces:
586      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
587      fMesh=merged_finley_meshes->getFinley_Mesh();      fMesh=merged_finley_meshes->getFinley_Mesh();
588      Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance);      Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
589    
590      //      //
591      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
592      checkFinleyError();      checkFinleyError();
593      return merged_meshes;      return merged_meshes;
 #else  
     {  
       stringstream temp;  
       temp << "glueFaces() not available in MPI yet...";  
       setFinleyError(VALUE_ERROR,temp.str().c_str());  
     }  
   
     //  
     // Convert any finley errors into a C++ exception  
     checkFinleyError();  
     return (AbstractContinuousDomain*)0;  
 #endif  
   
594    }    }
595    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,
596              double safety_factor,              double safety_factor,
597              double tolerance)              double tolerance,
598                            int optimize)
599    {    {
600      Finley_Mesh* fMesh=0;      Finley_Mesh* fMesh=0;
601      //      //
602      // merge the meshes:      // merge the meshes:
 #ifndef PASO_MPI  
603      AbstractContinuousDomain* merged_meshes=meshMerge(meshList);      AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
604      //      //
605      // join the faces:      // join the faces:
606      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);      const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
607      fMesh=merged_finley_meshes->getFinley_Mesh();      fMesh=merged_finley_meshes->getFinley_Mesh();
608      Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance);      Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
609      //      //
610      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
611      checkFinleyError();      checkFinleyError();
612      return merged_meshes;      return merged_meshes;
 #else  
     {  
       stringstream temp;  
       temp << "joinFaces() not available in MPI yet...";  
       setFinleyError(VALUE_ERROR,temp.str().c_str());  
     }  
     //  
     // Convert any finley errors into a C++ exception  
     checkFinleyError();  
     return (AbstractContinuousDomain*)0;  
   
 #endif  
613    }    }
614    
615  }  // end of namespace    // end of namespace
616    
617    }

Legend:
Removed from v.757  
changed lines
  Added in v.1346

  ViewVC Help
Powered by ViewVC 1.1.26