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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26