/[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 3317 - (show annotations)
Thu Oct 28 00:50:41 2010 UTC (8 years, 9 months ago) by caltinay
File size: 31640 byte(s)
Removed bogus mpi.h includes, replaced others by our esysUtils wrapper
and rearranged so that the wrapper always comes before netcdf which fixes
linking problems when disabling mpi on systems with netCDF 4.x.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26