/[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 3637 - (show annotations)
Fri Oct 21 04:29:33 2011 UTC (7 years, 9 months ago) by caltinay
File size: 38959 byte(s)
Fixed memory leaks in loadMesh() for dudley and finley.

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,"loadMesh: 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 inline void cleanupAndThrow(Finley_Mesh* mesh, Esys_MPIInfo* info, string msg)
51 {
52 Finley_Mesh_free(mesh);
53 Esys_MPIInfo_free(info);
54 string msgPrefix("loadMesh: NetCDF operation failed - ");
55 throw DataException(msgPrefix+msg);
56 }
57
58 // AbstractContinuousDomain* loadMesh(const std::string& fileName)
59 Domain_ptr loadMesh(const std::string& fileName)
60 {
61 #ifdef USE_NETCDF
62 Esys_MPIInfo *mpi_info = Esys_MPIInfo_alloc( MPI_COMM_WORLD );
63 Finley_Mesh *mesh_p=NULL;
64 char error_msg[LenErrorMsg_MAX];
65
66 char *fName = Esys_MPI_appendRankToFileName(fileName.c_str(),
67 mpi_info->size,
68 mpi_info->rank);
69
70 double blocktimer_start = blocktimer_time();
71 Finley_resetError();
72 int *first_DofComponent, *first_NodeComponent;
73
74 // Open NetCDF file for reading
75 NcAtt *attr;
76 NcVar *nc_var_temp;
77 // netCDF error handler
78 NcError err(NcError::silent_nonfatal);
79 // Create the NetCDF file.
80 NcFile dataFile(fName, NcFile::ReadOnly);
81 if (!dataFile.is_valid()) {
82 sprintf(error_msg,"loadMesh: Opening NetCDF file '%s' for reading failed.", fName);
83 Finley_setError(IO_ERROR,error_msg);
84 Esys_MPIInfo_free( mpi_info );
85 throw DataException(error_msg);
86 }
87
88 // Read NetCDF integer attributes
89 int mpi_size = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_size");
90 int mpi_rank = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_rank");
91 int numDim = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numDim");
92 int order = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"order");
93 int reduced_order = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"reduced_order");
94 int numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numNodes");
95 int num_Elements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements");
96 int num_FaceElements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements");
97 int num_ContactElements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_ContactElements");
98 int num_Points = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Points");
99 int num_Elements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements_numNodes");
100 int Elements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Elements_TypeId");
101 int num_FaceElements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements_numNodes");
102 int FaceElements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"FaceElements_TypeId");
103 int num_ContactElements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_ContactElements_numNodes");
104 int ContactElements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"ContactElements_TypeId");
105 int Points_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Points_TypeId");
106 int num_Tags = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Tags");
107
108 // Verify size and rank
109 if (mpi_info->size != mpi_size) {
110 sprintf(error_msg, "loadMesh: The NetCDF file '%s' can only be read on %d CPUs instead of %d", fName, mpi_size, mpi_info->size);
111 throw DataException(error_msg);
112 }
113 if (mpi_info->rank != mpi_rank) {
114 sprintf(error_msg, "loadMesh: The NetCDF file '%s' should be read on CPU #%d instead of %d", fName, mpi_rank, mpi_info->rank);
115 throw DataException(error_msg);
116 }
117
118 // Read mesh name
119 if (! (attr=dataFile.get_att("Name")) ) {
120 sprintf(error_msg,"loadMesh: Error retrieving mesh name from NetCDF file '%s'", fName);
121 throw DataException(error_msg);
122 }
123 char *name = attr->as_string(0);
124 delete attr;
125
126 TMPMEMFREE(fName);
127
128 /* allocate mesh */
129 mesh_p = Finley_Mesh_alloc(name,numDim,mpi_info);
130 if (Finley_noError()) {
131
132 /* read nodes */
133 Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
134 // Nodes_Id
135 if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
136 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_Id)");
137 if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) )
138 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_Id)");
139 // Nodes_Tag
140 if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
141 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_Tag)");
142 if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) )
143 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_Tag)");
144 // Nodes_gDOF
145 if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
146 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_gDOF)");
147 if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) )
148 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_gDOF)");
149 // Nodes_gNI
150 if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
151 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_gNI)");
152 if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) )
153 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_gNI)");
154 // Nodes_grDfI
155 if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
156 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_grDfI)");
157 if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) )
158 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_grDfI)");
159 // Nodes_grNI
160 if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
161 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_grNI)");
162 if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) )
163 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_grNI)");
164 // Nodes_Coordinates
165 if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates")))
166 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_Coordinates)");
167 if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) )
168 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_Coordinates)");
169
170 /* read elements */
171 if (Finley_noError()) {
172 Finley_ReferenceElementSet *refElements= Finley_ReferenceElementSet_alloc((Finley_ElementTypeId)Elements_TypeId,order, reduced_order);
173 if (Finley_noError()) {
174 mesh_p->Elements=Finley_ElementFile_alloc(refElements, mpi_info);
175 }
176 Finley_ReferenceElementSet_dealloc(refElements);
177 if (Finley_noError())
178 Finley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
179 if (Finley_noError()) {
180 mesh_p->Elements->minColor=0;
181 mesh_p->Elements->maxColor=num_Elements-1;
182 if (num_Elements>0) {
183 // Elements_Id
184 if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
185 cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Id)");
186 if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) )
187 cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Id)");
188 // Elements_Tag
189 if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
190 cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Tag)");
191 if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) )
192 cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Tag)");
193 // Elements_Owner
194 if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
195 cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Owner)");
196 if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) )
197 cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Owner)");
198 // Elements_Color
199 if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
200 cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Color)");
201 if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) )
202 cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Color)");
203 // Now we need to adjust maxColor
204 index_t mc=mesh_p->Elements->Color[0];
205 for (index_t i=1;i<num_Elements;++i) {
206 if (mc<mesh_p->Elements->Color[i]) {
207 mc = mesh_p->Elements->Color[i];
208 }
209 }
210 mesh_p->Elements->maxColor=mc;
211 // Elements_Nodes
212 int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
213 if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
214 TMPMEMFREE(Elements_Nodes);
215 cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Nodes)");
216 }
217 if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
218 TMPMEMFREE(Elements_Nodes);
219 cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Nodes)");
220 }
221
222 // Copy temp array into mesh_p->Elements->Nodes
223 for (int i=0; i<num_Elements; i++) {
224 for (int j=0; j<num_Elements_numNodes; j++) {
225 mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
226 = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
227 }
228 }
229 TMPMEMFREE(Elements_Nodes);
230 } /* num_Elements>0 */
231 }
232 }
233
234 /* get the face elements */
235 if (Finley_noError()) {
236 Finley_ReferenceElementSet *refFaceElements =
237 Finley_ReferenceElementSet_alloc((Finley_ElementTypeId)FaceElements_TypeId, order, reduced_order);
238 if (Finley_noError()) {
239 mesh_p->FaceElements=Finley_ElementFile_alloc(refFaceElements, mpi_info);
240 }
241 Finley_ReferenceElementSet_dealloc(refFaceElements);
242 if (Finley_noError())
243 Finley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
244 if (Finley_noError()) {
245 mesh_p->FaceElements->minColor=0;
246 mesh_p->FaceElements->maxColor=num_FaceElements-1;
247 if (num_FaceElements>0) {
248 // FaceElements_Id
249 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
250 cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Id)");
251 if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) )
252 cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Id)");
253 // FaceElements_Tag
254 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
255 cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Tag)");
256 if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) )
257 cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Tag)");
258 // FaceElements_Owner
259 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
260 cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Owner)");
261 if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) )
262 cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Owner)");
263 // FaceElements_Color
264 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
265 cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Color)");
266 if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) )
267 cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Color)");
268 // Now we need to adjust maxColor
269 index_t mc=mesh_p->FaceElements->Color[0];
270 for (index_t i=1;i<num_FaceElements;++i) {
271 if (mc<mesh_p->FaceElements->Color[i]) {
272 mc = mesh_p->FaceElements->Color[i];
273 }
274 }
275 mesh_p->FaceElements->maxColor=mc;
276 // FaceElements_Nodes
277 int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
278 if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
279 TMPMEMFREE(FaceElements_Nodes);
280 cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Nodes)");
281 }
282 if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
283 TMPMEMFREE(FaceElements_Nodes);
284 cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Nodes)");
285 }
286 // Copy temp array into mesh_p->FaceElements->Nodes
287 for (int i=0; i<num_FaceElements; i++) {
288 for (int j=0; j<num_FaceElements_numNodes; j++) {
289 mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
290 }
291 }
292 TMPMEMFREE(FaceElements_Nodes);
293 } /* num_FaceElements>0 */
294 }
295 }
296
297 /* get the Contact elements */
298 if (Finley_noError()) {
299 Finley_ReferenceElementSet *refContactElements =
300 Finley_ReferenceElementSet_alloc((Finley_ElementTypeId)ContactElements_TypeId, order, reduced_order);
301 if (Finley_noError()) {
302 mesh_p->ContactElements=Finley_ElementFile_alloc(refContactElements, mpi_info);
303 }
304 Finley_ReferenceElementSet_dealloc(refContactElements);
305 if (Finley_noError())
306 Finley_ElementFile_allocTable(mesh_p->ContactElements, num_ContactElements);
307 if (Finley_noError()) {
308 mesh_p->ContactElements->minColor=0;
309 mesh_p->ContactElements->maxColor=num_ContactElements-1;
310 if (num_ContactElements>0) {
311 // ContactElements_Id
312 if (! ( nc_var_temp = dataFile.get_var("ContactElements_Id")) )
313 cleanupAndThrow(mesh_p, mpi_info, "get_var(ContactElements_Id)");
314 if (! nc_var_temp->get(&mesh_p->ContactElements->Id[0], num_ContactElements) )
315 cleanupAndThrow(mesh_p, mpi_info, "get(ContactElements_Id)");
316 // ContactElements_Tag
317 if (! ( nc_var_temp = dataFile.get_var("ContactElements_Tag")) )
318 cleanupAndThrow(mesh_p, mpi_info, "get_var(ContactElements_Tag)");
319 if (! nc_var_temp->get(&mesh_p->ContactElements->Tag[0], num_ContactElements) )
320 cleanupAndThrow(mesh_p, mpi_info, "get(ContactElements_Tag)");
321 // ContactElements_Owner
322 if (! ( nc_var_temp = dataFile.get_var("ContactElements_Owner")) )
323 cleanupAndThrow(mesh_p, mpi_info, "get_var(ContactElements_Owner)");
324 if (! nc_var_temp->get(&mesh_p->ContactElements->Owner[0], num_ContactElements) )
325 cleanupAndThrow(mesh_p, mpi_info, "get(ContactElements_Owner)");
326 // ContactElements_Color
327 if (! ( nc_var_temp = dataFile.get_var("ContactElements_Color")) )
328 cleanupAndThrow(mesh_p, mpi_info, "get_var(ContactElements_Color)");
329 if (! nc_var_temp->get(&mesh_p->ContactElements->Color[0], num_ContactElements) )
330 cleanupAndThrow(mesh_p, mpi_info, "get(ContactElements_Color)");
331 // Now we need to adjust maxColor
332 index_t mc=mesh_p->ContactElements->Color[0];
333 for (index_t i=1;i<num_ContactElements;++i) {
334 if (mc<mesh_p->ContactElements->Color[i]) {
335 mc = mesh_p->ContactElements->Color[i];
336 }
337 }
338 mesh_p->ContactElements->maxColor=mc;
339 // ContactElements_Nodes
340 int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
341 if (!(nc_var_temp = dataFile.get_var("ContactElements_Nodes"))) {
342 TMPMEMFREE(ContactElements_Nodes);
343 cleanupAndThrow(mesh_p, mpi_info, "get_var(ContactElements_Nodes)");
344 }
345 if (! nc_var_temp->get(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes) ) {
346 TMPMEMFREE(ContactElements_Nodes);
347 cleanupAndThrow(mesh_p, mpi_info, "get(ContactElements_Nodes)");
348 }
349 // Copy temp array into mesh_p->ContactElements->Nodes
350 for (int i=0; i<num_ContactElements; i++) {
351 for (int j=0; j<num_ContactElements_numNodes; j++) {
352 mesh_p->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]= ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
353 }
354 }
355 TMPMEMFREE(ContactElements_Nodes);
356 } /* num_ContactElements>0 */
357 }
358 }
359
360 /* get the Points (nodal elements) */
361 if (Finley_noError()) {
362 Finley_ReferenceElementSet *refPoints =
363 Finley_ReferenceElementSet_alloc((Finley_ElementTypeId)Points_TypeId,order, reduced_order);
364 if (Finley_noError()) {
365 mesh_p->Points=Finley_ElementFile_alloc(refPoints, mpi_info);
366 }
367 Finley_ReferenceElementSet_dealloc(refPoints);
368 if (Finley_noError())
369 Finley_ElementFile_allocTable(mesh_p->Points, num_Points);
370 if (Finley_noError()) {
371 mesh_p->Points->minColor=0;
372 mesh_p->Points->maxColor=num_Points-1;
373 if (num_Points>0) {
374 // Points_Id
375 if (! ( nc_var_temp = dataFile.get_var("Points_Id")))
376 cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Id)");
377 if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points))
378 cleanupAndThrow(mesh_p, mpi_info, "get(Points_Id)");
379 // Points_Tag
380 if (! ( nc_var_temp = dataFile.get_var("Points_Tag")))
381 cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Tag)");
382 if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points))
383 cleanupAndThrow(mesh_p, mpi_info, "get(Points_Tag)");
384 // Points_Owner
385 if (! ( nc_var_temp = dataFile.get_var("Points_Owner")))
386 cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Owner)");
387 if (!nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points))
388 cleanupAndThrow(mesh_p, mpi_info, "get(Points_Owner)");
389 // Points_Color
390 if (! ( nc_var_temp = dataFile.get_var("Points_Color")))
391 cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Color)");
392 if (!nc_var_temp->get(&mesh_p->Points->Color[0], num_Points))
393 cleanupAndThrow(mesh_p, mpi_info, "get(Points_Color)");
394 // Now we need to adjust maxColor
395 index_t mc=mesh_p->Points->Color[0];
396 for (index_t i=1;i<num_Points;++i) {
397 if (mc<mesh_p->Points->Color[i]) {
398 mc = mesh_p->Points->Color[i];
399 }
400 }
401 mesh_p->Points->maxColor=mc;
402 // Points_Nodes
403 int *Points_Nodes = TMPMEMALLOC(num_Points,int);
404 if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
405 TMPMEMFREE(Points_Nodes);
406 cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Nodes)");
407 }
408 if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
409 TMPMEMFREE(Points_Nodes);
410 cleanupAndThrow(mesh_p, mpi_info, "get(Points_Nodes)");
411 }
412 // Copy temp array into mesh_p->Points->Nodes
413 for (int i=0; i<num_Points; i++) {
414 mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
415 }
416 TMPMEMFREE(Points_Nodes);
417 } /* num_Points>0 */
418 }
419 }
420
421 /* get the tags */
422 if (Finley_noError()) {
423 if (num_Tags>0) {
424 // Temp storage to gather node IDs
425 int *Tags_keys = TMPMEMALLOC(num_Tags, int);
426 char name_temp[4096];
427 int i;
428
429 // Tags_keys
430 if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) ) {
431 TMPMEMFREE(Tags_keys);
432 cleanupAndThrow(mesh_p, mpi_info, "get_var(Tags_keys)");
433 }
434 if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
435 TMPMEMFREE(Tags_keys);
436 cleanupAndThrow(mesh_p, mpi_info, "get(Tags_keys)");
437 }
438 for (i=0; i<num_Tags; i++) {
439 // Retrieve tag name
440 sprintf(name_temp, "Tags_name_%d", i);
441 if (! (attr=dataFile.get_att(name_temp)) ) {
442 TMPMEMFREE(Tags_keys);
443 sprintf(error_msg,"get_att(%s)", name_temp);
444 cleanupAndThrow(mesh_p, mpi_info, error_msg);
445 }
446 char *name = attr->as_string(0);
447 delete attr;
448 Finley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
449 }
450 TMPMEMFREE(Tags_keys);
451 }
452 }
453
454 if (Finley_noError()) {
455 // Nodes_DofDistribution
456 first_DofComponent = TMPMEMALLOC(mpi_size+1,index_t);
457 if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) ) {
458 TMPMEMFREE(first_DofComponent);
459 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_DofDistribution)");
460 }
461 if (! nc_var_temp->get(&first_DofComponent[0], mpi_size+1) ) {
462 TMPMEMFREE(first_DofComponent);
463 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_DofDistribution)");
464 }
465
466 // Nodes_NodeDistribution
467 first_NodeComponent = TMPMEMALLOC(mpi_size+1,index_t);
468 if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) ) {
469 TMPMEMFREE(first_DofComponent);
470 TMPMEMFREE(first_NodeComponent);
471 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_NodeDistribution)");
472 }
473 if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
474 TMPMEMFREE(first_DofComponent);
475 TMPMEMFREE(first_NodeComponent);
476 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_NodeDistribution)");
477 }
478 Finley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
479 TMPMEMFREE(first_DofComponent);
480 TMPMEMFREE(first_NodeComponent);
481 }
482
483 } /* Finley_noError() after Finley_Mesh_alloc() */
484
485 checkFinleyError();
486 AbstractContinuousDomain* dom=new MeshAdapter(mesh_p);
487
488 if (! Finley_noError()) {
489 Finley_Mesh_free(mesh_p);
490 }
491
492 blocktimer_increment("LoadMesh()", blocktimer_start);
493 return dom->getPtr();
494 #else
495 throw DataException("loadMesh: not compiled with NetCDF. Please contact your installation manager.");
496 #endif /* USE_NETCDF */
497 }
498
499 Domain_ptr readMesh(const std::string& fileName,
500 int integrationOrder,
501 int reducedIntegrationOrder,
502 int optimize)
503 {
504 //
505 // create a copy of the filename to overcome the non-constness of call
506 // to Finley_Mesh_read
507 Finley_Mesh* fMesh=0;
508 // Win32 refactor
509 if( fileName.size() == 0 )
510 {
511 throw DataException("Null file name!");
512 }
513
514 char *fName = TMPMEMALLOC(fileName.size()+1,char);
515
516 strcpy(fName,fileName.c_str());
517 double blocktimer_start = blocktimer_time();
518
519 fMesh=Finley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
520 checkFinleyError();
521 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
522
523 /* win32 refactor */
524 TMPMEMFREE(fName);
525
526 blocktimer_increment("ReadMesh()", blocktimer_start);
527 return temp->getPtr();
528 }
529
530 Domain_ptr readGmsh(const std::string& fileName,
531 int numDim,
532 int integrationOrder,
533 int reducedIntegrationOrder,
534 int optimize,
535 int useMacroElements)
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_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE), (useMacroElements ? TRUE : FALSE));
553 checkFinleyError();
554 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
555
556 /* win32 refactor */
557 TMPMEMFREE(fName);
558
559 blocktimer_increment("ReadGmsh()", blocktimer_start);
560 return temp->getPtr();
561 }
562
563 /* AbstractContinuousDomain* brick(int n0,int n1,int n2,int order,*/
564 Domain_ptr brick(int n0,int n1,int n2,int order,
565 double l0,double l1,double l2,
566 int periodic0,int periodic1,
567 int periodic2,
568 int integrationOrder,
569 int reducedIntegrationOrder,
570 int useElementsOnFace,
571 int useFullElementOrder,
572 int optimize,
573 const std::vector<double>& points,
574 const std::vector<int>& tags, const std::map<std::string, int>& tagnamestonums)
575 {
576 int numElements[]={n0,n1,n2};
577 double length[]={l0,l1,l2};
578 int periodic[]={periodic0, periodic1, periodic2};
579
580 //
581 // linearInterpolation
582 Finley_Mesh* fMesh=NULL;
583
584 if (order==1) {
585 fMesh=Finley_RectangularMesh_Hex8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
586 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE)) ;
587 } else if (order==2) {
588 fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
589 useElementsOnFace,useFullElementOrder,FALSE, (optimize ? TRUE : FALSE)) ;
590 } else if (order==-1) {
591 fMesh=Finley_RectangularMesh_Hex20(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
592 useElementsOnFace,useFullElementOrder,TRUE,(optimize ? TRUE : FALSE)) ;
593 } else {
594 stringstream temp;
595 temp << "Illegal interpolation order: " << order;
596 setFinleyError(VALUE_ERROR,temp.str().c_str());
597 }
598 //
599 // Convert any finley errors into a C++ exception
600 checkFinleyError();
601 MeshAdapter* temp=new MeshAdapter(fMesh);
602 temp->addDiracPoints(points, tags);
603 Finley_Mesh* out=temp->getMesh().get();
604 for (map<string, int>::const_iterator it=tagnamestonums.begin();it!=tagnamestonums.end();++it)
605 {
606 Finley_Mesh_addTagMap(out, it->first.c_str(), it->second);
607 }
608 return temp->getPtr();
609 }
610
611 Domain_ptr brick_driver(const boost::python::list& args)
612 {
613 using boost::python::extract;
614
615 // we need to convert lists to stl vectors
616 boost::python::list pypoints=extract<boost::python::list>(args[15]);
617 boost::python::list pytags=extract<boost::python::list>(args[16]);
618 int numpts=extract<int>(pypoints.attr("__len__")());
619 int numtags=extract<int>(pytags.attr("__len__")());
620 vector<double> points;
621 vector<int> tags;
622 tags.resize(numtags, -1);
623 for (int i=0;i<numpts;++i) {
624 boost::python::object temp=pypoints[i];
625 int l=extract<int>(temp.attr("__len__")());
626 for (int k=0;k<l;++k) {
627 points.push_back(extract<double>(temp[k]));
628 }
629 }
630 map<string, int> namestonums;
631 int curmax=40; // bricks use up to 30
632 for (int i=0;i<numtags;++i) {
633 extract<int> ex_int(pytags[i]);
634 extract<string> ex_str(pytags[i]);
635 if (ex_int.check()) {
636 tags[i]=ex_int();
637 if (tags[i]>= curmax) {
638 curmax=tags[i]+1;
639 }
640 } else if (ex_str.check()) {
641 string s=ex_str();
642 map<string, int>::iterator it=namestonums.find(s);
643 if (it!=namestonums.end()) {
644 // we have the tag already so look it up
645 tags[i]=it->second;
646 } else {
647 namestonums[s]=curmax;
648 tags[i]=curmax;
649 curmax++;
650 }
651 } else {
652 throw FinleyAdapterException("Error - Unable to extract tag value.");
653 }
654
655 }
656
657 Domain_ptr res=brick(extract<int>(args[0]), extract<int>(args[1]),
658 extract<int>(args[2]), extract<int>(args[3]),
659 extract<double>(args[4]), extract<double>(args[5]),
660 extract<double>(args[6]), extract<int>(args[7]),
661 extract<int>(args[8]), extract<int>(args[9]),
662 extract<int>(args[10]), extract<int>(args[11]),
663 extract<int>(args[12]), extract<int>(args[13]),
664 extract<int>(args[14]), points, tags, namestonums);
665
666 return res;
667 }
668
669 /* AbstractContinuousDomain* rectangle(int n0,int n1,int order,*/
670 Domain_ptr rectangle(int n0,int n1,int order,
671 double l0, double l1,
672 int periodic0,int periodic1,
673 int integrationOrder,
674 int reducedIntegrationOrder,
675 int useElementsOnFace,
676 int useFullElementOrder,
677 int optimize,
678 const vector<double>& points,
679 const vector<int>& tags,
680 const std::map<std::string, int>& tagnamestonums)
681 {
682 int numElements[]={n0,n1};
683 double length[]={l0,l1};
684 int periodic[]={periodic0, periodic1};
685
686 Finley_Mesh* fMesh=0;
687 if (order==1) {
688 fMesh=Finley_RectangularMesh_Rec4(numElements, length,periodic,integrationOrder,reducedIntegrationOrder,
689 useElementsOnFace,useFullElementOrder,(optimize ? TRUE : FALSE));
690 } else if (order==2) {
691 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
692 useElementsOnFace,useFullElementOrder,FALSE,(optimize ? TRUE : FALSE));
693 } else if (order==-1) {
694 fMesh=Finley_RectangularMesh_Rec8(numElements,length,periodic,integrationOrder,reducedIntegrationOrder,
695 useElementsOnFace,useFullElementOrder,TRUE,(optimize ? TRUE : FALSE));
696 } else {
697 stringstream temp;
698 temp << "Illegal interpolation order: " << order;
699 setFinleyError(VALUE_ERROR, temp.str().c_str());
700 }
701 //
702 // Convert any finley errors into a C++ exception
703 checkFinleyError();
704 MeshAdapter* temp=new MeshAdapter(fMesh);
705 temp->addDiracPoints(points, tags);
706 Finley_Mesh* out=temp->getMesh().get();
707 for (map<string, int>::const_iterator it=tagnamestonums.begin();it!=tagnamestonums.end();++it)
708 {
709 Finley_Mesh_addTagMap(out, it->first.c_str(), it->second);
710 }
711 Finley_ElementFile_setTagsInUse(out->Points);
712 return temp->getPtr();
713 }
714
715 Domain_ptr meshMerge(const boost::python::list& meshList)
716 {
717 Finley_Mesh* fMesh=0;
718 //
719 // extract the meshes from meshList
720 int numMsh=boost::python::extract<int>(meshList.attr("__len__")());
721 Finley_Mesh **mshes = (numMsh) ? TMPMEMALLOC(numMsh,Finley_Mesh*) : (Finley_Mesh**)NULL;
722 for (int i=0;i<numMsh;++i) {
723 AbstractContinuousDomain& meshListMember=boost::python::extract<AbstractContinuousDomain&>(meshList[i]);
724 const MeshAdapter* finley_meshListMember=static_cast<const MeshAdapter*>(&meshListMember);
725 mshes[i]=finley_meshListMember->getFinley_Mesh();
726 }
727 //
728 // merge the meshes:
729 fMesh=Finley_Mesh_merge(numMsh,mshes);
730 TMPMEMFREE(mshes);
731 //
732 // Convert any finley errors into a C++ exception
733 checkFinleyError();
734 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
735
736 return temp->getPtr();
737 }
738
739 Domain_ptr rectangle_driver(const boost::python::list& args)
740 {
741 using boost::python::extract;
742
743 // we need to convert lists to stl vectors
744 boost::python::list pypoints=extract<boost::python::list>(args[12]);
745 boost::python::list pytags=extract<boost::python::list>(args[13]);
746 int numpts=extract<int>(pypoints.attr("__len__")());
747 int numtags=extract<int>(pytags.attr("__len__")());
748 vector<double> points;
749 vector<int> tags;
750 tags.resize(numtags, -1);
751 for (int i=0;i<numpts;++i)
752 {
753 boost::python::object temp=pypoints[i];
754 int l=extract<int>(temp.attr("__len__")());
755 for (int k=0;k<l;++k)
756 {
757 points.push_back(extract<double>(temp[k]));
758 }
759 }
760 map<string, int> tagstonames;
761 int curmax=40;
762 // but which order to assign tags to names?????
763 for (int i=0;i<numtags;++i)
764 {
765 extract<int> ex_int(pytags[i]);
766 extract<string> ex_str(pytags[i]);
767 if (ex_int.check())
768 {
769 tags[i]=ex_int();
770 if (tags[i]>= curmax)
771 {
772 curmax=tags[i]+1;
773 }
774 }
775 else if (ex_str.check())
776 {
777 string s=ex_str();
778 map<string, int>::iterator it=tagstonames.find(s);
779 if (it!=tagstonames.end())
780 {
781 // we have the tag already so look it up
782 tags[i]=it->second;
783 }
784 else
785 {
786 tagstonames[s]=curmax;
787 tags[i]=curmax;
788 curmax++;
789 }
790 }
791 else
792 {
793 throw FinleyAdapterException("Error - Unable to extract tag value.");
794 }
795 }
796
797 return rectangle(extract<int>(args[0]), extract<int>(args[1]),
798 extract<int>(args[2]), extract<double>(args[3]),
799 extract<double>(args[4]), extract<int>(args[5]),
800 extract<int>(args[6]), extract<int>(args[7]),
801 extract<int>(args[8]), extract<int>(args[9]),
802 extract<int>(args[10]), extract<int>(args[11]),
803 points, tags, tagstonames);
804 }
805
806
807 Domain_ptr glueFaces(const boost::python::list& meshList,
808 double safety_factor,
809 double tolerance,
810 int optimize)
811 {
812 Finley_Mesh* fMesh=0;
813 //
814 // merge the meshes:
815 Domain_ptr merged_meshes=meshMerge(meshList);
816
817 //
818 // glue the faces:
819 const MeshAdapter* merged_finley_meshes=dynamic_cast<const MeshAdapter*>(merged_meshes.get());
820 fMesh=merged_finley_meshes->getFinley_Mesh();
821 Finley_Mesh_glueFaces(fMesh,safety_factor,tolerance,(optimize ? TRUE : FALSE));
822
823 //
824 // Convert any finley errors into a C++ exception
825 checkFinleyError();
826 return merged_meshes->getPtr();
827 }
828
829 Domain_ptr joinFaces(const boost::python::list& meshList,
830 double safety_factor,
831 double tolerance,
832 int optimize)
833 {
834 Finley_Mesh* fMesh=0;
835 //
836 // merge the meshes:
837 Domain_ptr merged_meshes=meshMerge(meshList);
838 //
839 // join the faces:
840 const MeshAdapter* merged_finley_meshes=static_cast<const MeshAdapter*>(merged_meshes.get());
841 fMesh=merged_finley_meshes->getFinley_Mesh();
842 Finley_Mesh_joinFaces(fMesh,safety_factor,tolerance, (optimize ? TRUE : FALSE));
843 //
844 // Convert any finley errors into a C++ exception
845 checkFinleyError();
846 return merged_meshes->getPtr();
847 }
848
849 // end of namespace
850
851 }
852

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26