/[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 6144 - (show annotations)
Wed Apr 6 05:25:13 2016 UTC (2 years, 7 months ago) by caltinay
File size: 37984 byte(s)
last round of namespacing defines.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26