/[escript]/branches/more_shared_ptrs_from_1812/finley/src/CPPAdapter/MeshAdapterFactory.cpp
ViewVC logotype

Annotation of /branches/more_shared_ptrs_from_1812/finley/src/CPPAdapter/MeshAdapterFactory.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1776 - (hide annotations)
Tue Sep 9 06:03:53 2008 UTC (11 years, 6 months ago) by ksteube
Original Path: trunk/finley/src/CPPAdapter/MeshAdapterFactory.cpp
File size: 33350 byte(s)
Cleared away some debugging prints

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26