/[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 1755 - (hide annotations)
Mon Sep 8 01:34:40 2008 UTC (11 years, 1 month ago) by ksteube
File size: 33630 byte(s)
Added new field to NetCDF dump of mesh

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 1755 { int i; printf("ksteube read Nodes_DofDistribution:"); for(i=0; i<mpi_size+1; i++) printf(" %d", first_DofComponent[i]); printf("\n"); }
186 ksteube 1345
187 ksteube 1755 // Nodes_NodeDistribution
188     first_NodeComponent = TMPMEMALLOC(mpi_size+1,index_t);
189     if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) )
190     throw DataException("Error - loadMesh:: unable to read Nodes_NodeDistribution from netCDF file: " + *fName);
191     if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
192     throw DataException("Error - loadMesh:: unable to recover Nodes_NodeDistribution from NetCDF file: " + *fName);
193     }
194     { int i; printf("ksteube read Nodes_NodeDistribution:"); for(i=0; i<mpi_size+1; i++) printf(" %d", first_NodeComponent[i]); printf("\n"); }
195    
196 ksteube 1345 /* read elements */
197     if (Finley_noError()) {
198 ksteube 1347 mesh_p->Elements=Finley_ElementFile_alloc((ElementTypeId)Elements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
199     if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
200     mesh_p->Elements->minColor=0;
201     mesh_p->Elements->maxColor=num_Elements-1;
202 ksteube 1346 if (num_Elements>0) {
203 ksteube 1345 if (Finley_noError()) {
204 ksteube 1346 // Elements_Id
205     if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
206     throw DataException("Error - loadMesh:: unable to read Elements_Id from netCDF file: " + *fName);
207     if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) ) {
208 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Id);
209 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Elements_Id from NetCDF file: " + *fName);
210     }
211     // Elements_Tag
212     if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
213     throw DataException("Error - loadMesh:: unable to read Elements_Tag from netCDF file: " + *fName);
214     if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) ) {
215 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Tag);
216 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Elements_Tag from NetCDF file: " + *fName);
217     }
218     // Elements_Owner
219     if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
220     throw DataException("Error - loadMesh:: unable to read Elements_Owner from netCDF file: " + *fName);
221     if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) ) {
222 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Owner);
223 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Elements_Owner from NetCDF file: " + *fName);
224     }
225     // Elements_Color
226     if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
227     throw DataException("Error - loadMesh:: unable to read Elements_Color from netCDF file: " + *fName);
228     if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) ) {
229 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Color);
230 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Elements_Color from NetCDF file: " + *fName);
231     }
232     // Elements_Nodes
233     int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
234     if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
235 ksteube 1755 TMPMEMFREE(mesh_p->Elements->Nodes);
236 ksteube 1346 throw DataException("Error - loadMesh:: unable to read Elements_Nodes from netCDF file: " + *fName);
237     }
238     if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
239 ksteube 1755 TMPMEMFREE(Elements_Nodes);
240 ksteube 1346 throw DataException("Error - load:: unable to recover Elements_Nodes from netCDF file: " + *fName);
241     }
242     // Copy temp array into mesh_p->Elements->Nodes
243     for (int i=0; i<num_Elements; i++) {
244     for (int j=0; j<num_Elements_numNodes; j++) {
245     mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
246     = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
247     }
248     }
249     TMPMEMFREE(Elements_Nodes);
250 ksteube 1345 }
251 ksteube 1346 } /* num_Elements>0 */
252 ksteube 1345 }
253    
254     /* get the face elements */
255 ksteube 1346 if (Finley_noError()) {
256 ksteube 1347 mesh_p->FaceElements=Finley_ElementFile_alloc((ElementTypeId)FaceElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
257     if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
258     mesh_p->FaceElements->minColor=0;
259     mesh_p->FaceElements->maxColor=num_FaceElements-1;
260 ksteube 1346 if (num_FaceElements>0) {
261     if (Finley_noError()) {
262     // FaceElements_Id
263     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
264     throw DataException("Error - loadMesh:: unable to read FaceElements_Id from netCDF file: " + *fName);
265     if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) ) {
266 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Id);
267 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover FaceElements_Id from NetCDF file: " + *fName);
268     }
269     // FaceElements_Tag
270     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
271     throw DataException("Error - loadMesh:: unable to read FaceElements_Tag from netCDF file: " + *fName);
272     if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) ) {
273 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Tag);
274 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover FaceElements_Tag from NetCDF file: " + *fName);
275     }
276     // FaceElements_Owner
277     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
278     throw DataException("Error - loadMesh:: unable to read FaceElements_Owner from netCDF file: " + *fName);
279     if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) ) {
280 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Owner);
281 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover FaceElements_Owner from NetCDF file: " + *fName);
282     }
283     // FaceElements_Color
284     if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
285     throw DataException("Error - loadMesh:: unable to read FaceElements_Color from netCDF file: " + *fName);
286     if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) ) {
287 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Color);
288 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover FaceElements_Color from NetCDF file: " + *fName);
289     }
290     // FaceElements_Nodes
291     int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
292     if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
293 ksteube 1755 TMPMEMFREE(mesh_p->FaceElements->Nodes);
294 ksteube 1346 throw DataException("Error - loadMesh:: unable to read FaceElements_Nodes from netCDF file: " + *fName);
295     }
296     if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
297 ksteube 1755 TMPMEMFREE(FaceElements_Nodes);
298 ksteube 1346 throw DataException("Error - load:: unable to recover FaceElements_Nodes from netCDF file: " + *fName);
299     }
300     // Copy temp array into mesh_p->FaceElements->Nodes
301     for (int i=0; i<num_FaceElements; i++) {
302     for (int j=0; j<num_FaceElements_numNodes; j++) {
303     mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]
304     = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
305     }
306     }
307     TMPMEMFREE(FaceElements_Nodes);
308     }
309     } /* num_FaceElements>0 */
310     }
311 ksteube 1345
312 ksteube 1346 /* get the Contact elements */
313     if (Finley_noError()) {
314 ksteube 1347 mesh_p->ContactElements=Finley_ElementFile_alloc((ElementTypeId)ContactElements_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
315     if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->ContactElements, num_ContactElements);
316     mesh_p->ContactElements->minColor=0;
317     mesh_p->ContactElements->maxColor=num_ContactElements-1;
318 ksteube 1346 if (num_ContactElements>0) {
319     if (Finley_noError()) {
320     // ContactElements_Id
321     if (! ( nc_var_temp = dataFile.get_var("ContactElements_Id")) )
322     throw DataException("Error - loadMesh:: unable to read ContactElements_Id from netCDF file: " + *fName);
323     if (! nc_var_temp->get(&mesh_p->ContactElements->Id[0], num_ContactElements) ) {
324 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Id);
325 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover ContactElements_Id from NetCDF file: " + *fName);
326     }
327     // ContactElements_Tag
328     if (! ( nc_var_temp = dataFile.get_var("ContactElements_Tag")) )
329     throw DataException("Error - loadMesh:: unable to read ContactElements_Tag from netCDF file: " + *fName);
330     if (! nc_var_temp->get(&mesh_p->ContactElements->Tag[0], num_ContactElements) ) {
331 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Tag);
332 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover ContactElements_Tag from NetCDF file: " + *fName);
333     }
334     // ContactElements_Owner
335     if (! ( nc_var_temp = dataFile.get_var("ContactElements_Owner")) )
336     throw DataException("Error - loadMesh:: unable to read ContactElements_Owner from netCDF file: " + *fName);
337     if (! nc_var_temp->get(&mesh_p->ContactElements->Owner[0], num_ContactElements) ) {
338 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Owner);
339 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover ContactElements_Owner from NetCDF file: " + *fName);
340     }
341     // ContactElements_Color
342     if (! ( nc_var_temp = dataFile.get_var("ContactElements_Color")) )
343     throw DataException("Error - loadMesh:: unable to read ContactElements_Color from netCDF file: " + *fName);
344     if (! nc_var_temp->get(&mesh_p->ContactElements->Color[0], num_ContactElements) ) {
345 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Color);
346 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover ContactElements_Color from NetCDF file: " + *fName);
347     }
348     // ContactElements_Nodes
349     int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
350     if (!(nc_var_temp = dataFile.get_var("ContactElements_Nodes"))) {
351 ksteube 1755 TMPMEMFREE(mesh_p->ContactElements->Nodes);
352 ksteube 1346 throw DataException("Error - loadMesh:: unable to read ContactElements_Nodes from netCDF file: " + *fName);
353     }
354     if (! nc_var_temp->get(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes) ) {
355 ksteube 1755 TMPMEMFREE(ContactElements_Nodes);
356 ksteube 1346 throw DataException("Error - load:: unable to recover ContactElements_Nodes from netCDF file: " + *fName);
357     }
358     // Copy temp array into mesh_p->ContactElements->Nodes
359     for (int i=0; i<num_ContactElements; i++) {
360     for (int j=0; j<num_ContactElements_numNodes; j++) {
361     mesh_p->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]
362     = ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
363     }
364     }
365     TMPMEMFREE(ContactElements_Nodes);
366     }
367     } /* num_ContactElements>0 */
368     }
369 ksteube 1345
370 ksteube 1346 /* get the Points (nodal elements) */
371     if (Finley_noError()) {
372 ksteube 1347 mesh_p->Points=Finley_ElementFile_alloc((ElementTypeId)Points_TypeId,mesh_p->order, mesh_p->reduced_order, mpi_info);
373     if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Points, num_Points);
374     mesh_p->Points->minColor=0;
375     mesh_p->Points->maxColor=num_Points-1;
376 ksteube 1346 if (num_Points>0) {
377     if (Finley_noError()) {
378     // Points_Id
379     if (! ( nc_var_temp = dataFile.get_var("Points_Id")) )
380     throw DataException("Error - loadMesh:: unable to read Points_Id from netCDF file: " + *fName);
381     if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points) ) {
382 ksteube 1755 TMPMEMFREE(mesh_p->Points->Id);
383 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Points_Id from NetCDF file: " + *fName);
384     }
385     // Points_Tag
386     if (! ( nc_var_temp = dataFile.get_var("Points_Tag")) )
387     throw DataException("Error - loadMesh:: unable to read Points_Tag from netCDF file: " + *fName);
388     if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points) ) {
389 ksteube 1755 TMPMEMFREE(mesh_p->Points->Tag);
390 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Points_Tag from NetCDF file: " + *fName);
391     }
392     // Points_Owner
393     if (! ( nc_var_temp = dataFile.get_var("Points_Owner")) )
394     throw DataException("Error - loadMesh:: unable to read Points_Owner from netCDF file: " + *fName);
395     if (! nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points) ) {
396 ksteube 1755 TMPMEMFREE(mesh_p->Points->Owner);
397 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Points_Owner from NetCDF file: " + *fName);
398     }
399     // Points_Color
400     if (! ( nc_var_temp = dataFile.get_var("Points_Color")) )
401     throw DataException("Error - loadMesh:: unable to read Points_Color from netCDF file: " + *fName);
402     if (! nc_var_temp->get(&mesh_p->Points->Color[0], num_Points) ) {
403 ksteube 1755 TMPMEMFREE(mesh_p->Points->Color);
404 ksteube 1346 throw DataException("Error - loadMesh:: unable to recover Points_Color from NetCDF file: " + *fName);
405     }
406     // Points_Nodes
407     int *Points_Nodes = TMPMEMALLOC(num_Points,int);
408     if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
409 ksteube 1755 TMPMEMFREE(mesh_p->Points->Nodes);
410 ksteube 1346 throw DataException("Error - loadMesh:: unable to read Points_Nodes from netCDF file: " + *fName);
411     }
412     if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
413 ksteube 1755 TMPMEMFREE(Points_Nodes);
414 ksteube 1346 throw DataException("Error - load:: unable to recover Points_Nodes from netCDF file: " + *fName);
415     }
416     // Copy temp array into mesh_p->Points->Nodes
417     for (int i=0; i<num_Points; i++) {
418     mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
419     }
420     TMPMEMFREE(Points_Nodes);
421     }
422     } /* num_Points>0 */
423     }
424 ksteube 1345
425 ksteube 1347 /* get the tags */
426     if (Finley_noError()) {
427     if (num_Tags>0) {
428     // Temp storage to gather node IDs
429     int *Tags_keys = TMPMEMALLOC(num_Tags, int);
430     char name_temp[4096];
431     int i;
432 ksteube 1345
433 ksteube 1347 // Tags_keys
434     if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) )
435     throw DataException("Error - loadMesh:: unable to read Tags_keys from netCDF file: " + *fName);
436     if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
437 ksteube 1755 TMPMEMFREE(Tags_keys);
438 ksteube 1347 throw DataException("Error - loadMesh:: unable to recover Tags_keys from NetCDF file: " + *fName);
439     }
440     for (i=0; i<num_Tags; i++) {
441     // Retrieve tag name
442     sprintf(name_temp, "Tags_name_%d", i);
443     if (! (attr=dataFile.get_att(name_temp)) ) {
444     sprintf(error_msg,"Error retrieving tag name from NetCDF file '%s'", fName);
445     throw DataException(error_msg);
446     }
447     char *name = attr->as_string(0);
448     delete attr;
449     Finley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
450     }
451     }
452     }
453    
454 ksteube 1345 } /* Finley_noError() after Finley_Mesh_alloc() */
455    
456 ksteube 1755 if (Finley_noError()) Finley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
457     TMPMEMFREE(first_DofComponent);
458     TMPMEMFREE(first_NodeComponent);
459 ksteube 1345
460 ksteube 1312 checkFinleyError();
461 ksteube 1345 temp=new MeshAdapter(mesh_p);
462    
463     if (! Finley_noError()) {
464     Finley_Mesh_free(mesh_p);
465     }
466    
467 ksteube 1312 /* win32 refactor */
468     TMPMEMFREE(fName);
469 ksteube 1345
470     blocktimer_increment("LoadMesh()", blocktimer_start);
471 ksteube 1312 return temp;
472 ksteube 1345 #else
473     throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
474     #endif /* USE_NETCDF */
475 ksteube 1312 }
476    
477 jgs 82 AbstractContinuousDomain* readMesh(const std::string& fileName,
478 gross 1059 int integrationOrder,
479     int reducedIntegrationOrder,
480 ksteube 1312 int optimize)
481 jgs 82 {
482     //
483     // create a copy of the filename to overcome the non-constness of call
484     // to Finley_Mesh_read
485 bcumming 751 Finley_Mesh* fMesh=0;
486 woo409 757 // Win32 refactor
487 phornby 1628 if( fileName.size() == 0 )
488     {
489     throw DataException("Null file name!");
490     }
491    
492     char *fName = TMPMEMALLOC(fileName.size()+1,char);
493    
494 jgs 82 strcpy(fName,fileName.c_str());
495 ksteube 1345 double blocktimer_start = blocktimer_time();
496 bcumming 751
497 ksteube 1312 fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
498 jgs 82 checkFinleyError();
499     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
500 woo409 757
501     /* win32 refactor */
502     TMPMEMFREE(fName);
503    
504 ksteube 1345 blocktimer_increment("ReadMesh()", blocktimer_start);
505 jgs 82 return temp;
506     }
507    
508 ksteube 1360 AbstractContinuousDomain* readMeshMPI(const std::string& fileName,
509     int integrationOrder,
510     int reducedIntegrationOrder,
511     int optimize)
512     {
513     //
514     // create a copy of the filename to overcome the non-constness of call
515     // to Finley_Mesh_read
516     Finley_Mesh* fMesh=0;
517     // Win32 refactor
518 phornby 1628 if( fileName.size() == 0 )
519     {
520     throw DataException("Null file name!");
521     }
522    
523     char *fName = TMPMEMALLOC(fileName.size()+1,char);
524    
525 ksteube 1360 strcpy(fName,fileName.c_str());
526     double blocktimer_start = blocktimer_time();
527    
528     fMesh=Finley_Mesh_read_MPI(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
529     checkFinleyError();
530     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
531    
532     /* win32 refactor */
533     TMPMEMFREE(fName);
534    
535     blocktimer_increment("ReadMesh()", blocktimer_start);
536     return temp;
537     }
538    
539 gross 934 AbstractContinuousDomain* readGmsh(const std::string& fileName,
540     int numDim,
541     int integrationOrder,
542     int reducedIntegrationOrder,
543 ksteube 1312 int optimize)
544 gross 934 {
545     //
546     // create a copy of the filename to overcome the non-constness of call
547     // to Finley_Mesh_read
548     Finley_Mesh* fMesh=0;
549     // Win32 refactor
550 phornby 1628 if( fileName.size() == 0 )
551     {
552     throw DataException("Null file name!");
553     }
554    
555     char *fName = TMPMEMALLOC(fileName.size()+1,char);
556    
557 gross 934 strcpy(fName,fileName.c_str());
558 ksteube 1345 double blocktimer_start = blocktimer_time();
559 gross 934
560 ksteube 1312 fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
561 gross 934 checkFinleyError();
562     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
563    
564     /* win32 refactor */
565     TMPMEMFREE(fName);
566    
567 ksteube 1345 blocktimer_increment("ReadGmsh()", blocktimer_start);
568 gross 934 return temp;
569     }
570    
571 jgs 82 AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,
572     double l0,double l1,double l2,
573     int periodic0,int periodic1,
574     int periodic2,
575     int integrationOrder,
576 gross 1059 int reducedIntegrationOrder,
577 ksteube 1312 int useElementsOnFace,
578     int useFullElementOrder,
579     int optimize)
580 jgs 82 {
581     int numElements[]={n0,n1,n2};
582     double length[]={l0,l1,l2};
583     int periodic[]={periodic0, periodic1, periodic2};
584    
585     //
586     // linearInterpolation
587 bcumming 751 Finley_Mesh* fMesh=NULL;
588    
589 jgs 82 if (order==1) {
590 gross 1062 fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
591 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
592 bcumming 751 }
593     else if (order==2) {
594 gross 1062 fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
595 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
596 jgs 82 } else {
597     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     AbstractContinuousDomain* rectangle(int n0,int n1,int order,
608     double l0, double l1,
609     int periodic0,int periodic1,
610     int integrationOrder,
611 gross 1059 int reducedIntegrationOrder,
612 ksteube 1312 int useElementsOnFace,
613     int useFullElementOrder,
614     int optimize)
615 jgs 82 {
616     int numElements[]={n0,n1};
617     double length[]={l0,l1};
618     int periodic[]={periodic0, periodic1};
619    
620 bcumming 751 Finley_Mesh* fMesh=0;
621 jgs 82 if (order==1) {
622 gross 1062 fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
623 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
624 bcumming 751 }
625     else if (order==2) {
626 gross 1062 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
627 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
628 bcumming 751 }
629     else {
630 jgs 82 stringstream temp;
631     temp << "Illegal interpolation order: " << order;
632     setFinleyError(VALUE_ERROR,temp.str().c_str());
633     }
634     //
635     // Convert any finley errors into a C++ exception
636     checkFinleyError();
637     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
638     return temp;
639     }
640 gross 1062
641 jgs 110 AbstractContinuousDomain* meshMerge(const boost::python::list& meshList)
642 jgs 82 {
643 jgs 110 Finley_Mesh* fMesh=0;
644     //
645     // extract the meshes from meshList
646     int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
647 woo409 757 Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
648 jgs 110 for (int i=0;i<numMsh;++i) {
649     AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
650     const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
651     mshes[i]=finley_meshListMember->getFinley_Mesh();
652     }
653     //
654     // merge the meshes:
655     fMesh=Finley_Mesh_merge(numMsh,mshes);
656 bcumming 759 TMPMEMFREE(mshes);
657 jgs 110 //
658     // Convert any finley errors into a C++ exception
659     checkFinleyError();
660     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
661 woo409 757
662 jgs 82 return temp;
663     }
664     AbstractContinuousDomain* glueFaces(const boost::python::list& meshList,
665 ksteube 1312 double safety_factor,
666     double tolerance,
667     int optimize)
668 jgs 82 {
669 jgs 110 Finley_Mesh* fMesh=0;
670     //
671     // merge the meshes:
672     AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
673     //
674     // glue the faces:
675     const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
676     fMesh=merged_finley_meshes->getFinley_Mesh();
677 ksteube 1312 Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
678 bcumming 751
679 jgs 110 //
680     // Convert any finley errors into a C++ exception
681     checkFinleyError();
682     return merged_meshes;
683 jgs 82 }
684     AbstractContinuousDomain* joinFaces(const boost::python::list& meshList,
685     double safety_factor,
686 gross 1059 double tolerance,
687 ksteube 1312 int optimize)
688 jgs 82 {
689 jgs 110 Finley_Mesh* fMesh=0;
690     //
691     // merge the meshes:
692     AbstractContinuousDomain* merged_meshes=meshMerge(meshList);
693     //
694     // join the faces:
695     const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes);
696     fMesh=merged_finley_meshes->getFinley_Mesh();
697 ksteube 1312 Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
698 jgs 110 //
699     // Convert any finley errors into a C++ exception
700     checkFinleyError();
701     return merged_meshes;
702 jgs 82 }
703    
704 gross 1062 // end of namespace
705    
706     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26