1 |
|
2 |
/* $Id$ */ |
3 |
|
4 |
/******************************************************* |
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 |
#ifdef PASO_MPI |
17 |
#include <mpi.h> |
18 |
#endif |
19 |
#ifdef USE_NETCDF |
20 |
#include <netcdfcpp.h> |
21 |
#endif |
22 |
#include "MeshAdapterFactory.h" |
23 |
#include "FinleyError.h" |
24 |
extern "C" { |
25 |
#include "escript/blocktimer.h" |
26 |
} |
27 |
|
28 |
#include <boost/python/extract.hpp> |
29 |
|
30 |
#include <sstream> |
31 |
|
32 |
using namespace std; |
33 |
using namespace escript; |
34 |
|
35 |
namespace finley { |
36 |
|
37 |
#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 |
{ |
54 |
#ifdef USE_NETCDF |
55 |
bool optimize=FALSE; // Don't optimize since this would cause problems with Data().dump() |
56 |
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 |
|
61 |
char *fName = Paso_MPI_appendRankToFileName(fileName.c_str(), |
62 |
mpi_info->size, |
63 |
mpi_info->rank); |
64 |
|
65 |
double blocktimer_start = blocktimer_time(); |
66 |
Finley_resetError(); |
67 |
int *first_DofComponent, *first_NodeComponent; |
68 |
|
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 |
throw DataException(error_msg); |
81 |
} |
82 |
|
83 |
// Read NetCDF integer attributes |
84 |
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 |
|
103 |
// 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 |
// 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 |
TMPMEMFREE(mesh_p->Nodes->Id); |
132 |
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 |
TMPMEMFREE(mesh_p->Nodes->Tag); |
139 |
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 |
TMPMEMFREE(mesh_p->Nodes->globalDegreesOfFreedom); |
146 |
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 |
TMPMEMFREE(mesh_p->Nodes->globalNodesIndex); |
153 |
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 |
TMPMEMFREE(mesh_p->Nodes->globalReducedDOFIndex); |
160 |
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 |
TMPMEMFREE(mesh_p->Nodes->globalReducedNodesIndex); |
167 |
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 |
TMPMEMFREE(mesh_p->Nodes->Coordinates); |
172 |
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 |
TMPMEMFREE(mesh_p->Nodes->Coordinates); |
176 |
throw DataException("Error - load:: unable to recover Nodes_Coordinates from netCDF file: " + *fName); |
177 |
} |
178 |
// Nodes_DofDistribution |
179 |
first_DofComponent = TMPMEMALLOC(mpi_size+1,index_t); |
180 |
if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) ) |
181 |
throw DataException("Error - loadMesh:: unable to read Nodes_DofDistribution from netCDF file: " + *fName); |
182 |
if (! nc_var_temp->get(&first_DofComponent[0], mpi_size+1) ) { |
183 |
throw DataException("Error - loadMesh:: unable to recover Nodes_DofDistribution from NetCDF file: " + *fName); |
184 |
} |
185 |
{ int i; printf("ksteube read Nodes_DofDistribution:"); for(i=0; i<mpi_size+1; i++) printf(" %d", first_DofComponent[i]); printf("\n"); } |
186 |
|
187 |
// 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 |
/* read elements */ |
197 |
if (Finley_noError()) { |
198 |
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 |
if (num_Elements>0) { |
203 |
if (Finley_noError()) { |
204 |
// 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 |
TMPMEMFREE(mesh_p->Elements->Id); |
209 |
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 |
TMPMEMFREE(mesh_p->Elements->Tag); |
216 |
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 |
TMPMEMFREE(mesh_p->Elements->Owner); |
223 |
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 |
TMPMEMFREE(mesh_p->Elements->Color); |
230 |
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 |
TMPMEMFREE(mesh_p->Elements->Nodes); |
236 |
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 |
TMPMEMFREE(Elements_Nodes); |
240 |
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 |
} |
251 |
} /* num_Elements>0 */ |
252 |
} |
253 |
|
254 |
/* get the face elements */ |
255 |
if (Finley_noError()) { |
256 |
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 |
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 |
TMPMEMFREE(mesh_p->FaceElements->Id); |
267 |
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 |
TMPMEMFREE(mesh_p->FaceElements->Tag); |
274 |
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 |
TMPMEMFREE(mesh_p->FaceElements->Owner); |
281 |
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 |
TMPMEMFREE(mesh_p->FaceElements->Color); |
288 |
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 |
TMPMEMFREE(mesh_p->FaceElements->Nodes); |
294 |
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 |
TMPMEMFREE(FaceElements_Nodes); |
298 |
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 |
|
312 |
/* get the Contact elements */ |
313 |
if (Finley_noError()) { |
314 |
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 |
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 |
TMPMEMFREE(mesh_p->ContactElements->Id); |
325 |
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 |
TMPMEMFREE(mesh_p->ContactElements->Tag); |
332 |
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 |
TMPMEMFREE(mesh_p->ContactElements->Owner); |
339 |
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 |
TMPMEMFREE(mesh_p->ContactElements->Color); |
346 |
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 |
TMPMEMFREE(mesh_p->ContactElements->Nodes); |
352 |
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 |
TMPMEMFREE(ContactElements_Nodes); |
356 |
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 |
|
370 |
/* get the Points (nodal elements) */ |
371 |
if (Finley_noError()) { |
372 |
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 |
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 |
TMPMEMFREE(mesh_p->Points->Id); |
383 |
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 |
TMPMEMFREE(mesh_p->Points->Tag); |
390 |
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 |
TMPMEMFREE(mesh_p->Points->Owner); |
397 |
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 |
TMPMEMFREE(mesh_p->Points->Color); |
404 |
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 |
TMPMEMFREE(mesh_p->Points->Nodes); |
410 |
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 |
TMPMEMFREE(Points_Nodes); |
414 |
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 |
|
425 |
/* 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 |
|
433 |
// 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 |
TMPMEMFREE(Tags_keys); |
438 |
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 |
} /* Finley_noError() after Finley_Mesh_alloc() */ |
455 |
|
456 |
if (Finley_noError()) Finley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent); |
457 |
TMPMEMFREE(first_DofComponent); |
458 |
TMPMEMFREE(first_NodeComponent); |
459 |
|
460 |
checkFinleyError(); |
461 |
temp=new MeshAdapter(mesh_p); |
462 |
|
463 |
if (! Finley_noError()) { |
464 |
Finley_Mesh_free(mesh_p); |
465 |
} |
466 |
|
467 |
/* win32 refactor */ |
468 |
TMPMEMFREE(fName); |
469 |
|
470 |
blocktimer_increment("LoadMesh()", blocktimer_start); |
471 |
return temp; |
472 |
#else |
473 |
throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager."); |
474 |
#endif /* USE_NETCDF */ |
475 |
} |
476 |
|
477 |
AbstractContinuousDomain* readMesh(const std::string& fileName, |
478 |
int integrationOrder, |
479 |
int reducedIntegrationOrder, |
480 |
int optimize) |
481 |
{ |
482 |
// |
483 |
// create a copy of the filename to overcome the non-constness of call |
484 |
// to Finley_Mesh_read |
485 |
Finley_Mesh* fMesh=0; |
486 |
// Win32 refactor |
487 |
if( fileName.size() == 0 ) |
488 |
{ |
489 |
throw DataException("Null file name!"); |
490 |
} |
491 |
|
492 |
char *fName = TMPMEMALLOC(fileName.size()+1,char); |
493 |
|
494 |
strcpy(fName,fileName.c_str()); |
495 |
double blocktimer_start = blocktimer_time(); |
496 |
|
497 |
fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE)); |
498 |
checkFinleyError(); |
499 |
AbstractContinuousDomain* temp=new MeshAdapter(fMesh); |
500 |
|
501 |
/* win32 refactor */ |
502 |
TMPMEMFREE(fName); |
503 |
|
504 |
blocktimer_increment("ReadMesh()", blocktimer_start); |
505 |
return temp; |
506 |
} |
507 |
|
508 |
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 |
if( fileName.size() == 0 ) |
519 |
{ |
520 |
throw DataException("Null file name!"); |
521 |
} |
522 |
|
523 |
char *fName = TMPMEMALLOC(fileName.size()+1,char); |
524 |
|
525 |
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 |
AbstractContinuousDomain* readGmsh(const std::string& fileName, |
540 |
int numDim, |
541 |
int integrationOrder, |
542 |
int reducedIntegrationOrder, |
543 |
int optimize) |
544 |
{ |
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 |
if( fileName.size() == 0 ) |
551 |
{ |
552 |
throw DataException("Null file name!"); |
553 |
} |
554 |
|
555 |
char *fName = TMPMEMALLOC(fileName.size()+1,char); |
556 |
|
557 |
strcpy(fName,fileName.c_str()); |
558 |
double blocktimer_start = blocktimer_time(); |
559 |
|
560 |
fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE)); |
561 |
checkFinleyError(); |
562 |
AbstractContinuousDomain* temp=new MeshAdapter(fMesh); |
563 |
|
564 |
/* win32 refactor */ |
565 |
TMPMEMFREE(fName); |
566 |
|
567 |
blocktimer_increment("ReadGmsh()", blocktimer_start); |
568 |
return temp; |
569 |
} |
570 |
|
571 |
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 |
int reducedIntegrationOrder, |
577 |
int useElementsOnFace, |
578 |
int useFullElementOrder, |
579 |
int optimize) |
580 |
{ |
581 |
int numElements[]={n0,n1,n2}; |
582 |
double length[]={l0,l1,l2}; |
583 |
int periodic[]={periodic0, periodic1, periodic2}; |
584 |
|
585 |
// |
586 |
// linearInterpolation |
587 |
Finley_Mesh* fMesh=NULL; |
588 |
|
589 |
if (order==1) { |
590 |
fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder, |
591 |
useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ; |
592 |
} |
593 |
else if (order==2) { |
594 |
fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder, |
595 |
useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ; |
596 |
} 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 |
int reducedIntegrationOrder, |
612 |
int useElementsOnFace, |
613 |
int useFullElementOrder, |
614 |
int optimize) |
615 |
{ |
616 |
int numElements[]={n0,n1}; |
617 |
double length[]={l0,l1}; |
618 |
int periodic[]={periodic0, periodic1}; |
619 |
|
620 |
Finley_Mesh* fMesh=0; |
621 |
if (order==1) { |
622 |
fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder, |
623 |
useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)); |
624 |
} |
625 |
else if (order==2) { |
626 |
fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder, |
627 |
useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)); |
628 |
} |
629 |
else { |
630 |
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 |
|
641 |
AbstractContinuousDomain* meshMerge(const boost::python::list& meshList) |
642 |
{ |
643 |
Finley_Mesh* fMesh=0; |
644 |
// |
645 |
// extract the meshes from meshList |
646 |
int numMsh=boost::python::extract<int>(meshList.attr("__len__")()); |
647 |
Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL; |
648 |
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 |
TMPMEMFREE(mshes); |
657 |
// |
658 |
// Convert any finley errors into a C++ exception |
659 |
checkFinleyError(); |
660 |
AbstractContinuousDomain* temp=new MeshAdapter(fMesh); |
661 |
|
662 |
return temp; |
663 |
} |
664 |
AbstractContinuousDomain* glueFaces(const boost::python::list& meshList, |
665 |
double safety_factor, |
666 |
double tolerance, |
667 |
int optimize) |
668 |
{ |
669 |
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 |
Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE)); |
678 |
|
679 |
// |
680 |
// Convert any finley errors into a C++ exception |
681 |
checkFinleyError(); |
682 |
return merged_meshes; |
683 |
} |
684 |
AbstractContinuousDomain* joinFaces(const boost::python::list& meshList, |
685 |
double safety_factor, |
686 |
double tolerance, |
687 |
int optimize) |
688 |
{ |
689 |
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 |
Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE)); |
698 |
// |
699 |
// Convert any finley errors into a C++ exception |
700 |
checkFinleyError(); |
701 |
return merged_meshes; |
702 |
} |
703 |
|
704 |
// end of namespace |
705 |
|
706 |
} |