/[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 1872 - (hide annotations)
Mon Oct 13 00:18:55 2008 UTC (10 years, 11 months ago) by jfenwick
File size: 32529 byte(s)
Closing the moreshared branch

1 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 ksteube 1312
14 ksteube 1811
15 ksteube 817 #ifdef PASO_MPI
16     #include <mpi.h>
17     #endif
18 ksteube 1345 #ifdef USE_NETCDF
19     #include <netcdfcpp.h>
20     #endif
21 jgs 203 #include "MeshAdapterFactory.h"
22 jgs 480 #include "FinleyError.h"
23 ksteube 1345 extern "C" {
24     #include "escript/blocktimer.h"
25     }
26 jgs 82
27 jgs 480 #include <boost/python/extract.hpp>
28    
29     #include <sstream>
30    
31 jgs 82 using namespace std;
32     using namespace escript;
33    
34     namespace finley {
35    
36 ksteube 1345 #ifdef USE_NETCDF
37     // A convenience method to retrieve an integer attribute from a NetCDF file
38     int NetCDF_Get_Int_Attribute(NcFile *dataFile, char *fName, char *attr_name) {
39     NcAtt *attr;
40     char error_msg[LenErrorMsg_MAX];
41     if (! (attr=dataFile->get_att(attr_name)) ) {
42     sprintf(error_msg,"Error retrieving integer attribute '%s' from NetCDF file '%s'", attr_name, fName);
43     throw DataException(error_msg);
44     }
45     int temp = attr->as_int(0);
46     delete attr;
47     return(temp);
48     }
49     #endif
50    
51 jfenwick 1872 // AbstractContinuousDomain* loadMesh(const std::string& fileName)
52     Domain_ptr 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 jfenwick 1872 return temp->getPtr();
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 jfenwick 1872 Domain_ptr 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 jfenwick 1872 return temp->getPtr();
504 ksteube 1360 }
505    
506 jfenwick 1872 Domain_ptr readGmsh(const std::string& fileName,
507 gross 934 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 jfenwick 1872 return temp->getPtr();
536 gross 934 }
537    
538 jfenwick 1872 /* AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,*/
539     Domain_ptr brick(int n0,int n1,int n2,int order,
540 jgs 82 double l0,double l1,double l2,
541     int periodic0,int periodic1,
542     int periodic2,
543     int integrationOrder,
544 gross 1059 int reducedIntegrationOrder,
545 ksteube 1312 int useElementsOnFace,
546     int useFullElementOrder,
547     int optimize)
548 jgs 82 {
549     int numElements[]={n0,n1,n2};
550     double length[]={l0,l1,l2};
551     int periodic[]={periodic0, periodic1, periodic2};
552    
553     //
554     // linearInterpolation
555 bcumming 751 Finley_Mesh* fMesh=NULL;
556    
557 jgs 82 if (order==1) {
558 gross 1062 fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
559 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
560 bcumming 751 }
561     else if (order==2) {
562 gross 1062 fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
563 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
564 jgs 82 } else {
565     stringstream temp;
566     temp << "Illegal interpolation order: " << order;
567     setFinleyError(VALUE_ERROR,temp.str().c_str());
568     }
569     //
570     // Convert any finley errors into a C++ exception
571     checkFinleyError();
572     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
573 jfenwick 1872 return temp->getPtr();
574 jgs 82 }
575 jfenwick 1872
576     /* AbstractContinuousDomain* rectangle(int n0,int n1,int order,*/
577     Domain_ptr rectangle(int n0,int n1,int order,
578 jgs 82 double l0, double l1,
579     int periodic0,int periodic1,
580     int integrationOrder,
581 gross 1059 int reducedIntegrationOrder,
582 ksteube 1312 int useElementsOnFace,
583     int useFullElementOrder,
584     int optimize)
585 jgs 82 {
586     int numElements[]={n0,n1};
587     double length[]={l0,l1};
588     int periodic[]={periodic0, periodic1};
589    
590 bcumming 751 Finley_Mesh* fMesh=0;
591 jgs 82 if (order==1) {
592 gross 1062 fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
593 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
594 bcumming 751 }
595     else if (order==2) {
596 gross 1062 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
597 ksteube 1312 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
598 bcumming 751 }
599     else {
600 jgs 82 stringstream temp;
601     temp << "Illegal interpolation order: " << order;
602     setFinleyError(VALUE_ERROR,temp.str().c_str());
603     }
604     //
605     // Convert any finley errors into a C++ exception
606     checkFinleyError();
607     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
608 jfenwick 1872 return temp->getPtr();
609 jgs 82 }
610 gross 1062
611 jfenwick 1872 Domain_ptr meshMerge(const boost::python::list& meshList)
612 jgs 82 {
613 jgs 110 Finley_Mesh* fMesh=0;
614     //
615     // extract the meshes from meshList
616     int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
617 woo409 757 Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
618 jgs 110 for (int i=0;i<numMsh;++i) {
619     AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
620     const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
621     mshes[i]=finley_meshListMember->getFinley_Mesh();
622     }
623     //
624     // merge the meshes:
625     fMesh=Finley_Mesh_merge(numMsh,mshes);
626 bcumming 759 TMPMEMFREE(mshes);
627 jgs 110 //
628     // Convert any finley errors into a C++ exception
629     checkFinleyError();
630     AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
631 woo409 757
632 jfenwick 1872 return temp->getPtr();
633 jgs 82 }
634 jfenwick 1872
635     Domain_ptr glueFaces(const boost::python::list& meshList,
636 ksteube 1312 double safety_factor,
637     double tolerance,
638     int optimize)
639 jgs 82 {
640 jgs 110 Finley_Mesh* fMesh=0;
641     //
642     // merge the meshes:
643 jfenwick 1872 Domain_ptr merged_meshes=meshMerge(meshList);
644    
645 jgs 110 //
646     // glue the faces:
647 jfenwick 1872 const MeshAdapter* merged_finley_meshes=dynamic_cast<const MeshAdapter*>(merged_meshes.get());
648 jgs 110 fMesh=merged_finley_meshes->getFinley_Mesh();
649 ksteube 1312 Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
650 bcumming 751
651 jgs 110 //
652     // Convert any finley errors into a C++ exception
653     checkFinleyError();
654 jfenwick 1872 return merged_meshes->getPtr();
655 jgs 82 }
656 jfenwick 1872 Domain_ptr joinFaces(const boost::python::list& meshList,
657 jgs 82 double safety_factor,
658 gross 1059 double tolerance,
659 ksteube 1312 int optimize)
660 jgs 82 {
661 jgs 110 Finley_Mesh* fMesh=0;
662     //
663     // merge the meshes:
664 jfenwick 1872 Domain_ptr merged_meshes=meshMerge(meshList);
665 jgs 110 //
666     // join the faces:
667 jfenwick 1872 const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes.get());
668 jgs 110 fMesh=merged_finley_meshes->getFinley_Mesh();
669 ksteube 1312 Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
670 jgs 110 //
671     // Convert any finley errors into a C++ exception
672     checkFinleyError();
673 jfenwick 1872 return merged_meshes->getPtr();
674 jgs 82 }
675    
676 gross 1062 // end of namespace
677    
678     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26