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

Contents of /branches/domexper/finley/src/CPPAdapter/MeshAdapterFactory.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3247 - (show annotations)
Wed Oct 6 05:53:06 2010 UTC (8 years, 6 months ago) by caltinay
File size: 31626 byte(s)
Fixed name clashes between dudley and finley so both can be used
simultaneously.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26