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

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

  ViewVC Help
Powered by ViewVC 1.1.26