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

Annotation of /trunk/finley/src/CPPAdapter/MeshAdapterFactory.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1347 - (hide annotations)
Fri Nov 16 05:37:07 2007 UTC (11 years, 10 months ago) by ksteube
File size: 31498 byte(s)
Completed mesh.dump(file) and mesh=LoadMesh(file) by adding TagMap and
implementing MPI parallelism.
Now allocating ElementFile for ContactElements even if there are none.
Removed file Mesh_dump.c since dump/loadMesh are in CPPAdapter/MeshAdapter*.cpp.

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26