/[escript]/trunk/finley/src/DomainFactory.cpp
ViewVC logotype

Contents of /trunk/finley/src/DomainFactory.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 1 week ago) by uqaeller
File size: 52805 byte(s)
Updated the copyright header.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17
18 #include <finley/DomainFactory.h>
19
20 #include <escript/index.h>
21 #include <escript/SubWorld.h>
22
23 #ifdef ESYS_HAVE_NETCDF
24 #ifdef NETCDF4
25 #include <ncDim.h>
26 #include <ncVar.h>
27 #include <ncFile.h>
28
29 #include <escript/NCHelper.h>
30
31 #else
32 #include <netcdfcpp.h>
33 #endif
34 #endif
35
36
37 #include <boost/python/extract.hpp>
38 #include <boost/scoped_array.hpp>
39
40 #include <sstream>
41
42 using namespace std;
43 using namespace escript;
44 namespace bp = boost::python;
45
46 #ifdef NETCDF4
47 using namespace netCDF;
48 #endif
49
50 namespace finley {
51
52 #ifdef ESYS_HAVE_NETCDF
53 #ifdef NETCDF4
54
55 // A convenience method to retrieve an integer attribute from a NetCDF file
56 template<typename T>
57 T ncReadAtt(NcFile& dataFile, const string& fName, const string& attrName)
58 {
59 NcGroupAtt attr = dataFile.getAtt(attrName.c_str());
60 if (attr.isNull()) {
61 stringstream msg;
62 msg << "loadMesh: Error retrieving integer attribute '" << attrName
63 << "' from NetCDF file '" << fName << "'";
64 throw IOError(msg.str());
65 }
66 T value;
67 attr.getValues(&value);
68 return value;
69 }
70 #else
71
72 // A convenience method to retrieve an integer attribute from a NetCDF file
73 template<typename T>
74 T ncReadAtt(NcFile* dataFile, const string& fName, const string& attrName)
75 {
76 NcAtt* attr = dataFile->get_att(attrName.c_str());
77 if (!attr) {
78 stringstream msg;
79 msg << "loadMesh: Error retrieving integer attribute '" << attrName
80 << "' from NetCDF file '" << fName << "'";
81 throw IOError(msg.str());
82 }
83 T value = (sizeof(T) > 4 ? attr->as_long(0) : attr->as_int(0));
84 delete attr;
85 return value;
86 }
87
88 #endif
89 #endif
90
91 inline void cleanupAndThrow(FinleyDomain* dom, string msg)
92 {
93 delete dom;
94 string msgPrefix("loadMesh: NetCDF operation failed - ");
95 throw IOError(msgPrefix+msg);
96 }
97
98 #ifdef NETCDF4
99
100 Domain_ptr FinleyDomain::load(const string& fileName)
101 {
102 #ifdef ESYS_HAVE_NETCDF
103 JMPI mpiInfo = makeInfo(MPI_COMM_WORLD);
104 const string fName(mpiInfo->appendRankToFileName(fileName));
105
106 // Open NetCDF file for reading
107 NcGroupAtt attr;
108 NcVar nc_var_temp;
109 NcFile dataFile;
110 // Create the NetCDF file.
111 if (!openNcFile(dataFile, fileName))
112 {
113 throw IOError("load: opening of netCDF file for input failed.");
114 }
115 // Read NetCDF integer attributes
116
117 // index_size was only introduced with 64-bit index support so fall back
118 // to 32 bits if not found.
119 int index_size;
120 try {
121 index_size = ncReadAtt<int>(dataFile, fName, "index_size");
122 } catch (IOError& e) {
123 index_size = 4;
124 }
125 // technically we could cast if reading 32-bit data on 64-bit escript
126 // but cost-benefit analysis clearly favours this implementation for now
127 if (sizeof(index_t) != index_size) {
128 throw IOError("loadMesh: size of index types at runtime differ from dump file");
129 }
130
131 int mpi_size = ncReadAtt<int>(dataFile, fName, "mpi_size");
132 int mpi_rank = ncReadAtt<int>(dataFile, fName, "mpi_rank");
133 int numDim = ncReadAtt<int>(dataFile, fName, "numDim");
134 int order = ncReadAtt<int>(dataFile, fName, "order");
135 int reduced_order = ncReadAtt<int>(dataFile, fName, "reduced_order");
136 dim_t numNodes = ncReadAtt<dim_t>(dataFile, fName, "numNodes");
137 dim_t num_Elements = ncReadAtt<dim_t>(dataFile, fName, "num_Elements");
138 dim_t num_FaceElements = ncReadAtt<dim_t>(dataFile, fName, "num_FaceElements");
139 dim_t num_ContactElements = ncReadAtt<dim_t>(dataFile, fName, "num_ContactElements");
140 dim_t num_Points = ncReadAtt<dim_t>(dataFile, fName, "num_Points");
141 int num_Elements_numNodes = ncReadAtt<int>(dataFile, fName, "num_Elements_numNodes");
142 int Elements_TypeId = ncReadAtt<int>(dataFile, fName, "Elements_TypeId");
143 int num_FaceElements_numNodes = ncReadAtt<int>(dataFile, fName, "num_FaceElements_numNodes");
144 int FaceElements_TypeId = ncReadAtt<int>(dataFile, fName, "FaceElements_TypeId");
145 int num_ContactElements_numNodes = ncReadAtt<int>(dataFile, fName, "num_ContactElements_numNodes");
146 int ContactElements_TypeId = ncReadAtt<int>(dataFile, fName, "ContactElements_TypeId");
147 int Points_TypeId = ncReadAtt<int>(dataFile, fName, "Points_TypeId");
148 int num_Tags = ncReadAtt<int>(dataFile, fName, "num_Tags");
149
150 // Verify size and rank
151 if (mpiInfo->size != mpi_size) {
152 stringstream msg;
153 msg << "loadMesh: The NetCDF file '" << fName
154 << "' can only be read on " << mpi_size
155 << " CPUs. Currently running: " << mpiInfo->size;
156 throw FinleyException(msg.str());
157 }
158 if (mpiInfo->rank != mpi_rank) {
159 stringstream msg;
160 msg << "loadMesh: The NetCDF file '" << fName
161 << "' should be read on CPU #" << mpi_rank
162 << " and NOT on #" << mpiInfo->rank;
163 throw FinleyException(msg.str());
164 }
165
166 // Read mesh name
167 if ((attr=dataFile.getAtt("Name")).isNull() ) {
168 stringstream msg;
169 msg << "loadMesh: Error retrieving mesh name from NetCDF file '"
170 << fName << "'";
171 throw IOError(msg.str());
172 }
173 string name;
174 attr.getValues(name);
175
176 // allocate mesh
177 FinleyDomain* dom = new FinleyDomain(name, numDim, mpiInfo);
178
179 // read nodes
180 NodeFile* nodes = dom->getNodes();
181 nodes->allocTable(numNodes);
182
183 try
184 {
185 // Nodes_Id
186 nc_var_temp = dataFile.getVar("Nodes_Id");
187 nc_var_temp.getVar(&nodes->Id[0]); // numNodes values read
188 // Nodes_Tag
189 nc_var_temp = dataFile.getVar("Nodes_Tag");
190 nc_var_temp.getVar(&nodes->Tag[0]); // numNodes
191 // Nodes_gDOF
192 nc_var_temp = dataFile.getVar("Nodes_gDOF");
193 nc_var_temp.getVar(&nodes->globalDegreesOfFreedom[0]); //numNodes
194 // Nodes_gNI
195 nc_var_temp = dataFile.getVar("Nodes_gNI");
196 nc_var_temp.getVar(&nodes->globalNodesIndex[0]); // numNodes
197 // Nodes_grDfI
198 nc_var_temp = dataFile.getVar("Nodes_grDfI");
199 nc_var_temp.getVar(&nodes->globalReducedDOFIndex[0]); // numNodes
200 // Nodes_grNI
201 nc_var_temp = dataFile.getVar("Nodes_grNI");
202 nc_var_temp.getVar(&nodes->globalReducedNodesIndex[0]); // numNodes
203 // Nodes_Coordinates
204 nc_var_temp = dataFile.getVar("Nodes_Coordinates");
205 nc_var_temp.getVar(&nodes->Coordinates[0]); // (numNodes, numDim)
206 }
207 catch (exceptions::NcException& e)
208 {
209 cleanupAndThrow(dom, "Read vars from file");
210 }
211 nodes->updateTagList();
212
213 // read elements
214 const_ReferenceElementSet_ptr refElements(new ReferenceElementSet(
215 (ElementTypeId)Elements_TypeId, order, reduced_order));
216 ElementFile* elements = new ElementFile(refElements, mpiInfo);
217 dom->setElements(elements);
218 elements->allocTable(num_Elements);
219 elements->minColor = 0;
220 elements->maxColor = num_Elements-1;
221 if (num_Elements > 0) {
222 try
223 {
224 // Elements_Id
225 nc_var_temp = dataFile.getVar("Elements_Id");
226 nc_var_temp.getVar(&elements->Id[0]); // num_Elements
227 // Elements_Tag
228 nc_var_temp = dataFile.getVar("Elements_Tag");
229 nc_var_temp.getVar(&elements->Tag[0]); // num_Elements
230 // Elements_Owner
231 nc_var_temp = dataFile.getVar("Elements_Owner");
232 nc_var_temp.getVar(&elements->Owner[0]); // num_Elements
233 // Elements_Color
234 nc_var_temp = dataFile.getVar("Elements_Color");
235 nc_var_temp.getVar(&elements->Color[0]); // num_Elements
236 }
237 catch (exceptions::NcException& e)
238 {
239 cleanupAndThrow(dom, "Readig element vars");
240 }
241 // Now we need to adjust maxColor
242 index_t mc = elements->Color[0];
243 for (index_t i = 1; i < num_Elements; ++i) {
244 if (mc < elements->Color[i]) {
245 mc = elements->Color[i];
246 }
247 }
248 elements->maxColor = mc;
249 // Elements_Nodes
250 int* Elements_Nodes = new int[num_Elements*num_Elements_numNodes];
251 try
252 {
253 nc_var_temp = dataFile.getVar("Elements_Nodes");
254 nc_var_temp.getVar(&Elements_Nodes[0]); // (num_Elements, num_Elements_numNodes) )
255 }
256 catch (exceptions::NcException& e)
257 {
258 delete[] Elements_Nodes;
259 cleanupAndThrow(dom, "get_var(Elements_Nodes)");
260 }
261 // Copy temp array into elements->Nodes
262 for (index_t i = 0; i < num_Elements; i++) {
263 for (int j = 0; j < num_Elements_numNodes; j++) {
264 elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
265 = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
266 }
267 }
268 delete[] Elements_Nodes;
269 } // num_Elements > 0
270 elements->updateTagList();
271
272 // get the face elements
273 const_ReferenceElementSet_ptr refFaceElements(
274 new ReferenceElementSet((ElementTypeId)FaceElements_TypeId,
275 order, reduced_order));
276 ElementFile* faces = new ElementFile(refFaceElements, mpiInfo);
277 dom->setFaceElements(faces);
278 faces->allocTable(num_FaceElements);
279 faces->minColor = 0;
280 faces->maxColor = num_FaceElements-1;
281 if (num_FaceElements > 0) {
282 try
283 {
284 // FaceElements_Id
285 nc_var_temp = dataFile.getVar("FaceElements_Id");
286 nc_var_temp.getVar(&faces->Id[0]); // num_FaceElements
287 // FaceElements_Tag
288 nc_var_temp = dataFile.getVar("FaceElements_Tag");
289 nc_var_temp.getVar(&faces->Tag[0]); // num_FaceElements
290 // FaceElements_Owner
291 nc_var_temp = dataFile.getVar("FaceElements_Owner");
292 nc_var_temp.getVar(&faces->Owner[0]); // num_FaceElements
293 // FaceElements_Color
294 nc_var_temp = dataFile.getVar("FaceElements_Color");
295 nc_var_temp.getVar(&faces->Color[0]); // num_FaceElements
296 }
297 catch (exceptions::NcException& e)
298 {
299 cleanupAndThrow(dom, "read face variables");
300 }
301 // Now we need to adjust maxColor
302 index_t mc = faces->Color[0];
303 for (index_t i = 1; i < num_FaceElements; ++i) {
304 if (mc < faces->Color[i]) {
305 mc = faces->Color[i];
306 }
307 }
308 faces->maxColor = mc;
309 // FaceElements_Nodes
310 int* FaceElements_Nodes = new int[num_FaceElements*num_FaceElements_numNodes];
311 try
312 {
313 nc_var_temp = dataFile.getVar("FaceElements_Nodes");
314 nc_var_temp.getVar(&(FaceElements_Nodes[0])); // num_FaceElements, num_FaceElements_numNodes) )
315 }
316 catch (exceptions::NcException& e)
317 {
318 delete[] FaceElements_Nodes;
319 cleanupAndThrow(dom, "read face elements");
320 }
321 if ((nc_var_temp = dataFile.getVar("FaceElements_Nodes")), nc_var_temp.isNull()) {
322 delete[] FaceElements_Nodes;
323 cleanupAndThrow(dom, "get_var(FaceElements_Nodes)");
324 }
325 nc_var_temp.getVar(&(FaceElements_Nodes[0])); // num_FaceElements, num_FaceElements_numNodes
326 // Copy temp array into faces->Nodes
327 for (index_t i = 0; i < num_FaceElements; i++) {
328 for (int j = 0; j < num_FaceElements_numNodes; j++) {
329 faces->Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
330 }
331 }
332 delete[] FaceElements_Nodes;
333 } // num_FaceElements > 0
334 faces->updateTagList();
335
336 // get the Contact elements
337 const_ReferenceElementSet_ptr refContactElements(
338 new ReferenceElementSet((ElementTypeId)ContactElements_TypeId,
339 order, reduced_order));
340 ElementFile* contacts = new ElementFile(refContactElements, mpiInfo);
341 dom->setContactElements(contacts);
342 contacts->allocTable(num_ContactElements);
343 contacts->minColor = 0;
344 contacts->maxColor = num_ContactElements-1;
345 if (num_ContactElements > 0) {
346 // ContactElements_Id
347 if (( nc_var_temp = dataFile.getVar("ContactElements_Id")), nc_var_temp.isNull() )
348 cleanupAndThrow(dom, "get_var(ContactElements_Id)");
349 nc_var_temp.getVar(&contacts->Id[0]); // num_ContactElements
350 // ContactElements_Tag
351 if (( nc_var_temp = dataFile.getVar("ContactElements_Tag")), nc_var_temp.isNull() )
352 cleanupAndThrow(dom, "get_var(ContactElements_Tag)");
353 nc_var_temp.getVar(&contacts->Tag[0]); //, num_ContactElements
354 // ContactElements_Owner
355 if (( nc_var_temp = dataFile.getVar("ContactElements_Owner")), nc_var_temp.isNull() )
356 cleanupAndThrow(dom, "get_var(ContactElements_Owner)");
357 nc_var_temp.getVar(&contacts->Owner[0]); //, num_ContactElements)
358 // ContactElements_Color
359 if (( nc_var_temp = dataFile.getVar("ContactElements_Color")), nc_var_temp.isNull() )
360 cleanupAndThrow(dom, "get_var(ContactElements_Color)");
361 nc_var_temp.getVar(&contacts->Color[0]); //, num_ContactElements
362 // Now we need to adjust maxColor
363 index_t mc = contacts->Color[0];
364 for (index_t i = 1; i < num_ContactElements; ++i) {
365 if (mc < contacts->Color[i]) {
366 mc = contacts->Color[i];
367 }
368 }
369 contacts->maxColor = mc;
370 // ContactElements_Nodes
371 int* ContactElements_Nodes = new int[num_ContactElements*num_ContactElements_numNodes];
372 if ((nc_var_temp = dataFile.getVar("ContactElements_Nodes")), nc_var_temp.isNull()) {
373 delete[] ContactElements_Nodes;
374 cleanupAndThrow(dom, "get_var(ContactElements_Nodes)");
375 }
376 nc_var_temp.getVar(&ContactElements_Nodes[0]); // num_ContactElements, num_ContactElements_numNodes
377 // Copy temp array into contacts->Nodes
378 for (index_t i = 0; i < num_ContactElements; i++) {
379 for (int j = 0; j < num_ContactElements_numNodes; j++) {
380 contacts->Nodes[INDEX2(j,i,num_ContactElements_numNodes)] = ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
381 }
382 }
383 delete[] ContactElements_Nodes;
384 } // num_ContactElements > 0
385 contacts->updateTagList();
386
387 // get the Points (nodal elements)
388 const_ReferenceElementSet_ptr refPoints(new ReferenceElementSet(
389 (ElementTypeId)Points_TypeId, order, reduced_order));
390 ElementFile* points = new ElementFile(refPoints, mpiInfo);
391 dom->setPoints(points);
392 points->allocTable(num_Points);
393 points->minColor = 0;
394 points->maxColor = num_Points-1;
395 if (num_Points > 0) {
396 // Points_Id
397 if (( nc_var_temp = dataFile.getVar("Points_Id")), nc_var_temp.isNull())
398 cleanupAndThrow(dom, "get_var(Points_Id)");
399 nc_var_temp.getVar(&points->Id[0]); // num_Points
400 // Points_Tag
401 if (( nc_var_temp = dataFile.getVar("Points_Tag")), nc_var_temp.isNull())
402 cleanupAndThrow(dom, "get_var(Points_Tag)");
403 nc_var_temp.getVar(&points->Tag[0]); // num_Points
404 // Points_Owner
405 if (( nc_var_temp = dataFile.getVar("Points_Owner")), nc_var_temp.isNull())
406 cleanupAndThrow(dom, "get_var(Points_Owner)");
407 nc_var_temp.getVar(&points->Owner[0]); // num_Points
408 // Points_Color
409 if (( nc_var_temp = dataFile.getVar("Points_Color")), nc_var_temp.isNull())
410 cleanupAndThrow(dom, "get_var(Points_Color)");
411 nc_var_temp.getVar(&points->Color[0]); // num_Points
412 // Now we need to adjust maxColor
413 index_t mc = points->Color[0];
414 for (index_t i = 1; i < num_Points; ++i) {
415 if (mc < points->Color[i]) {
416 mc = points->Color[i];
417 }
418 }
419 points->maxColor = mc;
420 // Points_Nodes
421 int* Points_Nodes = new int[num_Points];
422 if ((nc_var_temp = dataFile.getVar("Points_Nodes")), nc_var_temp.isNull()) {
423 delete[] Points_Nodes;
424 cleanupAndThrow(dom, "get_var(Points_Nodes)");
425 }
426 nc_var_temp.getVar(&Points_Nodes[0]); // num_Points
427 // Copy temp array into points->Nodes
428 for (index_t i = 0; i < num_Points; i++) {
429 points->Id[points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
430 }
431 delete[] Points_Nodes;
432 } // num_Points > 0
433 points->updateTagList();
434
435 // get the tags
436 if (num_Tags > 0) {
437 // Temp storage to gather node IDs
438 int *Tags_keys = new int[num_Tags];
439 char name_temp[4096];
440 int i;
441
442 // Tags_keys
443 if (( nc_var_temp = dataFile.getVar("Tags_keys")), nc_var_temp.isNull() ) {
444 delete[] Tags_keys;
445 cleanupAndThrow(dom, "get_var(Tags_keys)");
446 }
447 nc_var_temp.getVar(&Tags_keys[0]); // num_Tags
448 for (i=0; i<num_Tags; i++) {
449 // Retrieve tag name
450 sprintf(name_temp, "Tags_name_%d", i);
451 if ((attr=dataFile.getAtt(name_temp)), attr.isNull() ) {
452 delete[] Tags_keys;
453 stringstream msg;
454 msg << "get_att(" << name_temp << ")";
455 cleanupAndThrow(dom, msg.str());
456 }
457 std::string name;
458 attr.getValues(name);
459 // boost::scoped_array<char> name(attr->as_string(0));
460 // delete attr;
461 // dom->setTagMap(name.get(), Tags_keys[i]);
462 dom->setTagMap(name.c_str(), Tags_keys[i]);
463 }
464 delete[] Tags_keys;
465 }
466
467 // Nodes_DofDistribution
468 IndexVector first_DofComponent(mpi_size+1);
469 if ((nc_var_temp = dataFile.getVar("Nodes_DofDistribution")), nc_var_temp.isNull() ) {
470 cleanupAndThrow(dom, "get_var(Nodes_DofDistribution)");
471 }
472 nc_var_temp.getVar(&first_DofComponent[0]); // mpi_size+1
473
474 // Nodes_NodeDistribution
475 IndexVector first_NodeComponent(mpi_size+1);
476 if ((nc_var_temp = dataFile.getVar("Nodes_NodeDistribution")), nc_var_temp.isNull() ) {
477 cleanupAndThrow(dom, "get_var(Nodes_NodeDistribution)");
478 }
479 nc_var_temp.getVar(&first_NodeComponent[0]); // mpi_size+1
480 dom->createMappings(first_DofComponent, first_NodeComponent);
481
482 return dom->getPtr();
483 #else
484 throw FinleyException("loadMesh: not compiled with NetCDF. Please contact your installation manager.");
485 #endif // ESYS_HAVE_NETCDF
486 }
487
488
489 #else
490
491 Domain_ptr FinleyDomain::load(const string& fileName)
492 {
493 #ifdef ESYS_HAVE_NETCDF
494 JMPI mpiInfo = makeInfo(MPI_COMM_WORLD);
495 const string fName(mpiInfo->appendRankToFileName(fileName));
496
497 // Open NetCDF file for reading
498 NcAtt *attr;
499 NcVar *nc_var_temp;
500 // netCDF error handler
501 NcError err(NcError::silent_nonfatal);
502 // Create the NetCDF file.
503 NcFile dataFile(fName.c_str(), NcFile::ReadOnly);
504 if (!dataFile.is_valid()) {
505 stringstream msg;
506 msg << "loadMesh: Opening NetCDF file '" << fName << "' for reading failed.";
507 throw IOError(msg.str());
508 }
509
510 // Read NetCDF integer attributes
511
512 // index_size was only introduced with 64-bit index support so fall back
513 // to 32 bits if not found.
514 int index_size;
515 try {
516 index_size = ncReadAtt<int>(&dataFile, fName, "index_size");
517 } catch (IOError& e) {
518 index_size = 4;
519 }
520 // technically we could cast if reading 32-bit data on 64-bit escript
521 // but cost-benefit analysis clearly favours this implementation for now
522 if (sizeof(index_t) != index_size) {
523 throw IOError("loadMesh: size of index types at runtime differ from dump file");
524 }
525
526 int mpi_size = ncReadAtt<int>(&dataFile, fName, "mpi_size");
527 int mpi_rank = ncReadAtt<int>(&dataFile, fName, "mpi_rank");
528 int numDim = ncReadAtt<int>(&dataFile, fName, "numDim");
529 int order = ncReadAtt<int>(&dataFile, fName, "order");
530 int reduced_order = ncReadAtt<int>(&dataFile, fName, "reduced_order");
531 dim_t numNodes = ncReadAtt<dim_t>(&dataFile, fName, "numNodes");
532 dim_t num_Elements = ncReadAtt<dim_t>(&dataFile, fName, "num_Elements");
533 dim_t num_FaceElements = ncReadAtt<dim_t>(&dataFile, fName, "num_FaceElements");
534 dim_t num_ContactElements = ncReadAtt<dim_t>(&dataFile, fName, "num_ContactElements");
535 dim_t num_Points = ncReadAtt<dim_t>(&dataFile, fName, "num_Points");
536 int num_Elements_numNodes = ncReadAtt<int>(&dataFile, fName, "num_Elements_numNodes");
537 int Elements_TypeId = ncReadAtt<int>(&dataFile, fName, "Elements_TypeId");
538 int num_FaceElements_numNodes = ncReadAtt<int>(&dataFile, fName, "num_FaceElements_numNodes");
539 int FaceElements_TypeId = ncReadAtt<int>(&dataFile, fName, "FaceElements_TypeId");
540 int num_ContactElements_numNodes = ncReadAtt<int>(&dataFile, fName, "num_ContactElements_numNodes");
541 int ContactElements_TypeId = ncReadAtt<int>(&dataFile, fName, "ContactElements_TypeId");
542 int Points_TypeId = ncReadAtt<int>(&dataFile, fName, "Points_TypeId");
543 int num_Tags = ncReadAtt<int>(&dataFile, fName, "num_Tags");
544
545 // Verify size and rank
546 if (mpiInfo->size != mpi_size) {
547 stringstream msg;
548 msg << "loadMesh: The NetCDF file '" << fName
549 << "' can only be read on " << mpi_size
550 << " CPUs. Currently running: " << mpiInfo->size;
551 throw FinleyException(msg.str());
552 }
553 if (mpiInfo->rank != mpi_rank) {
554 stringstream msg;
555 msg << "loadMesh: The NetCDF file '" << fName
556 << "' should be read on CPU #" << mpi_rank
557 << " and NOT on #" << mpiInfo->rank;
558 throw FinleyException(msg.str());
559 }
560
561 // Read mesh name
562 if (! (attr=dataFile.get_att("Name")) ) {
563 stringstream msg;
564 msg << "loadMesh: Error retrieving mesh name from NetCDF file '"
565 << fName << "'";
566 throw IOError(msg.str());
567 }
568 boost::scoped_array<char> name(attr->as_string(0));
569 delete attr;
570
571 // allocate mesh
572 FinleyDomain* dom = new FinleyDomain(name.get(), numDim, mpiInfo);
573
574 // read nodes
575 NodeFile* nodes = dom->getNodes();
576 nodes->allocTable(numNodes);
577 // Nodes_Id
578 if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
579 cleanupAndThrow(dom, "get_var(Nodes_Id)");
580 if (! nc_var_temp->get(&nodes->Id[0], numNodes) )
581 cleanupAndThrow(dom, "get(Nodes_Id)");
582 // Nodes_Tag
583 if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
584 cleanupAndThrow(dom, "get_var(Nodes_Tag)");
585 if (! nc_var_temp->get(&nodes->Tag[0], numNodes) )
586 cleanupAndThrow(dom, "get(Nodes_Tag)");
587 // Nodes_gDOF
588 if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
589 cleanupAndThrow(dom, "get_var(Nodes_gDOF)");
590 if (! nc_var_temp->get(&nodes->globalDegreesOfFreedom[0], numNodes) )
591 cleanupAndThrow(dom, "get(Nodes_gDOF)");
592 // Nodes_gNI
593 if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
594 cleanupAndThrow(dom, "get_var(Nodes_gNI)");
595 if (! nc_var_temp->get(&nodes->globalNodesIndex[0], numNodes) )
596 cleanupAndThrow(dom, "get(Nodes_gNI)");
597 // Nodes_grDfI
598 if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
599 cleanupAndThrow(dom, "get_var(Nodes_grDfI)");
600 if (! nc_var_temp->get(&nodes->globalReducedDOFIndex[0], numNodes) )
601 cleanupAndThrow(dom, "get(Nodes_grDfI)");
602 // Nodes_grNI
603 if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
604 cleanupAndThrow(dom, "get_var(Nodes_grNI)");
605 if (! nc_var_temp->get(&nodes->globalReducedNodesIndex[0], numNodes) )
606 cleanupAndThrow(dom, "get(Nodes_grNI)");
607 // Nodes_Coordinates
608 if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates")))
609 cleanupAndThrow(dom, "get_var(Nodes_Coordinates)");
610 if (! nc_var_temp->get(&nodes->Coordinates[0], numNodes, numDim) )
611 cleanupAndThrow(dom, "get(Nodes_Coordinates)");
612
613 nodes->updateTagList();
614
615 // read elements
616 const_ReferenceElementSet_ptr refElements(new ReferenceElementSet(
617 (ElementTypeId)Elements_TypeId, order, reduced_order));
618 ElementFile* elements = new ElementFile(refElements, mpiInfo);
619 dom->setElements(elements);
620 elements->allocTable(num_Elements);
621 elements->minColor = 0;
622 elements->maxColor = num_Elements-1;
623 if (num_Elements > 0) {
624 // Elements_Id
625 if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
626 cleanupAndThrow(dom, "get_var(Elements_Id)");
627 if (! nc_var_temp->get(&elements->Id[0], num_Elements) )
628 cleanupAndThrow(dom, "get(Elements_Id)");
629 // Elements_Tag
630 if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
631 cleanupAndThrow(dom, "get_var(Elements_Tag)");
632 if (! nc_var_temp->get(&elements->Tag[0], num_Elements) )
633 cleanupAndThrow(dom, "get(Elements_Tag)");
634 // Elements_Owner
635 if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
636 cleanupAndThrow(dom, "get_var(Elements_Owner)");
637 if (! nc_var_temp->get(&elements->Owner[0], num_Elements) )
638 cleanupAndThrow(dom, "get(Elements_Owner)");
639 // Elements_Color
640 if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
641 cleanupAndThrow(dom, "get_var(Elements_Color)");
642 if (! nc_var_temp->get(&elements->Color[0], num_Elements) )
643 cleanupAndThrow(dom, "get(Elements_Color)");
644 // Now we need to adjust maxColor
645 index_t mc = elements->Color[0];
646 for (index_t i = 1; i < num_Elements; ++i) {
647 if (mc < elements->Color[i]) {
648 mc = elements->Color[i];
649 }
650 }
651 elements->maxColor = mc;
652 // Elements_Nodes
653 int* Elements_Nodes = new int[num_Elements*num_Elements_numNodes];
654 if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
655 delete[] Elements_Nodes;
656 cleanupAndThrow(dom, "get_var(Elements_Nodes)");
657 }
658 if (! nc_var_temp->get(&Elements_Nodes[0], num_Elements, num_Elements_numNodes) ) {
659 delete[] Elements_Nodes;
660 cleanupAndThrow(dom, "get(Elements_Nodes)");
661 }
662
663 // Copy temp array into elements->Nodes
664 for (index_t i = 0; i < num_Elements; i++) {
665 for (int j = 0; j < num_Elements_numNodes; j++) {
666 elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
667 = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
668 }
669 }
670 delete[] Elements_Nodes;
671 } // num_Elements > 0
672 elements->updateTagList();
673
674 // get the face elements
675 const_ReferenceElementSet_ptr refFaceElements(
676 new ReferenceElementSet((ElementTypeId)FaceElements_TypeId,
677 order, reduced_order));
678 ElementFile* faces = new ElementFile(refFaceElements, mpiInfo);
679 dom->setFaceElements(faces);
680 faces->allocTable(num_FaceElements);
681 faces->minColor = 0;
682 faces->maxColor = num_FaceElements-1;
683 if (num_FaceElements > 0) {
684 // FaceElements_Id
685 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
686 cleanupAndThrow(dom, "get_var(FaceElements_Id)");
687 if (! nc_var_temp->get(&faces->Id[0], num_FaceElements) )
688 cleanupAndThrow(dom, "get(FaceElements_Id)");
689 // FaceElements_Tag
690 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
691 cleanupAndThrow(dom, "get_var(FaceElements_Tag)");
692 if (! nc_var_temp->get(&faces->Tag[0], num_FaceElements) )
693 cleanupAndThrow(dom, "get(FaceElements_Tag)");
694 // FaceElements_Owner
695 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
696 cleanupAndThrow(dom, "get_var(FaceElements_Owner)");
697 if (! nc_var_temp->get(&faces->Owner[0], num_FaceElements) )
698 cleanupAndThrow(dom, "get(FaceElements_Owner)");
699 // FaceElements_Color
700 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
701 cleanupAndThrow(dom, "get_var(FaceElements_Color)");
702 if (! nc_var_temp->get(&faces->Color[0], num_FaceElements) )
703 cleanupAndThrow(dom, "get(FaceElements_Color)");
704 // Now we need to adjust maxColor
705 index_t mc = faces->Color[0];
706 for (index_t i = 1; i < num_FaceElements; ++i) {
707 if (mc < faces->Color[i]) {
708 mc = faces->Color[i];
709 }
710 }
711 faces->maxColor = mc;
712 // FaceElements_Nodes
713 int* FaceElements_Nodes = new int[num_FaceElements*num_FaceElements_numNodes];
714 if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
715 delete[] FaceElements_Nodes;
716 cleanupAndThrow(dom, "get_var(FaceElements_Nodes)");
717 }
718 if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
719 delete[] FaceElements_Nodes;
720 cleanupAndThrow(dom, "get(FaceElements_Nodes)");
721 }
722 // Copy temp array into faces->Nodes
723 for (index_t i = 0; i < num_FaceElements; i++) {
724 for (int j = 0; j < num_FaceElements_numNodes; j++) {
725 faces->Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
726 }
727 }
728 delete[] FaceElements_Nodes;
729 } // num_FaceElements > 0
730 faces->updateTagList();
731
732 // get the Contact elements
733 const_ReferenceElementSet_ptr refContactElements(
734 new ReferenceElementSet((ElementTypeId)ContactElements_TypeId,
735 order, reduced_order));
736 ElementFile* contacts = new ElementFile(refContactElements, mpiInfo);
737 dom->setContactElements(contacts);
738 contacts->allocTable(num_ContactElements);
739 contacts->minColor = 0;
740 contacts->maxColor = num_ContactElements-1;
741 if (num_ContactElements > 0) {
742 // ContactElements_Id
743 if (! ( nc_var_temp = dataFile.get_var("ContactElements_Id")) )
744 cleanupAndThrow(dom, "get_var(ContactElements_Id)");
745 if (! nc_var_temp->get(&contacts->Id[0], num_ContactElements) )
746 cleanupAndThrow(dom, "get(ContactElements_Id)");
747 // ContactElements_Tag
748 if (! ( nc_var_temp = dataFile.get_var("ContactElements_Tag")) )
749 cleanupAndThrow(dom, "get_var(ContactElements_Tag)");
750 if (! nc_var_temp->get(&contacts->Tag[0], num_ContactElements) )
751 cleanupAndThrow(dom, "get(ContactElements_Tag)");
752 // ContactElements_Owner
753 if (! ( nc_var_temp = dataFile.get_var("ContactElements_Owner")) )
754 cleanupAndThrow(dom, "get_var(ContactElements_Owner)");
755 if (! nc_var_temp->get(&contacts->Owner[0], num_ContactElements) )
756 cleanupAndThrow(dom, "get(ContactElements_Owner)");
757 // ContactElements_Color
758 if (! ( nc_var_temp = dataFile.get_var("ContactElements_Color")) )
759 cleanupAndThrow(dom, "get_var(ContactElements_Color)");
760 if (! nc_var_temp->get(&contacts->Color[0], num_ContactElements) )
761 cleanupAndThrow(dom, "get(ContactElements_Color)");
762 // Now we need to adjust maxColor
763 index_t mc = contacts->Color[0];
764 for (index_t i = 1; i < num_ContactElements; ++i) {
765 if (mc < contacts->Color[i]) {
766 mc = contacts->Color[i];
767 }
768 }
769 contacts->maxColor = mc;
770 // ContactElements_Nodes
771 int* ContactElements_Nodes = new int[num_ContactElements*num_ContactElements_numNodes];
772 if (!(nc_var_temp = dataFile.get_var("ContactElements_Nodes"))) {
773 delete[] ContactElements_Nodes;
774 cleanupAndThrow(dom, "get_var(ContactElements_Nodes)");
775 }
776 if (! nc_var_temp->get(&ContactElements_Nodes[0], num_ContactElements, num_ContactElements_numNodes) ) {
777 delete[] ContactElements_Nodes;
778 cleanupAndThrow(dom, "get(ContactElements_Nodes)");
779 }
780 // Copy temp array into contacts->Nodes
781 for (index_t i = 0; i < num_ContactElements; i++) {
782 for (int j = 0; j < num_ContactElements_numNodes; j++) {
783 contacts->Nodes[INDEX2(j,i,num_ContactElements_numNodes)] = ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)];
784 }
785 }
786 delete[] ContactElements_Nodes;
787 } // num_ContactElements > 0
788 contacts->updateTagList();
789
790 // get the Points (nodal elements)
791 const_ReferenceElementSet_ptr refPoints(new ReferenceElementSet(
792 (ElementTypeId)Points_TypeId, order, reduced_order));
793 ElementFile* points = new ElementFile(refPoints, mpiInfo);
794 dom->setPoints(points);
795 points->allocTable(num_Points);
796 points->minColor = 0;
797 points->maxColor = num_Points-1;
798 if (num_Points > 0) {
799 // Points_Id
800 if (! ( nc_var_temp = dataFile.get_var("Points_Id")))
801 cleanupAndThrow(dom, "get_var(Points_Id)");
802 if (! nc_var_temp->get(&points->Id[0], num_Points))
803 cleanupAndThrow(dom, "get(Points_Id)");
804 // Points_Tag
805 if (! ( nc_var_temp = dataFile.get_var("Points_Tag")))
806 cleanupAndThrow(dom, "get_var(Points_Tag)");
807 if (! nc_var_temp->get(&points->Tag[0], num_Points))
808 cleanupAndThrow(dom, "get(Points_Tag)");
809 // Points_Owner
810 if (! ( nc_var_temp = dataFile.get_var("Points_Owner")))
811 cleanupAndThrow(dom, "get_var(Points_Owner)");
812 if (!nc_var_temp->get(&points->Owner[0], num_Points))
813 cleanupAndThrow(dom, "get(Points_Owner)");
814 // Points_Color
815 if (! ( nc_var_temp = dataFile.get_var("Points_Color")))
816 cleanupAndThrow(dom, "get_var(Points_Color)");
817 if (!nc_var_temp->get(&points->Color[0], num_Points))
818 cleanupAndThrow(dom, "get(Points_Color)");
819 // Now we need to adjust maxColor
820 index_t mc = points->Color[0];
821 for (index_t i = 1; i < num_Points; ++i) {
822 if (mc < points->Color[i]) {
823 mc = points->Color[i];
824 }
825 }
826 points->maxColor = mc;
827 // Points_Nodes
828 int* Points_Nodes = new int[num_Points];
829 if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
830 delete[] Points_Nodes;
831 cleanupAndThrow(dom, "get_var(Points_Nodes)");
832 }
833 if (! nc_var_temp->get(&Points_Nodes[0], num_Points) ) {
834 delete[] Points_Nodes;
835 cleanupAndThrow(dom, "get(Points_Nodes)");
836 }
837 // Copy temp array into points->Nodes
838 for (index_t i = 0; i < num_Points; i++) {
839 points->Id[points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
840 }
841 delete[] Points_Nodes;
842 } // num_Points > 0
843 points->updateTagList();
844
845 // get the tags
846 if (num_Tags > 0) {
847 // Temp storage to gather node IDs
848 int *Tags_keys = new int[num_Tags];
849 char name_temp[4096];
850 int i;
851
852 // Tags_keys
853 if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) ) {
854 delete[] Tags_keys;
855 cleanupAndThrow(dom, "get_var(Tags_keys)");
856 }
857 if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
858 delete[] Tags_keys;
859 cleanupAndThrow(dom, "get(Tags_keys)");
860 }
861 for (i=0; i<num_Tags; i++) {
862 // Retrieve tag name
863 sprintf(name_temp, "Tags_name_%d", i);
864 if (! (attr=dataFile.get_att(name_temp)) ) {
865 delete[] Tags_keys;
866 stringstream msg;
867 msg << "get_att(" << name_temp << ")";
868 cleanupAndThrow(dom, msg.str());
869 }
870 boost::scoped_array<char> name(attr->as_string(0));
871 delete attr;
872 dom->setTagMap(name.get(), Tags_keys[i]);
873 }
874 delete[] Tags_keys;
875 }
876
877 // Nodes_DofDistribution
878 IndexVector first_DofComponent(mpi_size+1);
879 if (! (nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) ) {
880 cleanupAndThrow(dom, "get_var(Nodes_DofDistribution)");
881 }
882 if (!nc_var_temp->get(&first_DofComponent[0], mpi_size+1)) {
883 cleanupAndThrow(dom, "get(Nodes_DofDistribution)");
884 }
885
886 // Nodes_NodeDistribution
887 IndexVector first_NodeComponent(mpi_size+1);
888 if (! (nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) ) {
889 cleanupAndThrow(dom, "get_var(Nodes_NodeDistribution)");
890 }
891 if (!nc_var_temp->get(&first_NodeComponent[0], mpi_size+1)) {
892 cleanupAndThrow(dom, "get(Nodes_NodeDistribution)");
893 }
894 dom->createMappings(first_DofComponent, first_NodeComponent);
895
896 return dom->getPtr();
897 #else
898 throw FinleyException("loadMesh: not compiled with NetCDF. Please contact your installation manager.");
899 #endif // ESYS_HAVE_NETCDF
900 }
901
902 #endif
903
904 Domain_ptr readMesh_driver(const bp::list& args)
905 {
906 int l = len(args);
907 if (l < 7) {
908 throw ValueError("Insufficient arguments to readMesh_driver");
909 }
910 string fileName = bp::extract<string>(args[0])();
911 int integrationOrder = bp::extract<int>(args[1])();
912 int reducedIntegrationOrder = bp::extract<int>(args[2])();
913 bool optimize = bp::extract<bool>(args[3])();
914 vector<double> points;
915 vector<int> tags;
916
917 // we need to convert lists to stl vectors
918 bp::list pypoints = bp::extract<bp::list>(args[4]);
919 bp::list pytags = bp::extract<bp::list>(args[5]);
920 int numpts = bp::extract<int>(pypoints.attr("__len__")());
921 int numtags = bp::extract<int>(pytags.attr("__len__")());
922
923 bp::object pworld = args[6];
924 JMPI info;
925 if (!pworld.is_none()) {
926 bp::extract<SubWorld_ptr> ex(pworld);
927 if (!ex.check()) {
928 throw ValueError("Invalid escriptWorld parameter.");
929 }
930 info = ex()->getMPI();
931 } else {
932 info = makeInfo(MPI_COMM_WORLD);
933 }
934 Domain_ptr dom(FinleyDomain::read(info, fileName, integrationOrder,
935 reducedIntegrationOrder, optimize));
936
937 FinleyDomain* fd = dynamic_cast<FinleyDomain*>(dom.get());
938
939 for (int i = 0; i < numpts; ++i) {
940 bp::object temp = pypoints[i];
941 int l = bp::extract<int>(temp.attr("__len__")());
942 for (int k = 0; k < l; ++k) {
943 points.push_back(bp::extract<double>(temp[k]));
944 }
945 }
946 // bricks use up to 200 but the existing tag check will find that
947 int curmax = 40;
948 const TagMap& tagmap = fd->getTagMap();
949 // first we work out what tags are already in use
950 for (TagMap::const_iterator it = tagmap.begin(); it != tagmap.end(); ++it) {
951 if (it->second > curmax) {
952 curmax = it->second+1;
953 }
954 }
955
956 tags.resize(numtags, -1);
957 for (int i = 0; i < numtags; ++i) {
958 bp::extract<int> ex_int(pytags[i]);
959 bp::extract<string> ex_str(pytags[i]);
960 if (ex_int.check()) {
961 tags[i] = ex_int();
962 if (tags[i] >= curmax) {
963 curmax = tags[i]+1;
964 }
965 } else if (ex_str.check()) {
966 string s = ex_str();
967 TagMap::const_iterator it = tagmap.find(s);
968 if (it != tagmap.end()) {
969 // we have the tag already so look it up
970 tags[i] = it->second;
971 } else {
972 fd->setTagMap(s, curmax);
973 tags[i] = curmax;
974 curmax++;
975 }
976 } else {
977 throw FinleyException("Unable to extract tag value.");
978 }
979 }
980 // now we need to add the dirac points
981 fd->addDiracPoints(points, tags);
982 return dom;
983 }
984
985 Domain_ptr readGmsh_driver(const bp::list& args)
986 {
987 int l = len(args);
988 if (l < 7) {
989 throw ValueError("Insufficient arguments to readMesh_driver");
990 }
991 string fileName = bp::extract<string>(args[0])();
992 int numDim = bp::extract<int>(args[1])();
993 int integrationOrder = bp::extract<int>(args[2])();
994 int reducedIntegrationOrder = bp::extract<int>(args[3])();
995 bool optimize = bp::extract<bool>(args[4])();
996 bool useMacroElements = bp::extract<bool>(args[5])();
997 vector<double> points;
998 vector<int> tags;
999
1000 // we need to convert lists to stl vectors
1001 bp::list pypoints = bp::extract<bp::list>(args[6]);
1002 bp::list pytags = bp::extract<bp::list>(args[7]);
1003 int numpts = bp::extract<int>(pypoints.attr("__len__")());
1004 int numtags = bp::extract<int>(pytags.attr("__len__")());
1005 bp::object pworld = args[8];
1006 JMPI info;
1007 if (!pworld.is_none()) {
1008 bp::extract<SubWorld_ptr> ex(pworld);
1009 if (!ex.check()) {
1010 throw ValueError("Invalid escriptWorld parameter.");
1011 }
1012 info = ex()->getMPI();
1013 } else {
1014 info = makeInfo(MPI_COMM_WORLD);
1015 }
1016 Domain_ptr dom(FinleyDomain::readGmsh(info, fileName, numDim,
1017 integrationOrder, reducedIntegrationOrder,
1018 optimize, useMacroElements));
1019 FinleyDomain* fd = dynamic_cast<FinleyDomain*>(dom.get());
1020
1021 for (int i = 0; i < numpts; ++i) {
1022 bp::object temp = pypoints[i];
1023 int l = bp::extract<int>(temp.attr("__len__")());
1024 for (int k = 0; k < l; ++k) {
1025 points.push_back(bp::extract<double>(temp[k]));
1026 }
1027 }
1028 int curmax = 40; // bricks use up to 30
1029 const TagMap& tagmap = fd->getTagMap();
1030 // first we work out what tags are already in use
1031 for (TagMap::const_iterator it = tagmap.begin(); it != tagmap.end(); ++it) {
1032 if (it->second > curmax) {
1033 curmax = it->second+1;
1034 }
1035 }
1036
1037 tags.resize(numtags, -1);
1038 for (int i = 0; i < numtags; ++i) {
1039 bp::extract<int> ex_int(pytags[i]);
1040 bp::extract<string> ex_str(pytags[i]);
1041 if (ex_int.check()) {
1042 tags[i] = ex_int();
1043 if (tags[i] >= curmax) {
1044 curmax = tags[i]+1;
1045 }
1046 } else if (ex_str.check()) {
1047 string s = ex_str();
1048 TagMap::const_iterator it = tagmap.find(s);
1049 if (it != tagmap.end()) {
1050 // we have the tag already so look it up
1051 tags[i] = it->second;
1052 } else {
1053 fd->setTagMap(s, curmax);
1054 tags[i] = curmax;
1055 curmax++;
1056 }
1057 } else {
1058 throw FinleyException("Unable to extract tag value");
1059 }
1060 }
1061 // now we need to add the dirac points
1062 fd->addDiracPoints(points, tags);
1063 return dom;
1064 }
1065
1066 Domain_ptr brick(JMPI info, dim_t n0, dim_t n1, dim_t n2, int order,
1067 double l0, double l1, double l2,
1068 bool periodic0, bool periodic1, bool periodic2,
1069 int integrationOrder, int reducedIntegrationOrder,
1070 bool useElementsOnFace, bool useFullElementOrder,
1071 bool optimize, const std::vector<double>& points,
1072 const std::vector<int>& tags,
1073 const std::map<std::string, int>& tagNamesToNums)
1074 {
1075 Domain_ptr dom;
1076 if (order == 1) {
1077 dom = FinleyDomain::createHex8(n0, n1, n2, l0, l1, l2, periodic0,
1078 periodic1, periodic2, integrationOrder,
1079 reducedIntegrationOrder, useElementsOnFace, optimize, info);
1080 } else if (order == 2) {
1081 dom = FinleyDomain::createHex20(n0, n1, n2, l0, l1, l2, periodic0,
1082 periodic1, periodic2, integrationOrder,
1083 reducedIntegrationOrder, useElementsOnFace,
1084 useFullElementOrder, false, optimize, info);
1085 } else if (order == -1) {
1086 dom = FinleyDomain::createHex20(n0, n1, n2, l0, l1, l2, periodic0,
1087 periodic1, periodic2, integrationOrder,
1088 reducedIntegrationOrder, useElementsOnFace,
1089 useFullElementOrder, true, optimize, info);
1090 } else {
1091 stringstream message;
1092 message << "Illegal interpolation order " << order;
1093 throw ValueError(message.str());
1094 }
1095
1096 FinleyDomain* fd = dynamic_cast<FinleyDomain*>(dom.get());
1097 fd->addDiracPoints(points, tags);
1098 for (TagMap::const_iterator it = tagNamesToNums.begin(); it != tagNamesToNums.end(); ++it) {
1099 fd->setTagMap(it->first, it->second);
1100 }
1101 fd->getPoints()->updateTagList();
1102 return dom;
1103 }
1104
1105 Domain_ptr brick_driver(const bp::list& args)
1106 {
1107 // we need to convert lists to stl vectors
1108 bp::list pypoints = bp::extract<bp::list>(args[15]);
1109 bp::list pytags = bp::extract<bp::list>(args[16]);
1110 int numpts = bp::extract<int>(pypoints.attr("__len__")());
1111 int numtags = bp::extract<int>(pytags.attr("__len__")());
1112 vector<double> points;
1113 vector<int> tags;
1114 tags.resize(numtags, -1);
1115 for (int i = 0; i < numpts; ++i) {
1116 bp::object temp = pypoints[i];
1117 int l = bp::extract<int>(temp.attr("__len__")());
1118 for (int k = 0; k < l; ++k) {
1119 points.push_back(bp::extract<double>(temp[k]));
1120 }
1121 }
1122 map<string, int> namestonums;
1123 int curmax = 40; // bricks use up to 30
1124 for (int i = 0; i < numtags; ++i) {
1125 bp::extract<int> ex_int(pytags[i]);
1126 bp::extract<string> ex_str(pytags[i]);
1127 if (ex_int.check()) {
1128 tags[i] = ex_int();
1129 if (tags[i] >= curmax) {
1130 curmax = tags[i]+1;
1131 }
1132 } else if (ex_str.check()) {
1133 string s = ex_str();
1134 TagMap::iterator it = namestonums.find(s);
1135 if (it != namestonums.end()) {
1136 // we have the tag already so look it up
1137 tags[i] = it->second;
1138 } else {
1139 namestonums[s] = curmax;
1140 tags[i] = curmax;
1141 curmax++;
1142 }
1143 } else {
1144 throw FinleyException("Unable to extract tag value.");
1145 }
1146 }
1147 bp::object pworld = args[17];
1148 JMPI info;
1149 if (!pworld.is_none()) {
1150 bp::extract<SubWorld_ptr> ex(pworld);
1151 if (!ex.check()) {
1152 throw ValueError("Invalid escriptWorld parameter.");
1153 }
1154 info = ex()->getMPI();
1155 } else {
1156 info = makeInfo(MPI_COMM_WORLD);
1157 }
1158 return brick(info, static_cast<dim_t>(bp::extract<float>(args[0])),
1159 static_cast<dim_t>(bp::extract<float>(args[1])),
1160 static_cast<dim_t>(bp::extract<float>(args[2])),
1161 bp::extract<int>(args[3]), bp::extract<double>(args[4]),
1162 bp::extract<double>(args[5]), bp::extract<double>(args[6]),
1163 bp::extract<int>(args[7]), bp::extract<int>(args[8]),
1164 bp::extract<int>(args[9]), bp::extract<int>(args[10]),
1165 bp::extract<int>(args[11]), bp::extract<int>(args[12]),
1166 bp::extract<int>(args[13]), bp::extract<int>(args[14]),
1167 points, tags, namestonums);
1168 }
1169
1170 Domain_ptr rectangle(JMPI info, dim_t n0, dim_t n1, int order,
1171 double l0, double l1, bool periodic0, bool periodic1,
1172 int integrationOrder, int reducedIntegrationOrder,
1173 bool useElementsOnFace, bool useFullElementOrder,
1174 bool optimize, const vector<double>& points,
1175 const vector<int>& tags,
1176 const std::map<std::string, int>& tagNamesToNums)
1177 {
1178 Domain_ptr dom;
1179 if (order == 1) {
1180 dom = FinleyDomain::createRec4(n0, n1, l0, l1, periodic0, periodic1,
1181 integrationOrder, reducedIntegrationOrder,
1182 useElementsOnFace, optimize, info);
1183 } else if (order == 2) {
1184 dom = FinleyDomain::createRec8(n0, n1, l0, l1, periodic0, periodic1,
1185 integrationOrder, reducedIntegrationOrder,
1186 useElementsOnFace,useFullElementOrder, false, optimize, info);
1187 } else if (order == -1) {
1188 dom = FinleyDomain::createRec8(n0, n1, l0, l1, periodic0, periodic1,
1189 integrationOrder, reducedIntegrationOrder,
1190 useElementsOnFace, useFullElementOrder, true, optimize, info);
1191 } else {
1192 stringstream message;
1193 message << "Illegal interpolation order " << order;
1194 throw ValueError(message.str());
1195 }
1196
1197 FinleyDomain* fd = dynamic_cast<FinleyDomain*>(dom.get());
1198 fd->addDiracPoints(points, tags);
1199 for (TagMap::const_iterator it = tagNamesToNums.begin(); it != tagNamesToNums.end(); ++it)
1200 {
1201 fd->setTagMap(it->first, it->second);
1202 }
1203 fd->getPoints()->updateTagList();
1204 return dom;
1205 }
1206
1207 Domain_ptr rectangle_driver(const bp::list& args)
1208 {
1209 // we need to convert lists to stl vectors
1210 bp::list pypoints = bp::extract<bp::list>(args[12]);
1211 bp::list pytags = bp::extract<bp::list>(args[13]);
1212 int numpts = bp::extract<int>(pypoints.attr("__len__")());
1213 int numtags = bp::extract<int>(pytags.attr("__len__")());
1214 vector<double> points;
1215 vector<int> tags;
1216 tags.resize(numtags, -1);
1217 for (int i = 0; i < numpts; ++i) {
1218 bp::object temp = pypoints[i];
1219 int l = bp::extract<int>(temp.attr("__len__")());
1220 for (int k = 0; k < l; ++k) {
1221 points.push_back(bp::extract<double>(temp[k]));
1222 }
1223 }
1224 TagMap tagstonames;
1225 int curmax = 40;
1226 // but which order to assign tags to names?????
1227 for (int i = 0; i < numtags; ++i) {
1228 bp::extract<int> ex_int(pytags[i]);
1229 bp::extract<string> ex_str(pytags[i]);
1230 if (ex_int.check()) {
1231 tags[i] = ex_int();
1232 if (tags[i] >= curmax) {
1233 curmax = tags[i]+1;
1234 }
1235 } else if (ex_str.check()) {
1236 string s = ex_str();
1237 TagMap::iterator it = tagstonames.find(s);
1238 if (it != tagstonames.end()) {
1239 // we have the tag already so look it up
1240 tags[i] = it->second;
1241 } else {
1242 tagstonames[s] = curmax;
1243 tags[i] = curmax;
1244 curmax++;
1245 }
1246 } else {
1247 throw FinleyException("Unable to extract tag value.");
1248 }
1249 }
1250 bp::object pworld = args[14];
1251 JMPI info;
1252 if (!pworld.is_none()) {
1253 bp::extract<SubWorld_ptr> ex(pworld);
1254 if (!ex.check()) {
1255 throw ValueError("Invalid escriptWorld parameter.");
1256 }
1257 info = ex()->getMPI();
1258 } else {
1259 info = makeInfo(MPI_COMM_WORLD);
1260 }
1261
1262 return rectangle(info, static_cast<dim_t>(bp::extract<float>(args[0])),
1263 static_cast<dim_t>(bp::extract<float>(args[1])),
1264 bp::extract<int>(args[2]), bp::extract<double>(args[3]),
1265 bp::extract<double>(args[4]), bp::extract<int>(args[5]),
1266 bp::extract<int>(args[6]), bp::extract<int>(args[7]),
1267 bp::extract<int>(args[8]), bp::extract<int>(args[9]),
1268 bp::extract<int>(args[10]), bp::extract<int>(args[11]),
1269 points, tags, tagstonames);
1270 }
1271
1272 Domain_ptr meshMerge(const bp::list& meshList)
1273 {
1274 // extract the meshes from meshList
1275 int num = bp::extract<int>(meshList.attr("__len__")());
1276 vector<const FinleyDomain*> meshes(num);
1277 for (int i = 0; i < num; ++i) {
1278 AbstractContinuousDomain& meshListMember = bp::extract<AbstractContinuousDomain&>(meshList[i]);
1279 meshes[i] = dynamic_cast<const FinleyDomain*>(&meshListMember);
1280 }
1281
1282 // merge the meshes
1283 FinleyDomain* dom = FinleyDomain::merge(meshes);
1284
1285 return dom->getPtr();
1286 }
1287
1288 Domain_ptr glueFaces(const bp::list& meshList, double safetyFactor,
1289 double tolerance, bool optimize)
1290 {
1291 // merge the meshes
1292 Domain_ptr merged_meshes = meshMerge(meshList);
1293
1294 // glue the faces
1295 FinleyDomain* merged = dynamic_cast<FinleyDomain*>(merged_meshes.get());
1296 merged->glueFaces(safetyFactor, tolerance, optimize);
1297 return merged_meshes;
1298 }
1299
1300 Domain_ptr joinFaces(const bp::list& meshList, double safetyFactor,
1301 double tolerance, bool optimize)
1302 {
1303 // merge the meshes
1304 Domain_ptr merged_meshes = meshMerge(meshList);
1305
1306 // join the faces
1307 FinleyDomain* merged = dynamic_cast<FinleyDomain*>(merged_meshes.get());
1308 merged->joinFaces(safetyFactor, tolerance, optimize);
1309 return merged_meshes;
1310 }
1311
1312 } // namespace finley

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26