/[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 1807 - (hide annotations)
Thu Sep 25 01:04:51 2008 UTC (10 years, 11 months ago) by ksteube
File size: 32411 byte(s)
The new MPI parallel ReadMesh is now the default after having passed
all tests both with and without MPI.
If you have been using ReadMesh you should notice no difference.
If you were using ReadMeshMPI then change to ReadMesh.

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 char *fName = Paso_MPI_appendRankToFileName(fileName.c_str(),
62     mpi_info->size,
63     mpi_info->rank);
64    
65 ksteube 1345 double blocktimer_start = blocktimer_time();
66     Finley_resetError();
67 ksteube 1755 int *first_DofComponent, *first_NodeComponent;
68 ksteube 1345
69     // Open NetCDF file for reading
70     NcAtt *attr;
71     NcVar *nc_var_temp;
72     // netCDF error handler
73     NcError err(NcError::silent_nonfatal);
74     // Create the NetCDF file.
75     NcFile dataFile(fName, NcFile::ReadOnly);
76     if (!dataFile.is_valid()) {
77     sprintf(error_msg,"loadMesh: Opening file NetCDF %s for reading failed.", fName);
78     Finley_setError(IO_ERROR,error_msg);
79     Paso_MPIInfo_free( mpi_info );
80 ksteube 1347 throw DataException(error_msg);
81 ksteube 1345 }
82    
83     // Read NetCDF integer attributes
84 phornby 1634 int mpi_size = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_size");
85     int mpi_rank = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_rank");
86     int numDim = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numDim");
87     int order = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"order");
88     int reduced_order = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"reduced_order");
89     int numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numNodes");
90     int num_Elements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements");
91     int num_FaceElements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements");
92     int num_ContactElements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_ContactElements");
93     int num_Points = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Points");
94     int num_Elements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements_numNodes");
95     int Elements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Elements_TypeId");
96     int num_FaceElements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements_numNodes");
97     int FaceElements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"FaceElements_TypeId");
98     int num_ContactElements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_ContactElements_numNodes");
99     int ContactElements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"ContactElements_TypeId");
100     int Points_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Points_TypeId");
101     int num_Tags = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Tags");
102 ksteube 1345
103 ksteube 1347 // Verify size and rank
104     if (mpi_info->size != mpi_size) {
105     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);
106     throw DataException(error_msg);
107     }
108     if (mpi_info->rank != mpi_rank) {
109     sprintf(error_msg, "Error loadMesh - The NetCDF file '%s' should be read on CPU #%d instead of %d", fName, mpi_rank, mpi_info->rank);
110     throw DataException(error_msg);
111     }
112    
113 ksteube 1345 // Read mesh name
114     if (! (attr=dataFile.get_att("Name")) ) {
115     sprintf(error_msg,"Error retrieving mesh name from NetCDF file '%s'", fName);
116     throw DataException(error_msg);
117     }
118     char *name = attr->as_string(0);
119     delete attr;
120    
121     /* allocate mesh */
122     mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
123     if (Finley_noError()) {
124    
125     /* read nodes */
126     Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
127     // Nodes_Id
128     if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
129     throw DataException("Error - loadMesh:: unable to read Nodes_Id from netCDF file: " + *fName);
130     if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {
131 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->Id);
132 ksteube 1345 throw DataException("Error - loadMesh:: unable to recover Nodes_Id from NetCDF file: " + *fName);
133     }
134     // Nodes_Tag
135     if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
136     throw DataException("Error - loadMesh:: unable to read Nodes_Tag from netCDF file: " + *fName);
137     if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) ) {
138 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->Tag);
139 ksteube 1345 throw DataException("Error - loadMesh:: unable to recover Nodes_Tag from NetCDF file: " + *fName);
140     }
141     // Nodes_gDOF
142     if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
143     throw DataException("Error - loadMesh:: unable to read Nodes_gDOF from netCDF file: " + *fName);
144     if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) ) {
145 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->globalDegreesOfFreedom);
146 ksteube 1345 throw DataException("Error - loadMesh:: unable to recover Nodes_gDOF from NetCDF file: " + *fName);
147     }
148     // Nodes_gNI
149     if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
150     throw DataException("Error - loadMesh:: unable to read Nodes_gNI from netCDF file: " + *fName);
151     if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) ) {
152 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->globalNodesIndex);
153 ksteube 1345 throw DataException("Error - loadMesh:: unable to recover Nodes_gNI from NetCDF file: " + *fName);
154     }
155     // Nodes_grDfI
156     if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
157     throw DataException("Error - loadMesh:: unable to read Nodes_grDfI from netCDF file: " + *fName);
158     if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) ) {
159 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->globalReducedDOFIndex);
160 ksteube 1345 throw DataException("Error - loadMesh:: unable to recover Nodes_grDfI from NetCDF file: " + *fName);
161     }
162     // Nodes_grNI
163     if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
164     throw DataException("Error - loadMesh:: unable to read Nodes_grNI from netCDF file: " + *fName);
165     if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) ) {
166 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->globalReducedNodesIndex);
167 ksteube 1345 throw DataException("Error - loadMesh:: unable to recover Nodes_grNI from NetCDF file: " + *fName);
168     }
169     // Nodes_Coordinates
170     if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates"))) {
171 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->Coordinates);
172 ksteube 1345 throw DataException("Error - loadMesh:: unable to read Nodes_Coordinates from netCDF file: " + *fName);
173     }
174     if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) ) {
175 ksteube 1755 TMPMEMFREE(mesh_p->Nodes->Coordinates);
176 ksteube 1345 throw DataException("Error - load:: unable to recover Nodes_Coordinates from netCDF file: " + *fName);
177     }
178 ksteube 1347 // Nodes_DofDistribution
179 ksteube 1755 first_DofComponent = TMPMEMALLOC(mpi_size+1,index_t);
180 ksteube 1347 if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) )
181     throw DataException("Error - loadMesh:: unable to read Nodes_DofDistribution from netCDF file: " + *fName);
182 ksteube 1755 if (! nc_var_temp->get(&first_DofComponent[0], mpi_size+1) ) {
183 ksteube 1347 throw DataException("Error - loadMesh:: unable to recover Nodes_DofDistribution from NetCDF file: " + *fName);
184     }
185 ksteube 1345
186 ksteube 1755 // Nodes_NodeDistribution
187     first_NodeComponent = TMPMEMALLOC(mpi_size+1,index_t);
188     if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) )
189     throw DataException("Error - loadMesh:: unable to read Nodes_NodeDistribution from netCDF file: " + *fName);
190     if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
191     throw DataException("Error - loadMesh:: unable to recover Nodes_NodeDistribution from NetCDF file: " + *fName);
192     }
193    
194 ksteube 1345 /* 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 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Id);
207 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Tag);
214 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Owner);
221 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Color);
228 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Nodes);
234 ksteube 1346 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 ksteube 1755 TMPMEMFREE(Elements_Nodes);
238 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Id);
265 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Tag);
272 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Owner);
279 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Color);
286 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Nodes);
292 ksteube 1346 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 ksteube 1755 TMPMEMFREE(FaceElements_Nodes);
296 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Id);
323 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Tag);
330 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Owner);
337 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Color);
344 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Nodes);
350 ksteube 1346 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 ksteube 1755 TMPMEMFREE(ContactElements_Nodes);
354 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->Points->Id);
381 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->Points->Tag);
388 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->Points->Owner);
395 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->Points->Color);
402 ksteube 1346 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 ksteube 1755 TMPMEMFREE(mesh_p->Points->Nodes);
408 ksteube 1346 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 ksteube 1755 TMPMEMFREE(Points_Nodes);
412 ksteube 1346 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 ksteube 1755 TMPMEMFREE(Tags_keys);
436 ksteube 1347 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 1755 if (Finley_noError()) Finley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
455     TMPMEMFREE(first_DofComponent);
456     TMPMEMFREE(first_NodeComponent);
457 ksteube 1345
458 ksteube 1312 checkFinleyError();
459 ksteube 1345 temp=new MeshAdapter(mesh_p);
460    
461     if (! Finley_noError()) {
462     Finley_Mesh_free(mesh_p);
463     }
464    
465 ksteube 1312 /* win32 refactor */
466     TMPMEMFREE(fName);
467 ksteube 1345
468     blocktimer_increment("LoadMesh()", blocktimer_start);
469 ksteube 1312 return temp;
470 ksteube 1345 #else
471     throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
472     #endif /* USE_NETCDF */
473 ksteube 1312 }
474    
475 jgs 82 AbstractContinuousDomain* readMesh(const std::string& fileName,
476 gross 1059 int integrationOrder,
477     int reducedIntegrationOrder,
478 ksteube 1312 int optimize)
479 jgs 82 {
480     //
481     // create a copy of the filename to overcome the non-constness of call
482     // to Finley_Mesh_read
483 bcumming 751 Finley_Mesh* fMesh=0;
484 woo409 757 // Win32 refactor
485 phornby 1628 if( fileName.size() == 0 )
486     {
487     throw DataException("Null file name!");
488     }
489    
490     char *fName = TMPMEMALLOC(fileName.size()+1,char);
491    
492 jgs 82 strcpy(fName,fileName.c_str());
493 ksteube 1345 double blocktimer_start = blocktimer_time();
494 bcumming 751
495 ksteube 1360 fMesh=Finley_Mesh_read_MPI(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
496     checkFinleyError();
497     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
498    
499     /* win32 refactor */
500     TMPMEMFREE(fName);
501    
502     blocktimer_increment("ReadMesh()", blocktimer_start);
503     return temp;
504     }
505    
506 gross 934 AbstractContinuousDomain* readGmsh(const std::string& fileName,
507     int numDim,
508     int integrationOrder,
509     int reducedIntegrationOrder,
510 ksteube 1312 int optimize)
511 gross 934 {
512     //
513     // create a copy of the filename to overcome the non-constness of call
514     // to Finley_Mesh_read
515     Finley_Mesh* fMesh=0;
516     // Win32 refactor
517 phornby 1628 if( fileName.size() == 0 )
518     {
519     throw DataException("Null file name!");
520     }
521    
522     char *fName = TMPMEMALLOC(fileName.size()+1,char);
523    
524 gross 934 strcpy(fName,fileName.c_str());
525 ksteube 1345 double blocktimer_start = blocktimer_time();
526 gross 934
527 ksteube 1312 fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
528 gross 934 checkFinleyError();
529     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
530    
531     /* win32 refactor */
532     TMPMEMFREE(fName);
533    
534 ksteube 1345 blocktimer_increment("ReadGmsh()", blocktimer_start);
535 gross 934 return temp;
536     }
537    
538 jgs 82 AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,
539     double l0,double l1,double l2,
540     int periodic0,int periodic1,
541     int periodic2,
542     int integrationOrder,
543 gross 1059 int reducedIntegrationOrder,
544 ksteube 1312 int useElementsOnFace,
545     int useFullElementOrder,
546     int optimize)
547 jgs 82 {
548     int numElements[]={n0,n1,n2};
549     double length[]={l0,l1,l2};
550     int periodic[]={periodic0, periodic1, periodic2};
551    
552     //
553     // linearInterpolation
554 bcumming 751 Finley_Mesh* fMesh=NULL;
555    
556 jgs 82 if (order==1) {
557 gross 1062 fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
558 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
559 bcumming 751 }
560     else if (order==2) {
561 gross 1062 fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
562 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
563 jgs 82 } else {
564     stringstream temp;
565     temp << "Illegal interpolation order: " << order;
566     setFinleyError(VALUE_ERROR,temp.str().c_str());
567     }
568     //
569     // Convert any finley errors into a C++ exception
570     checkFinleyError();
571     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
572     return temp;
573     }
574     AbstractContinuousDomain* rectangle(int n0,int n1,int order,
575     double l0, double l1,
576     int periodic0,int periodic1,
577     int integrationOrder,
578 gross 1059 int reducedIntegrationOrder,
579 ksteube 1312 int useElementsOnFace,
580     int useFullElementOrder,
581     int optimize)
582 jgs 82 {
583     int numElements[]={n0,n1};
584     double length[]={l0,l1};
585     int periodic[]={periodic0, periodic1};
586    
587 bcumming 751 Finley_Mesh* fMesh=0;
588 jgs 82 if (order==1) {
589 gross 1062 fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
590 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
591 bcumming 751 }
592     else if (order==2) {
593 gross 1062 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
594 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
595 bcumming 751 }
596     else {
597 jgs 82 stringstream temp;
598     temp << "Illegal interpolation order: " << order;
599     setFinleyError(VALUE_ERROR,temp.str().c_str());
600     }
601     //
602     // Convert any finley errors into a C++ exception
603     checkFinleyError();
604     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
605     return temp;
606     }
607 gross 1062
608 jgs 110 AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
609 jgs 82 {
610 jgs 110 Finley_Mesh* fMesh=0;
611     //
612     // extract the meshes from meshList
613     int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
614 woo409 757 Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
615 jgs 110 for (int i=0;i<numMsh;++i) {
616     AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
617     const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
618     mshes[i]=finley_meshListMember->getFinley_Mesh();
619     }
620     //
621     // merge the meshes:
622     fMesh=Finley_Mesh_merge(numMsh,mshes);
623 bcumming 759 TMPMEMFREE(mshes);
624 jgs 110 //
625     // Convert any finley errors into a C++ exception
626     checkFinleyError();
627     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
628 woo409 757
629 jgs 82 return temp;
630     }
631     AbstractContinuousDomain* glueFaces(const boost::python::list& meshList,
632 ksteube 1312 double safety_factor,
633     double tolerance,
634     int optimize)
635 jgs 82 {
636 jgs 110 Finley_Mesh* fMesh=0;
637     //
638     // merge the meshes:
639     AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
640     //
641     // glue the faces:
642     const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
643     fMesh=merged_finley_meshes->getFinley_Mesh();
644 ksteube 1312 Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
645 bcumming 751
646 jgs 110 //
647     // Convert any finley errors into a C++ exception
648     checkFinleyError();
649     return merged_meshes;
650 jgs 82 }
651     AbstractContinuousDomain* joinFaces(const boost::python::list& meshList,
652     double safety_factor,
653 gross 1059 double tolerance,
654 ksteube 1312 int optimize)
655 jgs 82 {
656 jgs 110 Finley_Mesh* fMesh=0;
657     //
658     // merge the meshes:
659     AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
660     //
661     // join the faces:
662     const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
663     fMesh=merged_finley_meshes->getFinley_Mesh();
664 ksteube 1312 Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
665 jgs 110 //
666     // Convert any finley errors into a C++ exception
667     checkFinleyError();
668     return merged_meshes;
669 jgs 82 }
670    
671 gross 1062 // end of namespace
672    
673     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26