/[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 1384 - (hide annotations)
Fri Jan 11 02:29:38 2008 UTC (11 years, 8 months ago) by phornby
Original Path: temp_trunk_copy/finley/src/CPPAdapter/MeshAdapterFactory.cpp
File size: 32390 byte(s)
Make a temp copy of the trunk before checking in the windows changes


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 ksteube 1360 AbstractContinuousDomain* readMeshMPI(const std::string& fileName,
490     int integrationOrder,
491     int reducedIntegrationOrder,
492     int optimize)
493     {
494     //
495     // create a copy of the filename to overcome the non-constness of call
496     // to Finley_Mesh_read
497     Finley_Mesh* fMesh=0;
498     // Win32 refactor
499     char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;
500     strcpy(fName,fileName.c_str());
501     double blocktimer_start = blocktimer_time();
502    
503     fMesh=Finley_Mesh_read_MPI(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
504     checkFinleyError();
505     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
506    
507     /* win32 refactor */
508     TMPMEMFREE(fName);
509    
510     blocktimer_increment("ReadMesh()", blocktimer_start);
511     return temp;
512     }
513    
514 gross 934 AbstractContinuousDomain* readGmsh(const std::string& fileName,
515     int numDim,
516     int integrationOrder,
517     int reducedIntegrationOrder,
518 ksteube 1312 int optimize)
519 gross 934 {
520     //
521     // create a copy of the filename to overcome the non-constness of call
522     // to Finley_Mesh_read
523     Finley_Mesh* fMesh=0;
524     // Win32 refactor
525     char *fName = ((fileName.size()+1)>0) ? TMPMEMALLOC((fileName.size()+1),char) : (char*)NULL;
526     strcpy(fName,fileName.c_str());
527 ksteube 1345 double blocktimer_start = blocktimer_time();
528 gross 934
529 ksteube 1312 fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
530 gross 934 checkFinleyError();
531     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
532    
533     /* win32 refactor */
534     TMPMEMFREE(fName);
535    
536 ksteube 1345 blocktimer_increment("ReadGmsh()", blocktimer_start);
537 gross 934 return temp;
538     }
539    
540 jgs 82 AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,
541     double l0,double l1,double l2,
542     int periodic0,int periodic1,
543     int periodic2,
544     int integrationOrder,
545 gross 1059 int reducedIntegrationOrder,
546 ksteube 1312 int useElementsOnFace,
547     int useFullElementOrder,
548     int optimize)
549 jgs 82 {
550     int numElements[]={n0,n1,n2};
551     double length[]={l0,l1,l2};
552     int periodic[]={periodic0, periodic1, periodic2};
553    
554     //
555     // linearInterpolation
556 bcumming 751 Finley_Mesh* fMesh=NULL;
557    
558 jgs 82 if (order==1) {
559 gross 1062 fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
560 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
561 bcumming 751 }
562     else if (order==2) {
563 gross 1062 fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
564 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
565 jgs 82 } else {
566     stringstream temp;
567     temp << "Illegal interpolation order: " << order;
568     setFinleyError(VALUE_ERROR,temp.str().c_str());
569     }
570     //
571     // Convert any finley errors into a C++ exception
572     checkFinleyError();
573     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
574     return temp;
575     }
576     AbstractContinuousDomain* rectangle(int n0,int n1,int order,
577     double l0, double l1,
578     int periodic0,int periodic1,
579     int integrationOrder,
580 gross 1059 int reducedIntegrationOrder,
581 ksteube 1312 int useElementsOnFace,
582     int useFullElementOrder,
583     int optimize)
584 jgs 82 {
585     int numElements[]={n0,n1};
586     double length[]={l0,l1};
587     int periodic[]={periodic0, periodic1};
588    
589 bcumming 751 Finley_Mesh* fMesh=0;
590 jgs 82 if (order==1) {
591 gross 1062 fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
592 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
593 bcumming 751 }
594     else if (order==2) {
595 gross 1062 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
596 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
597 bcumming 751 }
598     else {
599 jgs 82 stringstream temp;
600     temp << "Illegal interpolation order: " << order;
601     setFinleyError(VALUE_ERROR,temp.str().c_str());
602     }
603     //
604     // Convert any finley errors into a C++ exception
605     checkFinleyError();
606     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
607     return temp;
608     }
609 gross 1062
610 jgs 110 AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
611 jgs 82 {
612 jgs 110 Finley_Mesh* fMesh=0;
613     //
614     // extract the meshes from meshList
615     int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
616 woo409 757 Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
617 jgs 110 for (int i=0;i<numMsh;++i) {
618     AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
619     const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
620     mshes[i]=finley_meshListMember->getFinley_Mesh();
621     }
622     //
623     // merge the meshes:
624     fMesh=Finley_Mesh_merge(numMsh,mshes);
625 bcumming 759 TMPMEMFREE(mshes);
626 jgs 110 //
627     // Convert any finley errors into a C++ exception
628     checkFinleyError();
629     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
630 woo409 757
631 jgs 82 return temp;
632     }
633     AbstractContinuousDomain* glueFaces(const boost::python::list& meshList,
634 ksteube 1312 double safety_factor,
635     double tolerance,
636     int optimize)
637 jgs 82 {
638 jgs 110 Finley_Mesh* fMesh=0;
639     //
640     // merge the meshes:
641     AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
642     //
643     // glue the faces:
644     const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
645     fMesh=merged_finley_meshes->getFinley_Mesh();
646 ksteube 1312 Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
647 bcumming 751
648 jgs 110 //
649     // Convert any finley errors into a C++ exception
650     checkFinleyError();
651     return merged_meshes;
652 jgs 82 }
653     AbstractContinuousDomain* joinFaces(const boost::python::list& meshList,
654     double safety_factor,
655 gross 1059 double tolerance,
656 ksteube 1312 int optimize)
657 jgs 82 {
658 jgs 110 Finley_Mesh* fMesh=0;
659     //
660     // merge the meshes:
661     AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
662     //
663     // join the faces:
664     const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
665     fMesh=merged_finley_meshes->getFinley_Mesh();
666 ksteube 1312 Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
667 jgs 110 //
668     // Convert any finley errors into a C++ exception
669     checkFinleyError();
670     return merged_meshes;
671 jgs 82 }
672    
673 gross 1062 // end of namespace
674    
675     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26