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

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

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

trunk/esys2/finley/src/CPPAdapter/MeshAdapterFactory.cpp revision 102 by jgs, Wed Dec 15 07:08:39 2004 UTC trunk/finley/src/CPPAdapter/MeshAdapterFactory.cpp revision 1345 by ksteube, Wed Nov 14 07:53:34 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"
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        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        // create a copy of the filename to overcome the non-constness of call
60        // to Finley_Mesh_read
61        // Win32 refactor
62        char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;
63        strcpy(fName,fileName.c_str());
64    
65        printf("ksteube finley::loadMesh %s\n", fName);
66    
67        double blocktimer_start = blocktimer_time();
68        Finley_resetError();
69    
70        // Open NetCDF file for reading
71        NcAtt *attr;
72        NcVar *nc_var_temp;
73        // netCDF error handler
74        NcError err(NcError::silent_nonfatal);
75        // Create the NetCDF file.
76        NcFile dataFile(fName, NcFile::ReadOnly);
77        if (!dataFile.is_valid()) {
78          sprintf(error_msg,"loadMesh: Opening file NetCDF %s for reading failed.", fName);
79          Finley_setError(IO_ERROR,error_msg);
80          Paso_MPIInfo_free( mpi_info );
81          throw DataException("Error - loadMesh:: Could not read NetCDF file.");
82        }
83    
84        // Read NetCDF integer attributes
85        int mpi_size        = NetCDF_Get_Int_Attribute(&dataFile, fName, "mpi_size");
86        int mpi_rank        = NetCDF_Get_Int_Attribute(&dataFile, fName, "mpi_rank");
87        int numDim          = NetCDF_Get_Int_Attribute(&dataFile, fName, "numDim");
88        int order           = NetCDF_Get_Int_Attribute(&dataFile, fName, "order");
89        int reduced_order       = NetCDF_Get_Int_Attribute(&dataFile, fName, "reduced_order");
90        int numNodes        = NetCDF_Get_Int_Attribute(&dataFile, fName, "numNodes");
91        int num_Elements        = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_Elements");
92        int num_FaceElements    = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_FaceElements");
93        int num_ContactElements = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_ContactElements");
94        int num_Points      = NetCDF_Get_Int_Attribute(&dataFile, fName, "num_Points");
95    
96        // Read mesh name
97        if (! (attr=dataFile.get_att("Name")) ) {
98          sprintf(error_msg,"Error retrieving mesh name from NetCDF file '%s'", fName);
99          throw DataException(error_msg);
100        }
101        char *name = attr->as_string(0);
102        delete attr;
103    
104        /* allocate mesh */
105        mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
106        if (Finley_noError()) {
107    
108            /* read nodes */
109            Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
110        // Nodes_Id
111            if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
112              throw DataException("Error - loadMesh:: unable to read Nodes_Id from netCDF file: " + *fName);
113            if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {
114              free(&mesh_p->Nodes->Id);
115              throw DataException("Error - loadMesh:: unable to recover Nodes_Id from NetCDF file: " + *fName);
116            }
117    // printf("ksteube Nodes_Id: "); for (int i=0; i<numNodes; i++) { printf(" %d", mesh_p->Nodes->Id[i]); } printf("\n");
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    #if 0 /* Not yet...finish the above first */
164            /* read elements */
165            if (Finley_noError()) {
166                 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
167                 if (Finley_noError()) {
168                     Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
169                     mesh_p->Elements->minColor=0;
170                     mesh_p->Elements->maxColor=numEle-1;
171                     if (Finley_noError()) {
172             }
173             }
174        }
175    #endif
176    
177            /* get the face elements */
178    
179            /* get the Contact face element */
180    
181            /* get the nodal element */
182    
183            /* get the name tags */
184    
185        } /* Finley_noError() after Finley_Mesh_alloc() */
186      
187    #if 0 /* Not yet...finish the above first */
188        if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
189        if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
190    #endif
191    
192        checkFinleyError();
193        temp=new MeshAdapter(mesh_p);
194    
195        if (! Finley_noError()) {
196          Finley_Mesh_free(mesh_p);
197        }
198    
199        /* win32 refactor */
200        TMPMEMFREE(fName);
201    
202        blocktimer_increment("LoadMesh()", blocktimer_start);
203        return temp;
204    #else
205        throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
206    #endif /* USE_NETCDF */
207      }
208    
209    AbstractContinuousDomain* readMesh(const std::string& fileName,    AbstractContinuousDomain* readMesh(const std::string& fileName,
210                       int integrationOrder)                       int integrationOrder,
211                                         int reducedIntegrationOrder,
212                                         int optimize)
213    {    {
214      //      //
215      // 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
216      // to Finley_Mesh_read      // to Finley_Mesh_read
217      char fName[fileName.size()+1];      Finley_Mesh* fMesh=0;
218        // Win32 refactor
219        char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;
220      strcpy(fName,fileName.c_str());      strcpy(fName,fileName.c_str());
221      Finley_Mesh* fMesh=Finley_Mesh_read(fName,integrationOrder);      double blocktimer_start = blocktimer_time();
222    
223        fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
224      checkFinleyError();      checkFinleyError();
225      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
226        
227        /* win32 refactor */
228        TMPMEMFREE(fName);
229        
230        blocktimer_increment("ReadMesh()", blocktimer_start);
231        return temp;
232      }
233    
234      AbstractContinuousDomain* readGmsh(const std::string& fileName,
235                                         int numDim,
236                                         int integrationOrder,
237                                         int reducedIntegrationOrder,
238                                         int optimize)
239      {
240        //
241        // create a copy of the filename to overcome the non-constness of call
242        // to Finley_Mesh_read
243        Finley_Mesh* fMesh=0;
244        // Win32 refactor
245        char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;
246        strcpy(fName,fileName.c_str());
247        double blocktimer_start = blocktimer_time();
248    
249        fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
250        checkFinleyError();
251        AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
252        
253        /* win32 refactor */
254        TMPMEMFREE(fName);
255        
256        blocktimer_increment("ReadGmsh()", blocktimer_start);
257      return temp;      return temp;
258    }    }
259    
# Line 50  namespace finley { Line 262  namespace finley {
262              int periodic0,int periodic1,              int periodic0,int periodic1,
263              int periodic2,              int periodic2,
264              int integrationOrder,              int integrationOrder,
265              int useElementsOnFace)                      int reducedIntegrationOrder,
266                int useElementsOnFace,
267                int useFullElementOrder,
268                        int optimize)
269    {    {
 //     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;  
         
270      int numElements[]={n0,n1,n2};      int numElements[]={n0,n1,n2};
271      double length[]={l0,l1,l2};      double length[]={l0,l1,l2};
272      int periodic[]={periodic0, periodic1, periodic2};      int periodic[]={periodic0, periodic1, periodic2};
273    
274      //      //
275      // linearInterpolation      // linearInterpolation
276      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=NULL;
277    
278      if (order==1) {      if (order==1) {
279        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
280                      useElementsOnFace) ;                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
281      } else if (order==2) {      }
282        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,          else if (order==2) {
283                       useElementsOnFace) ;        fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
284                         useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
285      } else {      } else {
286        stringstream temp;        stringstream temp;
287        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
# Line 89  namespace finley { Line 297  namespace finley {
297              double l0, double l1,              double l0, double l1,
298              int periodic0,int periodic1,              int periodic0,int periodic1,
299              int integrationOrder,              int integrationOrder,
300              int useElementsOnFace)                          int reducedIntegrationOrder,
301                int useElementsOnFace,
302                    int useFullElementOrder,
303                            int optimize)
304    {    {
305      int numElements[]={n0,n1};      int numElements[]={n0,n1};
306      double length[]={l0,l1};      double length[]={l0,l1};
307      int periodic[]={periodic0, periodic1};      int periodic[]={periodic0, periodic1};
308    
309      Finley_Mesh* fMesh;      Finley_Mesh* fMesh=0;
310      if (order==1) {      if (order==1) {
311        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,        fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
312                      useElementsOnFace);                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
313      } else if (order==2) {      }
314        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,      else if (order==2) {
315                      useElementsOnFace);        fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
316      } else {                      useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
317        }
318        else {
319        stringstream temp;        stringstream temp;
320        temp << "Illegal interpolation order: " << order;        temp << "Illegal interpolation order: " << order;
321        setFinleyError(VALUE_ERROR,temp.str().c_str());        setFinleyError(VALUE_ERROR,temp.str().c_str());
# Line 113  namespace finley { Line 326  namespace finley {
326      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
327      return temp;      return temp;
328    }    }
329    AbstractContinuousDomain*  interval(int n0,int order,double l0,int periodic0,  
330                 int integrationOrder,    AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
331                 int useElementsOnFace)    {
332    {      Finley_Mesh* fMesh=0;
333      int numElements[]={n0};      //
334      double length[]={l0};      // extract the meshes from meshList
335      int periodic[]={periodic0};      int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
336      Finley_Mesh* fMesh;      Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
337      if (order==1) {      for (int i=0;i<numMsh;++i) {
338        fMesh=Finley_RectangularMesh_Line2(numElements, length,periodic,integrationOrder,           AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
339                       useElementsOnFace);           const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
340      } 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());  
341      }      }
342      //      //
343        // merge the meshes:
344        fMesh=Finley_Mesh_merge(numMsh,mshes);
345          TMPMEMFREE(mshes);
346        //
347      // Convert any finley errors into a C++ exception      // Convert any finley errors into a C++ exception
348      checkFinleyError();      checkFinleyError();
349      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);      AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
350      return temp;  
   }  
   AbstractContinuousDomain*  meshMerge(const boost::python::list& meshList)  
   {  
     AbstractContinuousDomain* temp=new MeshAdapter(0);  
351      return temp;      return temp;
352    }    }
353    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  glueFaces(const boost::python::list& meshList,
354              double safetyFactor,                                     double safety_factor,
355              double tolerance)                             double tolerance,
356                                           int optimize)
357    {    {
358      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
359      return temp;      //
360        // merge the meshes:
361        AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
362        //
363        // glue the faces:
364        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
365        fMesh=merged_finley_meshes->getFinley_Mesh();
366        Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
367    
368        //
369        // Convert any finley errors into a C++ exception
370        checkFinleyError();
371        return merged_meshes;
372    }    }
373    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,    AbstractContinuousDomain*  joinFaces(const boost::python::list& meshList,
374              double safety_factor,              double safety_factor,
375              double tolerance)              double tolerance,
376                            int optimize)
377    {    {
378      AbstractContinuousDomain* temp=new MeshAdapter(0);      Finley_Mesh* fMesh=0;
379      return temp;      //
380        // merge the meshes:
381        AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
382        //
383        // join the faces:
384        const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
385        fMesh=merged_finley_meshes->getFinley_Mesh();
386        Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
387        //
388        // Convert any finley errors into a C++ exception
389        checkFinleyError();
390        return merged_meshes;
391    }    }
392    
393  }  // end of namespace    // end of namespace
394    
395    }

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

  ViewVC Help
Powered by ViewVC 1.1.26