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

Contents of /trunk/finley/src/CPPAdapter/MeshAdapterFactory.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2748 - (show annotations)
Tue Nov 17 07:32:59 2009 UTC (9 years, 10 months ago) by gross
File size: 33657 byte(s)
Macro elements are implemented now. VTK writer for macro elements still needs testing.
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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
14
15 #ifdef PASO_MPI
16 #include <mpi.h>
17 #endif
18 #ifdef USE_NETCDF
19 #include <netcdfcpp.h>
20 #endif
21 #include "MeshAdapterFactory.h"
22 #include "FinleyError.h"
23 extern "C" {
24 #include "esysUtils/blocktimer.h"
25 }
26
27 #include <boost/python/extract.hpp>
28
29 #include <sstream>
30
31 using namespace std;
32 using namespace escript;
33
34 namespace finley {
35
36 #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 // AbstractContinuousDomain* loadMesh(const std::string& fileName)
52 Domain_ptr loadMesh(const std::string& fileName)
53 {
54 #ifdef USE_NETCDF
55 Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
56 AbstractContinuousDomain* temp;
57 Finley_Mesh *mesh_p=NULL;
58 char error_msg[LenErrorMsg_MAX];
59
60 char *fName = Paso_MPI_appendRankToFileName(fileName.c_str(),
61 mpi_info->size,
62 mpi_info->rank);
63
64 double blocktimer_start = blocktimer_time();
65 Finley_resetError();
66 int *first_DofComponent, *first_NodeComponent;
67
68 // Open NetCDF file for reading
69 NcAtt *attr;
70 NcVar *nc_var_temp;
71 // netCDF error handler
72 NcError err(NcError::silent_nonfatal);
73 // Create the NetCDF file.
74 NcFile dataFile(fName, NcFile::ReadOnly);
75 if (!dataFile.is_valid()) {
76 sprintf(error_msg,"loadMesh: Opening file NetCDF %s for reading failed.", fName);
77 Finley_setError(IO_ERROR,error_msg);
78 Paso_MPIInfo_free( mpi_info );
79 throw DataException(error_msg);
80 }
81
82 // Read NetCDF integer attributes
83 int mpi_size = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_size");
84 int mpi_rank = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_rank");
85 int numDim = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numDim");
86 int order = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"order");
87 int reduced_order = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"reduced_order");
88 int numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numNodes");
89 int num_Elements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements");
90 int num_FaceElements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements");
91 int num_ContactElements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_ContactElements");
92 int num_Points = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Points");
93 int num_Elements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements_numNodes");
94 int Elements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Elements_TypeId");
95 int num_FaceElements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements_numNodes");
96 int FaceElements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"FaceElements_TypeId");
97 int num_ContactElements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_ContactElements_numNodes");
98 int ContactElements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"ContactElements_TypeId");
99 int Points_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Points_TypeId");
100 int num_Tags = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Tags");
101
102 // Verify size and rank
103 if (mpi_info->size != mpi_size) {
104 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);
105 throw DataException(error_msg);
106 }
107 if (mpi_info->rank != mpi_rank) {
108 sprintf(error_msg, "Error loadMesh - The NetCDF file '%s' should be read on CPU #%d instead of %d", fName, mpi_rank, mpi_info->rank);
109 throw DataException(error_msg);
110 }
111
112 // Read mesh name
113 if (! (attr=dataFile.get_att("Name")) ) {
114 sprintf(error_msg,"Error retrieving mesh name from NetCDF file '%s'", fName);
115 throw DataException(error_msg);
116 }
117 char *name = attr->as_string(0);
118 delete attr;
119
120 /* allocate mesh */
121 mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
122 if (Finley_noError()) {
123
124 /* read nodes */
125 Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
126 // Nodes_Id
127 if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
128 throw DataException("Error - loadMesh:: unable to read Nodes_Id from netCDF file: " + *fName);
129 if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) ) {
130 TMPMEMFREE(mesh_p->Nodes->Id);
131 throw DataException("Error - loadMesh:: unable to recover Nodes_Id from NetCDF file: " + *fName);
132 }
133 // Nodes_Tag
134 if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
135 throw DataException("Error - loadMesh:: unable to read Nodes_Tag from netCDF file: " + *fName);
136 if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) ) {
137 TMPMEMFREE(mesh_p->Nodes->Tag);
138 throw DataException("Error - loadMesh:: unable to recover Nodes_Tag from NetCDF file: " + *fName);
139 }
140 // Nodes_gDOF
141 if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
142 throw DataException("Error - loadMesh:: unable to read Nodes_gDOF from netCDF file: " + *fName);
143 if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) ) {
144 TMPMEMFREE(mesh_p->Nodes->globalDegreesOfFreedom);
145 throw DataException("Error - loadMesh:: unable to recover Nodes_gDOF from NetCDF file: " + *fName);
146 }
147 // Nodes_gNI
148 if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
149 throw DataException("Error - loadMesh:: unable to read Nodes_gNI from netCDF file: " + *fName);
150 if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) ) {
151 TMPMEMFREE(mesh_p->Nodes->globalNodesIndex);
152 throw DataException("Error - loadMesh:: unable to recover Nodes_gNI from NetCDF file: " + *fName);
153 }
154 // Nodes_grDfI
155 if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
156 throw DataException("Error - loadMesh:: unable to read Nodes_grDfI from netCDF file: " + *fName);
157 if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) ) {
158 TMPMEMFREE(mesh_p->Nodes->globalReducedDOFIndex);
159 throw DataException("Error - loadMesh:: unable to recover Nodes_grDfI from NetCDF file: " + *fName);
160 }
161 // Nodes_grNI
162 if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
163 throw DataException("Error - loadMesh:: unable to read Nodes_grNI from netCDF file: " + *fName);
164 if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) ) {
165 TMPMEMFREE(mesh_p->Nodes->globalReducedNodesIndex);
166 throw DataException("Error - loadMesh:: unable to recover Nodes_grNI from NetCDF file: " + *fName);
167 }
168 // Nodes_Coordinates
169 if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates"))) {
170 TMPMEMFREE(mesh_p->Nodes->Coordinates);
171 throw DataException("Error - loadMesh:: unable to read Nodes_Coordinates from netCDF file: " + *fName);
172 }
173 if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) ) {
174 TMPMEMFREE(mesh_p->Nodes->Coordinates);
175 throw DataException("Error - load:: unable to recover Nodes_Coordinates from netCDF file: " + *fName);
176 }
177 // Nodes_DofDistribution
178 first_DofComponent = TMPMEMALLOC(mpi_size+1,index_t);
179 if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) )
180 throw DataException("Error - loadMesh:: unable to read Nodes_DofDistribution from netCDF file: " + *fName);
181 if (! nc_var_temp->get(&first_DofComponent[0], mpi_size+1) ) {
182 throw DataException("Error - loadMesh:: unable to recover Nodes_DofDistribution from NetCDF file: " + *fName);
183 }
184
185 // Nodes_NodeDistribution
186 first_NodeComponent = TMPMEMALLOC(mpi_size+1,index_t);
187 if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) )
188 throw DataException("Error - loadMesh:: unable to read Nodes_NodeDistribution from netCDF file: " + *fName);
189 if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
190 throw DataException("Error - loadMesh:: unable to recover Nodes_NodeDistribution from NetCDF file: " + *fName);
191 }
192
193 /* read elements */
194 if (Finley_noError()) {
195 Finley_ReferenceElementSet *refElements= Finley_ReferenceElementSet_alloc((ElementTypeId)Elements_TypeId,mesh_p->order, mesh_p->reduced_order);
196 if (Finley_noError()) {
197 mesh_p->Elements=Finley_ElementFile_alloc(refElements, mpi_info);
198 }
199 Finley_ReferenceElementSet_dealloc(refElements);
200 if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
201 if (Finley_noError()) {
202 mesh_p->Elements->minColor=0;
203 mesh_p->Elements->maxColor=num_Elements-1;
204 if (num_Elements>0) {
205 // Elements_Id
206 if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
207 throw DataException("Error - loadMesh:: unable to read Elements_Id from netCDF file: " + *fName);
208 if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) ) {
209 TMPMEMFREE(mesh_p->Elements->Id);
210 throw DataException("Error - loadMesh:: unable to recover Elements_Id from NetCDF file: " + *fName);
211 }
212 // Elements_Tag
213 if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
214 throw DataException("Error - loadMesh:: unable to read Elements_Tag from netCDF file: " + *fName);
215 if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) ) {
216 TMPMEMFREE(mesh_p->Elements->Tag);
217 throw DataException("Error - loadMesh:: unable to recover Elements_Tag from NetCDF file: " + *fName);
218 }
219 // Elements_Owner
220 if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
221 throw DataException("Error - loadMesh:: unable to read Elements_Owner from netCDF file: " + *fName);
222 if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) ) {
223 TMPMEMFREE(mesh_p->Elements->Owner);
224 throw DataException("Error - loadMesh:: unable to recover Elements_Owner from NetCDF file: " + *fName);
225 }
226 // Elements_Color
227 if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
228 throw DataException("Error - loadMesh:: unable to read Elements_Color from netCDF file: " + *fName);
229 if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) ) {
230 TMPMEMFREE(mesh_p->Elements->Color);
231 throw DataException("Error - loadMesh:: unable to recover Elements_Color from NetCDF file: " + *fName);
232 }
233 // Elements_Nodes
234 int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
235 if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
236 TMPMEMFREE(mesh_p->Elements->Nodes);
237 throw DataException("Error - loadMesh:: unable to read Elements_Nodes from netCDF file: " + *fName);
238 }
239 if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
240 TMPMEMFREE(Elements_Nodes);
241 throw DataException("Error - load:: unable to recover Elements_Nodes from netCDF file: " + *fName);
242 }
243 // Copy temp array into mesh_p->Elements->Nodes
244 for (int i=0; i<num_Elements; i++) {
245 for (int j=0; j<num_Elements_numNodes; j++) {
246 mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
247 = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
248 }
249 }
250 TMPMEMFREE(Elements_Nodes);
251
252 } /* num_Elements>0 */
253 }
254 }
255
256 /* get the face elements */
257 if (Finley_noError()) {
258 Finley_ReferenceElementSet *refFaceElements= Finley_ReferenceElementSet_alloc((ElementTypeId)FaceElements_TypeId ,mesh_p->order, mesh_p->reduced_order);
259 if (Finley_noError()) {
260 mesh_p->FaceElements=Finley_ElementFile_alloc(refFaceElements, mpi_info);
261 }
262 Finley_ReferenceElementSet_dealloc(refFaceElements);
263 if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
264 if (Finley_noError()) {
265 mesh_p->FaceElements->minColor=0;
266 mesh_p->FaceElements->maxColor=num_FaceElements-1;
267 if (num_FaceElements>0) {
268 // FaceElements_Id
269 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
270 throw DataException("Error - loadMesh:: unable to read FaceElements_Id from netCDF file: " + *fName);
271 if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) ) {
272 TMPMEMFREE(mesh_p->FaceElements->Id);
273 throw DataException("Error - loadMesh:: unable to recover FaceElements_Id from NetCDF file: " + *fName);
274 }
275 // FaceElements_Tag
276 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
277 throw DataException("Error - loadMesh:: unable to read FaceElements_Tag from netCDF file: " + *fName);
278 if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) ) {
279 TMPMEMFREE(mesh_p->FaceElements->Tag);
280 throw DataException("Error - loadMesh:: unable to recover FaceElements_Tag from NetCDF file: " + *fName);
281 }
282 // FaceElements_Owner
283 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
284 throw DataException("Error - loadMesh:: unable to read FaceElements_Owner from netCDF file: " + *fName);
285 if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) ) {
286 TMPMEMFREE(mesh_p->FaceElements->Owner);
287 throw DataException("Error - loadMesh:: unable to recover FaceElements_Owner from NetCDF file: " + *fName);
288 }
289 // FaceElements_Color
290 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
291 throw DataException("Error - loadMesh:: unable to read FaceElements_Color from netCDF file: " + *fName);
292 if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) ) {
293 TMPMEMFREE(mesh_p->FaceElements->Color);
294 throw DataException("Error - loadMesh:: unable to recover FaceElements_Color from NetCDF file: " + *fName);
295 }
296 // FaceElements_Nodes
297 int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
298 if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
299 TMPMEMFREE(mesh_p->FaceElements->Nodes);
300 throw DataException("Error - loadMesh:: unable to read FaceElements_Nodes from netCDF file: " + *fName);
301 }
302 if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
303 TMPMEMFREE(FaceElements_Nodes);
304 throw DataException("Error - load:: unable to recover FaceElements_Nodes from netCDF file: " + *fName);
305 }
306 // Copy temp array into mesh_p->FaceElements->Nodes
307 for (int i=0; i<num_FaceElements; i++) {
308 for (int j=0; j<num_FaceElements_numNodes; j++) {
309 mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
310 }
311 }
312 TMPMEMFREE(FaceElements_Nodes);
313 } /* num_FaceElements>0 */
314 }
315 }
316
317 /* get the Contact elements */
318 if (Finley_noError()) {
319 Finley_ReferenceElementSet *refContactElements= Finley_ReferenceElementSet_alloc((ElementTypeId)ContactElements_TypeId,mesh_p->order, mesh_p->reduced_order);
320 if (Finley_noError()) {
321 mesh_p->ContactElements=Finley_ElementFile_alloc(refContactElements, mpi_info);
322 }
323 Finley_ReferenceElementSet_dealloc(refContactElements);
324 if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->ContactElements, num_ContactElements);
325 if (Finley_noError()) {
326 mesh_p->ContactElements->minColor=0;
327 mesh_p->ContactElements->maxColor=num_ContactElements-1;
328 if (num_ContactElements>0) {
329 // ContactElements_Id
330 if (! ( nc_var_temp = dataFile.get_var("ContactElements_Id")) )
331 throw DataException("Error - loadMesh:: unable to read ContactElements_Id from netCDF file: " + *fName);
332 if (! nc_var_temp->get(&mesh_p->ContactElements->Id[0], num_ContactElements) ) {
333 TMPMEMFREE(mesh_p->ContactElements->Id);
334 throw DataException("Error - loadMesh:: unable to recover ContactElements_Id from NetCDF file: " + *fName);
335 }
336 // ContactElements_Tag
337 if (! ( nc_var_temp = dataFile.get_var("ContactElements_Tag")) )
338 throw DataException("Error - loadMesh:: unable to read ContactElements_Tag from netCDF file: " + *fName);
339 if (! nc_var_temp->get(&mesh_p->ContactElements->Tag[0], num_ContactElements) ) {
340 TMPMEMFREE(mesh_p->ContactElements->Tag);
341 throw DataException("Error - loadMesh:: unable to recover ContactElements_Tag from NetCDF file: " + *fName);
342 }
343 // ContactElements_Owner
344 if (! ( nc_var_temp = dataFile.get_var("ContactElements_Owner")) )
345 throw DataException("Error - loadMesh:: unable to read ContactElements_Owner from netCDF file: " + *fName);
346 if (! nc_var_temp->get(&mesh_p->ContactElements->Owner[0], num_ContactElements) ) {
347 TMPMEMFREE(mesh_p->ContactElements->Owner);
348 throw DataException("Error - loadMesh:: unable to recover ContactElements_Owner from NetCDF file: " + *fName);
349 }
350 // ContactElements_Color
351 if (! ( nc_var_temp = dataFile.get_var("ContactElements_Color")) )
352 throw DataException("Error - loadMesh:: unable to read ContactElements_Color from netCDF file: " + *fName);
353 if (! nc_var_temp->get(&mesh_p->ContactElements->Color[0], num_ContactElements) ) {
354 TMPMEMFREE(mesh_p->ContactElements->Color);
355 throw DataException("Error - loadMesh:: unable to recover ContactElements_Color from NetCDF file: " + *fName);
356 }
357 // ContactElements_Nodes
358 int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
359 if (!(nc_var_temp = dataFile.get_var("ContactElements_Nodes"))) {
360 TMPMEMFREE(mesh_p->ContactElements->Nodes);
361 throw DataException("Error - loadMesh:: unable to read ContactElements_Nodes from netCDF file: " + *fName);
362 }
363 if (! nc_var_temp->get(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes) ) {
364 TMPMEMFREE(ContactElements_Nodes);
365 throw DataException("Error - load:: unable to recover ContactElements_Nodes from netCDF file: " + *fName);
366 }
367 // Copy temp array into mesh_p->ContactElements->Nodes
368 for (int i=0; i<num_ContactElements; i++) {
369 for (int j=0; j<num_ContactElements_numNodes; j++) {
370 mesh_p->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]= ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
371 }
372 }
373 TMPMEMFREE(ContactElements_Nodes);
374 } /* num_ContactElements>0 */
375 }
376
377 }
378
379 /* get the Points (nodal elements) */
380 if (Finley_noError()) {
381 Finley_ReferenceElementSet *refPoints= Finley_ReferenceElementSet_alloc((ElementTypeId)Points_TypeId,mesh_p->order, mesh_p->reduced_order);
382 if (Finley_noError()) {
383 mesh_p->Points=Finley_ElementFile_alloc(refPoints, mpi_info);
384 }
385 Finley_ReferenceElementSet_dealloc(refPoints);
386 if (Finley_noError()) Finley_ElementFile_allocTable(mesh_p->Points, num_Points);
387 if (Finley_noError()) {
388 mesh_p->Points->minColor=0;
389 mesh_p->Points->maxColor=num_Points-1;
390 if (num_Points>0) {
391 // Points_Id
392 if (! ( nc_var_temp = dataFile.get_var("Points_Id")) )
393 throw DataException("Error - loadMesh:: unable to read Points_Id from netCDF file: " + *fName);
394 if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points) ) {
395 TMPMEMFREE(mesh_p->Points->Id);
396 throw DataException("Error - loadMesh:: unable to recover Points_Id from NetCDF file: " + *fName);
397 }
398 // Points_Tag
399 if (! ( nc_var_temp = dataFile.get_var("Points_Tag")) )
400 throw DataException("Error - loadMesh:: unable to read Points_Tag from netCDF file: " + *fName);
401 if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points) ) {
402 TMPMEMFREE(mesh_p->Points->Tag);
403 throw DataException("Error - loadMesh:: unable to recover Points_Tag from NetCDF file: " + *fName);
404 }
405 // Points_Owner
406 if (! ( nc_var_temp = dataFile.get_var("Points_Owner")) )
407 throw DataException("Error - loadMesh:: unable to read Points_Owner from netCDF file: " + *fName);
408 if (! nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points) ) {
409 TMPMEMFREE(mesh_p->Points->Owner);
410 throw DataException("Error - loadMesh:: unable to recover Points_Owner from NetCDF file: " + *fName);
411 }
412 // Points_Color
413 if (! ( nc_var_temp = dataFile.get_var("Points_Color")) )
414 throw DataException("Error - loadMesh:: unable to read Points_Color from netCDF file: " + *fName);
415 if (! nc_var_temp->get(&mesh_p->Points->Color[0], num_Points) ) {
416 TMPMEMFREE(mesh_p->Points->Color);
417 throw DataException("Error - loadMesh:: unable to recover Points_Color from NetCDF file: " + *fName);
418 }
419 // Points_Nodes
420 int *Points_Nodes = TMPMEMALLOC(num_Points,int);
421 if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
422 TMPMEMFREE(mesh_p->Points->Nodes);
423 throw DataException("Error - loadMesh:: unable to read Points_Nodes from netCDF file: " + *fName);
424 }
425 if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
426 TMPMEMFREE(Points_Nodes);
427 throw DataException("Error - load:: unable to recover Points_Nodes from netCDF file: " + *fName);
428 }
429 // Copy temp array into mesh_p->Points->Nodes
430 for (int i=0; i<num_Points; i++) {
431 mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
432 }
433 TMPMEMFREE(Points_Nodes);
434
435 } /* num_Points>0 */
436 }
437 }
438
439 /* get the tags */
440 if (Finley_noError()) {
441 if (num_Tags>0) {
442 // Temp storage to gather node IDs
443 int *Tags_keys = TMPMEMALLOC(num_Tags, int);
444 char name_temp[4096];
445 int i;
446
447 // Tags_keys
448 if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) )
449 throw DataException("Error - loadMesh:: unable to read Tags_keys from netCDF file: " + *fName);
450 if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
451 TMPMEMFREE(Tags_keys);
452 throw DataException("Error - loadMesh:: unable to recover Tags_keys from NetCDF file: " + *fName);
453 }
454 for (i=0; i<num_Tags; i++) {
455 // Retrieve tag name
456 sprintf(name_temp, "Tags_name_%d", i);
457 if (! (attr=dataFile.get_att(name_temp)) ) {
458 sprintf(error_msg,"Error retrieving tag name from NetCDF file '%s'", fName);
459 throw DataException(error_msg);
460 }
461 char *name = attr->as_string(0);
462 delete attr;
463 Finley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
464 }
465 }
466 }
467
468 if (Finley_noError()) Finley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
469 TMPMEMFREE(first_DofComponent);
470 TMPMEMFREE(first_NodeComponent);
471
472 } /* Finley_noError() after Finley_Mesh_alloc() */
473
474 checkFinleyError();
475 temp=new MeshAdapter(mesh_p);
476
477 if (! Finley_noError()) {
478 Finley_Mesh_free(mesh_p);
479 }
480
481 /* win32 refactor */
482 TMPMEMFREE(fName);
483
484 blocktimer_increment("LoadMesh()", blocktimer_start);
485 return temp->getPtr();
486 #else
487 throw DataException("Error - loadMesh: is not compiled with NetCFD. Please contact your installation manager.");
488 #endif /* USE_NETCDF */
489 }
490
491 Domain_ptr readMesh(const std::string& fileName,
492 int integrationOrder,
493 int reducedIntegrationOrder,
494 int optimize)
495 {
496 //
497 // create a copy of the filename to overcome the non-constness of call
498 // to Finley_Mesh_read
499 Finley_Mesh* fMesh=0;
500 // Win32 refactor
501 if( fileName.size() == 0 )
502 {
503 throw DataException("Null file name!");
504 }
505
506 char *fName = TMPMEMALLOC(fileName.size()+1,char);
507
508 strcpy(fName,fileName.c_str());
509 double blocktimer_start = blocktimer_time();
510
511 fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
512 checkFinleyError();
513 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
514
515 /* win32 refactor */
516 TMPMEMFREE(fName);
517
518 blocktimer_increment("ReadMesh()", blocktimer_start);
519 return temp->getPtr();
520 }
521
522 Domain_ptr readGmsh(const std::string& fileName,
523 int numDim,
524 int integrationOrder,
525 int reducedIntegrationOrder,
526 int optimize,
527 int useMacroElements)
528 {
529 //
530 // create a copy of the filename to overcome the non-constness of call
531 // to Finley_Mesh_read
532 Finley_Mesh* fMesh=0;
533 // Win32 refactor
534 if( fileName.size() == 0 )
535 {
536 throw DataException("Null file name!");
537 }
538
539 char *fName = TMPMEMALLOC(fileName.size()+1,char);
540
541 strcpy(fName,fileName.c_str());
542 double blocktimer_start = blocktimer_time();
543
544 fMesh=Finley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE), (useMacroElements ? TRUE : FALSE));
545 checkFinleyError();
546 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
547
548 /* win32 refactor */
549 TMPMEMFREE(fName);
550
551 blocktimer_increment("ReadGmsh()", blocktimer_start);
552 return temp->getPtr();
553 }
554
555 /* AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,*/
556 Domain_ptr brick(int n0,int n1,int n2,int order,
557 double l0,double l1,double l2,
558 int periodic0,int periodic1,
559 int periodic2,
560 int integrationOrder,
561 int reducedIntegrationOrder,
562 int useElementsOnFace,
563 int useFullElementOrder,
564 int optimize)
565 {
566 int numElements[]={n0,n1,n2};
567 double length[]={l0,l1,l2};
568 int periodic[]={periodic0, periodic1, periodic2};
569
570 //
571 // linearInterpolation
572 Finley_Mesh* fMesh=NULL;
573
574 if (order==1) {
575 fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
576 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
577 } else if (order==2) {
578 fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
579 useElementsOnFace,useFullElementOrder,FALSE, (optimize ? TRUE : FALSE)) ;
580 } else if (order==-1) {
581 fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
582 useElementsOnFace,useFullElementOrder,TRUE,(optimize ? TRUE : FALSE)) ;
583 } else {
584 stringstream temp;
585 temp << "Illegal interpolation order: " << order;
586 setFinleyError(VALUE_ERROR,temp.str().c_str());
587 }
588 //
589 // Convert any finley errors into a C++ exception
590 checkFinleyError();
591 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
592 return temp->getPtr();
593 }
594
595 /* AbstractContinuousDomain* rectangle(int n0,int n1,int order,*/
596 Domain_ptr rectangle(int n0,int n1,int order,
597 double l0, double l1,
598 int periodic0,int periodic1,
599 int integrationOrder,
600 int reducedIntegrationOrder,
601 int useElementsOnFace,
602 int useFullElementOrder,
603 int optimize)
604 {
605 int numElements[]={n0,n1};
606 double length[]={l0,l1};
607 int periodic[]={periodic0, periodic1};
608
609 Finley_Mesh* fMesh=0;
610 if (order==1) {
611 fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
612 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
613 } else if (order==2) {
614 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
615 useElementsOnFace,useFullElementOrder,FALSE,(optimize ? TRUE : FALSE));
616 } else if (order==-1) {
617 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
618 useElementsOnFace,useFullElementOrder,TRUE,(optimize ? TRUE : FALSE));
619 } else {
620 stringstream temp;
621 temp << "Illegal interpolation order: " << order;
622 setFinleyError(VALUE_ERROR,temp.str().c_str());
623 }
624 //
625 // Convert any finley errors into a C++ exception
626 checkFinleyError();
627 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
628 return temp->getPtr();
629 }
630
631 Domain_ptr meshMerge(const boost::python::list& meshList)
632 {
633 Finley_Mesh* fMesh=0;
634 //
635 // extract the meshes from meshList
636 int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
637 Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
638 for (int i=0;i<numMsh;++i) {
639 AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
640 const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
641 mshes[i]=finley_meshListMember->getFinley_Mesh();
642 }
643 //
644 // merge the meshes:
645 fMesh=Finley_Mesh_merge(numMsh,mshes);
646 TMPMEMFREE(mshes);
647 //
648 // Convert any finley errors into a C++ exception
649 checkFinleyError();
650 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
651
652 return temp->getPtr();
653 }
654
655 Domain_ptr glueFaces(const boost::python::list& meshList,
656 double safety_factor,
657 double tolerance,
658 int optimize)
659 {
660 Finley_Mesh* fMesh=0;
661 //
662 // merge the meshes:
663 Domain_ptr merged_meshes=meshMerge(meshList);
664
665 //
666 // glue the faces:
667 const MeshAdapter* merged_finley_meshes=dynamic_cast<const MeshAdapter*>(merged_meshes.get());
668 fMesh=merged_finley_meshes->getFinley_Mesh();
669 Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
670
671 //
672 // Convert any finley errors into a C++ exception
673 checkFinleyError();
674 return merged_meshes->getPtr();
675 }
676 Domain_ptr joinFaces(const boost::python::list& meshList,
677 double safety_factor,
678 double tolerance,
679 int optimize)
680 {
681 Finley_Mesh* fMesh=0;
682 //
683 // merge the meshes:
684 Domain_ptr merged_meshes=meshMerge(meshList);
685 //
686 // join the faces:
687 const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes.get());
688 fMesh=merged_finley_meshes->getFinley_Mesh();
689 Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
690 //
691 // Convert any finley errors into a C++ exception
692 checkFinleyError();
693 return merged_meshes->getPtr();
694 }
695
696 // end of namespace
697
698 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26