/[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 1634 - (hide annotations)
Sat Jul 12 09:08:33 2008 UTC (11 years, 2 months ago) by phornby
File size: 33100 byte(s)
linux_gcc_eg_options.py:
remove the std99 option, it is no longer needed as the code compiles without
C 1999 extension (need for these extensions elinminated in windows port).
Turn on all warnings except unknown pragmas. Should catch a lot of stuff.

SConstruct:
Impassioned plea

system_dep.h:
Add the standard incantation for dealing with const declarations
in C code called from C and C++

blocktimer:
Get the calling interface right for C code called from C and C++
and use __const as defined in system_dep.h
(Should be re-factored into compiler_dep.h file).

MeshAdapterFactory.cpp:
Since we have (effectively) no control over netCDF policy,
cast const char *'s to char *'s


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26