/[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 3555 - (show annotations)
Thu Aug 25 08:54:10 2011 UTC (7 years, 9 months ago) by jfenwick
File size: 36286 byte(s)
Committing for home transfer

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 std::vector<double>& points,
607 const std::vector<int>& 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 using boost::python::extract;
645
646 // we need to convert lists to stl vectors
647 boost::python::list pypoints=boost::python::extract<boost::python::list>(args[15]);
648 boost::python::list pytags=boost::python::extract<boost::python::list>(args[16]);
649 int numpts=boost::python::extract<int>(pypoints.attr("__len__")());
650 int numtags=boost::python::extract<int>(pytags.attr("__len__")());
651 vector<double> points;
652 vector<int> tags;
653 tags.resize(numtags, -1);
654 for (int i=0;i<numpts;++i)
655 {
656 boost::python::object temp=pypoints[i];
657 int l=boost::python::extract<int>(temp.attr("__len__")());
658 for (int k=0;k<l;++k)
659 {
660 points.push_back(extract<double>(temp[k]));
661 }
662 }
663 map<string, int> tagstonames;
664 int curmax=-1;
665 // but which order to assign tags to names?????
666 for (int i=0;i<numtags;++i)
667 {
668 extract<int> ex_int(pytags[i]);
669 extract<string> ex_str(pytags[i]);
670 if (ex_int.check())
671 {
672 tags[i]=ex_int();
673 if (tags[i]> curmax)
674 {
675 curmax=tags[i];
676 }
677 }
678 else if (ex_str.check())
679 {
680 string s=ex_str();
681 map<string, int>::iterator it=tagstonames.find(s);
682 if (it!=tagstonames.end())
683 {
684 // we have the tag already so look it up
685 tags[i]=it->second;
686 }
687 else
688 {
689 tagstonames[s]=curmax;
690 tags[i]=curmax;
691 curmax++;
692 }
693 }
694 else
695 {
696 throw FinleyAdapterException("Error - Unable to extract tag value.");
697 }
698
699 }
700
701 return brick(boost::python::extract<int>(args[0]), boost::python::extract<int>(args[1]), boost::python::extract<int>(args[2]),
702 boost::python::extract<int>(args[3]), boost::python::extract<double>(args[4]),
703 boost::python::extract<double>(args[5]),
704 boost::python::extract<double>(args[6]), boost::python::extract<int>(args[7]),
705 boost::python::extract<int>(args[8]),
706 boost::python::extract<int>(args[9]),
707 boost::python::extract<int>(args[10]), boost::python::extract<int>(args[11]),
708 boost::python::extract<int>(args[12]),
709 boost::python::extract<int>(args[13]), boost::python::extract<int>(args[14]),
710 points,
711 tags);
712 }
713
714 /* AbstractContinuousDomain* rectangle(int n0,int n1,int order,*/
715 Domain_ptr rectangle(int n0,int n1,int order,
716 double l0, double l1,
717 int periodic0,int periodic1,
718 int integrationOrder,
719 int reducedIntegrationOrder,
720 int useElementsOnFace,
721 int useFullElementOrder,
722 int optimize,
723 const vector<double>& points,
724 const vector<int>& tags)
725 {
726 int numElements[]={n0,n1};
727 double length[]={l0,l1};
728 int periodic[]={periodic0, periodic1};
729
730 Finley_Mesh* fMesh=0;
731 if (order==1) {
732 fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
733 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
734 } else if (order==2) {
735 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
736 useElementsOnFace,useFullElementOrder,FALSE,(optimize ? TRUE : FALSE));
737 } else if (order==-1) {
738 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
739 useElementsOnFace,useFullElementOrder,TRUE,(optimize ? TRUE : FALSE));
740 } else {
741 stringstream temp;
742 temp << "Illegal interpolation order: " << order;
743 setFinleyError(VALUE_ERROR,temp.str().c_str());
744 }
745 //
746 // Convert any finley errors into a C++ exception
747 checkFinleyError();
748 MeshAdapter* temp=new MeshAdapter(fMesh);
749 temp->addDiracPoints(points, tags);
750 return temp->getPtr();
751 }
752
753 Domain_ptr meshMerge(const boost::python::list& meshList)
754 {
755 Finley_Mesh* fMesh=0;
756 //
757 // extract the meshes from meshList
758 int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
759 Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
760 for (int i=0;i<numMsh;++i) {
761 AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
762 const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
763 mshes[i]=finley_meshListMember->getFinley_Mesh();
764 }
765 //
766 // merge the meshes:
767 fMesh=Finley_Mesh_merge(numMsh,mshes);
768 TMPMEMFREE(mshes);
769 //
770 // Convert any finley errors into a C++ exception
771 checkFinleyError();
772 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
773
774 return temp->getPtr();
775 }
776
777 Domain_ptr rectangle_driver(const boost::python::list& args)
778 {
779 using boost::python::extract;
780
781 // we need to convert lists to stl vectors
782 boost::python::list pypoints=boost::python::extract<boost::python::list>(args[12]);
783 boost::python::list pytags=boost::python::extract<boost::python::list>(args[13]);
784 int numpts=boost::python::extract<int>(pypoints.attr("__len__")());
785 int numtags=boost::python::extract<int>(pytags.attr("__len__")());
786 vector<double> points;
787 vector<int> tags;
788 tags.resize(numtags, -1);
789 for (int i=0;i<numpts;++i)
790 {
791 boost::python::object temp=pypoints[i];
792 int l=boost::python::extract<int>(temp.attr("__len__")());
793 for (int k=0;k<l;++k)
794 {
795 points.push_back(extract<double>(temp[k]));
796 }
797 }
798 map<string, int> tagstonames;
799 int curmax=-1;
800 // but which order to assign tags to names?????
801 for (int i=0;i<numtags;++i)
802 {
803 extract<int> ex_int(pytags[i]);
804 extract<string> ex_str(pytags[i]);
805 if (ex_int.check())
806 {
807 tags[i]=ex_int();
808 if (tags[i]> curmax)
809 {
810 curmax=tags[i];
811 }
812 }
813 else if (ex_str.check())
814 {
815 string s=ex_str();
816 map<string, int>::iterator it=tagstonames.find(s);
817 if (it!=tagstonames.end())
818 {
819 // we have the tag already so look it up
820 tags[i]=it->second;
821 }
822 else
823 {
824 tagstonames[s]=curmax;
825 tags[i]=curmax;
826 curmax++;
827 }
828 }
829 else
830 {
831 throw FinleyAdapterException("Error - Unable to extract tag value.");
832 }
833
834 }
835
836 return rectangle(boost::python::extract<int>(args[0]), boost::python::extract<int>(args[1]), boost::python::extract<int>(args[2]),
837 boost::python::extract<double>(args[3]),
838 boost::python::extract<double>(args[4]),
839 boost::python::extract<int>(args[5]), boost::python::extract<int>(args[6]),
840 boost::python::extract<int>(args[7]),
841 boost::python::extract<int>(args[8]),
842 boost::python::extract<int>(args[9]), boost::python::extract<int>(args[10]),
843 boost::python::extract<int>(args[11]),
844 points,
845 tags);
846 }
847
848
849 Domain_ptr glueFaces(const boost::python::list& meshList,
850 double safety_factor,
851 double tolerance,
852 int optimize)
853 {
854 Finley_Mesh* fMesh=0;
855 //
856 // merge the meshes:
857 Domain_ptr merged_meshes=meshMerge(meshList);
858
859 //
860 // glue the faces:
861 const MeshAdapter* merged_finley_meshes=dynamic_cast<const MeshAdapter*>(merged_meshes.get());
862 fMesh=merged_finley_meshes->getFinley_Mesh();
863 Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
864
865 //
866 // Convert any finley errors into a C++ exception
867 checkFinleyError();
868 return merged_meshes->getPtr();
869 }
870 Domain_ptr joinFaces(const boost::python::list& meshList,
871 double safety_factor,
872 double tolerance,
873 int optimize)
874 {
875 Finley_Mesh* fMesh=0;
876 //
877 // merge the meshes:
878 Domain_ptr merged_meshes=meshMerge(meshList);
879 //
880 // join the faces:
881 const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes.get());
882 fMesh=merged_finley_meshes->getFinley_Mesh();
883 Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
884 //
885 // Convert any finley errors into a C++ exception
886 checkFinleyError();
887 return merged_meshes->getPtr();
888 }
889
890 // end of namespace
891
892 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26