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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3554 - (show annotations)
Tue Aug 23 00:23:53 2011 UTC (7 years, 9 months ago) by jfenwick
File size: 32699 byte(s)
First step towards moving dirac

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 const boost::python::list& points,
607 const boost::python::list& tags
608 )
609 {
610 int numElements[]={n0,n1,n2};
611 double length[]={l0,l1,l2};
612 int periodic[]={periodic0, periodic1, periodic2};
613
614 //
615 // linearInterpolation
616 Finley_Mesh* fMesh=NULL;
617
618 if (order==1) {
619 fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
620 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
621 } else if (order==2) {
622 fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
623 useElementsOnFace,useFullElementOrder,FALSE, (optimize ? TRUE : FALSE)) ;
624 } else if (order==-1) {
625 fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
626 useElementsOnFace,useFullElementOrder,TRUE,(optimize ? TRUE : FALSE)) ;
627 } else {
628 stringstream temp;
629 temp << "Illegal interpolation order: " << order;
630 setFinleyError(VALUE_ERROR,temp.str().c_str());
631 }
632 //
633 // Convert any finley errors into a C++ exception
634 checkFinleyError();
635 MeshAdapter* temp=new MeshAdapter(fMesh);
636 temp->addDiracPoints(points, tags);
637 return temp->getPtr();
638 }
639
640
641
642 Domain_ptr brick_driver(const boost::python::list& args)
643 {
644 return brick(boost::python::extract<int>(args[0]), boost::python::extract<int>(args[1]), boost::python::extract<int>(args[2]),
645 boost::python::extract<int>(args[3]), boost::python::extract<double>(args[4]),
646 boost::python::extract<double>(args[5]),
647 boost::python::extract<double>(args[6]), boost::python::extract<int>(args[7]),
648 boost::python::extract<int>(args[8]),
649 boost::python::extract<int>(args[9]),
650 boost::python::extract<int>(args[10]), boost::python::extract<int>(args[11]),
651 boost::python::extract<int>(args[12]),
652 boost::python::extract<int>(args[13]), boost::python::extract<int>(args[14]),
653 boost::python::extract<boost::python::list>(args[15]),
654 boost::python::extract<boost::python::list>(args[16]));
655 }
656
657 /* AbstractContinuousDomain* rectangle(int n0,int n1,int order,*/
658 Domain_ptr rectangle(int n0,int n1,int order,
659 double l0, double l1,
660 int periodic0,int periodic1,
661 int integrationOrder,
662 int reducedIntegrationOrder,
663 int useElementsOnFace,
664 int useFullElementOrder,
665 int optimize,
666 const boost::python::list& points,
667 const boost::python::list& tags)
668 {
669 int numElements[]={n0,n1};
670 double length[]={l0,l1};
671 int periodic[]={periodic0, periodic1};
672
673 Finley_Mesh* fMesh=0;
674 if (order==1) {
675 fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
676 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
677 } else if (order==2) {
678 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
679 useElementsOnFace,useFullElementOrder,FALSE,(optimize ? TRUE : FALSE));
680 } else if (order==-1) {
681 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
682 useElementsOnFace,useFullElementOrder,TRUE,(optimize ? TRUE : FALSE));
683 } else {
684 stringstream temp;
685 temp << "Illegal interpolation order: " << order;
686 setFinleyError(VALUE_ERROR,temp.str().c_str());
687 }
688 //
689 // Convert any finley errors into a C++ exception
690 checkFinleyError();
691 MeshAdapter* temp=new MeshAdapter(fMesh);
692 temp->addDiracPoints(points, tags);
693 return temp->getPtr();
694 }
695
696 Domain_ptr meshMerge(const boost::python::list& meshList)
697 {
698 Finley_Mesh* fMesh=0;
699 //
700 // extract the meshes from meshList
701 int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
702 Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
703 for (int i=0;i<numMsh;++i) {
704 AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
705 const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
706 mshes[i]=finley_meshListMember->getFinley_Mesh();
707 }
708 //
709 // merge the meshes:
710 fMesh=Finley_Mesh_merge(numMsh,mshes);
711 TMPMEMFREE(mshes);
712 //
713 // Convert any finley errors into a C++ exception
714 checkFinleyError();
715 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
716
717 return temp->getPtr();
718 }
719
720 Domain_ptr glueFaces(const boost::python::list& meshList,
721 double safety_factor,
722 double tolerance,
723 int optimize)
724 {
725 Finley_Mesh* fMesh=0;
726 //
727 // merge the meshes:
728 Domain_ptr merged_meshes=meshMerge(meshList);
729
730 //
731 // glue the faces:
732 const MeshAdapter* merged_finley_meshes=dynamic_cast<const MeshAdapter*>(merged_meshes.get());
733 fMesh=merged_finley_meshes->getFinley_Mesh();
734 Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
735
736 //
737 // Convert any finley errors into a C++ exception
738 checkFinleyError();
739 return merged_meshes->getPtr();
740 }
741 Domain_ptr joinFaces(const boost::python::list& meshList,
742 double safety_factor,
743 double tolerance,
744 int optimize)
745 {
746 Finley_Mesh* fMesh=0;
747 //
748 // merge the meshes:
749 Domain_ptr merged_meshes=meshMerge(meshList);
750 //
751 // join the faces:
752 const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes.get());
753 fMesh=merged_finley_meshes->getFinley_Mesh();
754 Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
755 //
756 // Convert any finley errors into a C++ exception
757 checkFinleyError();
758 return merged_meshes->getPtr();
759 }
760
761 // end of namespace
762
763 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26