/[escript]/branches/trilinos_from_5897/dudley/src/CPPAdapter/MeshAdapterFactory.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/dudley/src/CPPAdapter/MeshAdapterFactory.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6079 - (show annotations)
Mon Mar 21 12:22:38 2016 UTC (2 years, 11 months ago) by caltinay
File size: 20409 byte(s)
Big commit - making dudley much more like finley to make it more
managable. Fixed quite a few issues that had been fixed in finley.
Disposed of all ReducedNode/ReducedDOF entities that dudley never supported.
Compiles and passes tests.

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 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 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 <dudley/Dudley.h>
19 #include <dudley/Mesh.h>
20 #include <dudley/TriangularMesh.h>
21
22 #include <escript/SubWorld.h>
23
24 #ifdef USE_NETCDF
25 #include <netcdfcpp.h>
26 #endif
27
28 #include <boost/python/extract.hpp>
29 #include <boost/scoped_array.hpp>
30
31 #include <sstream>
32
33 using namespace std;
34 using namespace escript;
35
36 namespace dudley {
37
38 #ifdef USE_NETCDF
39 // A convenience method to retrieve an integer attribute from a NetCDF file
40 template<typename T>
41 T ncReadAtt(NcFile* dataFile, const string& fName, const string& attrName)
42 {
43 NcAtt* attr = dataFile->get_att(attrName.c_str());
44 if (!attr) {
45 stringstream msg;
46 msg << "loadMesh: Error retrieving integer attribute '" << attrName
47 << "' from NetCDF file '" << fName << "'";
48 throw escript::IOError(msg.str());
49 }
50 T value = (sizeof(T) > 4 ? attr->as_long(0) : attr->as_int(0));
51 delete attr;
52 return value;
53 }
54 #endif
55
56 inline void cleanupAndThrow(Mesh* mesh, string msg)
57 {
58 delete mesh;
59 string msgPrefix("loadMesh: NetCDF operation failed - ");
60 throw escript::IOError(msgPrefix+msg);
61 }
62
63 Domain_ptr loadMesh(const std::string& fileName)
64 {
65 #ifdef USE_NETCDF
66 escript::JMPI mpi_info = escript::makeInfo(MPI_COMM_WORLD);
67 const string fName(mpi_info->appendRankToFileName(fileName));
68
69 // Open NetCDF file for reading
70 NcAtt *attr;
71 NcVar *nc_var_temp;
72 // netCDF error handler
73 NcError err(NcError::silent_nonfatal);
74 // Create the NetCDF file.
75 NcFile dataFile(fName.c_str(), NcFile::ReadOnly);
76 if (!dataFile.is_valid()) {
77 stringstream msg;
78 msg << "loadMesh: Opening NetCDF file '" << fName << "' for reading failed.";
79 throw escript::IOError(msg.str());
80 }
81
82 // Read NetCDF integer attributes
83
84 // index_size was only introduced with 64-bit index support so fall back
85 // to 32 bits if not found.
86 int index_size;
87 try {
88 index_size = ncReadAtt<int>(&dataFile, fName, "index_size");
89 } catch (escript::IOError& e) {
90 index_size = 4;
91 }
92 // technically we could cast if reading 32-bit data on 64-bit escript
93 // but cost-benefit analysis clearly favours this implementation for now
94 if (sizeof(index_t) != index_size) {
95 throw escript::IOError("loadMesh: size of index types at runtime differ from dump file");
96 }
97
98 int mpi_size = ncReadAtt<int>(&dataFile, fName, "mpi_size");
99 int mpi_rank = ncReadAtt<int>(&dataFile, fName, "mpi_rank");
100 int numDim = ncReadAtt<int>(&dataFile, fName, "numDim");
101 dim_t numNodes = ncReadAtt<dim_t>(&dataFile, fName, "numNodes");
102 dim_t num_Elements = ncReadAtt<dim_t>(&dataFile, fName, "num_Elements");
103 dim_t num_FaceElements = ncReadAtt<dim_t>(&dataFile, fName, "num_FaceElements");
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 Points_TypeId = ncReadAtt<int>(&dataFile, fName, "Points_TypeId");
110 int num_Tags = ncReadAtt<int>(&dataFile, fName, "num_Tags");
111
112 // Verify size and rank
113 if (mpi_info->size != mpi_size) {
114 stringstream msg;
115 msg << "loadMesh: The NetCDF file '" << fName
116 << "' can only be read on " << mpi_size
117 << " CPUs. Currently running: " << mpi_info->size;
118 throw DudleyException(msg.str());
119 }
120 if (mpi_info->rank != mpi_rank) {
121 stringstream msg;
122 msg << "loadMesh: The NetCDF file '" << fName
123 << "' should be read on CPU #" << mpi_rank
124 << " and NOT on #" << mpi_info->rank;
125 throw DudleyException(msg.str());
126 }
127
128 // Read mesh name
129 if (! (attr=dataFile.get_att("Name")) ) {
130 stringstream msg;
131 msg << "loadMesh: Error retrieving mesh name from NetCDF file '"
132 << fName << "'";
133 throw escript::IOError(msg.str());
134 }
135 boost::scoped_array<char> name(attr->as_string(0));
136 delete attr;
137
138 // allocate mesh
139 Mesh* mesh = new Mesh(name.get(), numDim, mpi_info);
140 // read nodes
141 mesh->Nodes->allocTable(numNodes);
142 // Nodes_Id
143 if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
144 cleanupAndThrow(mesh, "get_var(Nodes_Id)");
145 if (! nc_var_temp->get(&mesh->Nodes->Id[0], numNodes) )
146 cleanupAndThrow(mesh, "get(Nodes_Id)");
147 // Nodes_Tag
148 if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
149 cleanupAndThrow(mesh, "get_var(Nodes_Tag)");
150 if (! nc_var_temp->get(&mesh->Nodes->Tag[0], numNodes) )
151 cleanupAndThrow(mesh, "get(Nodes_Tag)");
152 // Nodes_gDOF
153 if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
154 cleanupAndThrow(mesh, "get_var(Nodes_gDOF)");
155 if (! nc_var_temp->get(&mesh->Nodes->globalDegreesOfFreedom[0], numNodes) )
156 cleanupAndThrow(mesh, "get(Nodes_gDOF)");
157 // Nodes_gNI
158 if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
159 cleanupAndThrow(mesh, "get_var(Nodes_gNI)");
160 if (! nc_var_temp->get(&mesh->Nodes->globalNodesIndex[0], numNodes) )
161 cleanupAndThrow(mesh, "get(Nodes_gNI)");
162 // Nodes_Coordinates
163 if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates")))
164 cleanupAndThrow(mesh, "get_var(Nodes_Coordinates)");
165 if (! nc_var_temp->get(&(mesh->Nodes->Coordinates[0]), numNodes, numDim) )
166 cleanupAndThrow(mesh, "get(Nodes_Coordinates)");
167
168 mesh->Nodes->updateTagList();
169
170 /* read elements */
171 mesh->Elements = new ElementFile((ElementTypeId)Elements_TypeId, mpi_info);
172 mesh->Elements->allocTable(num_Elements);
173 mesh->Elements->minColor = 0;
174 mesh->Elements->maxColor = num_Elements-1;
175 if (num_Elements > 0) {
176 // Elements_Id
177 if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
178 cleanupAndThrow(mesh, "get_var(Elements_Id)");
179 if (! nc_var_temp->get(&mesh->Elements->Id[0], num_Elements) )
180 cleanupAndThrow(mesh, "get(Elements_Id)");
181 // Elements_Tag
182 if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
183 cleanupAndThrow(mesh, "get_var(Elements_Tag)");
184 if (! nc_var_temp->get(&mesh->Elements->Tag[0], num_Elements) )
185 cleanupAndThrow(mesh, "get(Elements_Tag)");
186 // Elements_Owner
187 if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
188 cleanupAndThrow(mesh, "get_var(Elements_Owner)");
189 if (! nc_var_temp->get(&mesh->Elements->Owner[0], num_Elements) )
190 cleanupAndThrow(mesh, "get(Elements_Owner)");
191 // Elements_Color
192 if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
193 cleanupAndThrow(mesh, "get_var(Elements_Color)");
194 if (! nc_var_temp->get(&mesh->Elements->Color[0], num_Elements) )
195 cleanupAndThrow(mesh, "get(Elements_Color)");
196 // Elements_Nodes
197 int *Elements_Nodes = new int[num_Elements*num_Elements_numNodes];
198 if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
199 delete[] Elements_Nodes;
200 cleanupAndThrow(mesh, "get_var(Elements_Nodes)");
201 }
202 if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
203 delete[] Elements_Nodes;
204 cleanupAndThrow(mesh, "get(Elements_Nodes)");
205 }
206
207 // Copy temp array into mesh->Elements->Nodes
208 for (int i=0; i<num_Elements; i++) {
209 for (int j=0; j<num_Elements_numNodes; j++) {
210 mesh->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
211 = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
212 }
213 }
214 delete[] Elements_Nodes;
215 } // num_Elements > 0
216 mesh->Elements->updateTagList();
217
218 /* get the face elements */
219 mesh->FaceElements = new ElementFile((ElementTypeId)FaceElements_TypeId, mpi_info);
220 mesh->FaceElements->allocTable(num_FaceElements);
221 mesh->FaceElements->minColor = 0;
222 mesh->FaceElements->maxColor = num_FaceElements-1;
223 if (num_FaceElements > 0) {
224 // FaceElements_Id
225 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
226 cleanupAndThrow(mesh, "get_var(FaceElements_Id)");
227 if (! nc_var_temp->get(&mesh->FaceElements->Id[0], num_FaceElements) )
228 cleanupAndThrow(mesh, "get(FaceElements_Id)");
229 // FaceElements_Tag
230 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
231 cleanupAndThrow(mesh, "get_var(FaceElements_Tag)");
232 if (! nc_var_temp->get(&mesh->FaceElements->Tag[0], num_FaceElements) )
233 cleanupAndThrow(mesh, "get(FaceElements_Tag)");
234 // FaceElements_Owner
235 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
236 cleanupAndThrow(mesh, "get_var(FaceElements_Owner)");
237 if (! nc_var_temp->get(&mesh->FaceElements->Owner[0], num_FaceElements) )
238 cleanupAndThrow(mesh, "get(FaceElements_Owner)");
239 // FaceElements_Color
240 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
241 cleanupAndThrow(mesh, "get_var(FaceElements_Color)");
242 if (! nc_var_temp->get(&mesh->FaceElements->Color[0], num_FaceElements) )
243 cleanupAndThrow(mesh, "get(FaceElements_Color)");
244 // FaceElements_Nodes
245 int *FaceElements_Nodes = new int[num_FaceElements*num_FaceElements_numNodes];
246 if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
247 delete[] FaceElements_Nodes;
248 cleanupAndThrow(mesh, "get_var(FaceElements_Nodes)");
249 }
250 if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
251 delete[] FaceElements_Nodes;
252 cleanupAndThrow(mesh, "get(FaceElements_Nodes)");
253 }
254 // Copy temp array into mesh->FaceElements->Nodes
255 for (int i=0; i<num_FaceElements; i++) {
256 for (int j=0; j<num_FaceElements_numNodes; j++) {
257 mesh->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
258 }
259 }
260 delete[] FaceElements_Nodes;
261 } // num_FaceElements > 0
262 mesh->FaceElements->updateTagList();
263
264 // get the Points (nodal elements)
265 mesh->Points = new ElementFile((ElementTypeId)Points_TypeId, mpi_info);
266 mesh->Points->allocTable(num_Points);
267 mesh->Points->minColor = 0;
268 mesh->Points->maxColor = num_Points-1;
269 if (num_Points > 0) {
270 // Points_Id
271 if (! ( nc_var_temp = dataFile.get_var("Points_Id")))
272 cleanupAndThrow(mesh, "get_var(Points_Id)");
273 if (! nc_var_temp->get(&mesh->Points->Id[0], num_Points))
274 cleanupAndThrow(mesh, "get(Points_Id)");
275 // Points_Tag
276 if (! ( nc_var_temp = dataFile.get_var("Points_Tag")))
277 cleanupAndThrow(mesh, "get_var(Points_Tag)");
278 if (! nc_var_temp->get(&mesh->Points->Tag[0], num_Points))
279 cleanupAndThrow(mesh, "get(Points_Tag)");
280 // Points_Owner
281 if (! ( nc_var_temp = dataFile.get_var("Points_Owner")))
282 cleanupAndThrow(mesh, "get_var(Points_Owner)");
283 if (!nc_var_temp->get(&mesh->Points->Owner[0], num_Points))
284 cleanupAndThrow(mesh, "get(Points_Owner)");
285 // Points_Color
286 if (! ( nc_var_temp = dataFile.get_var("Points_Color")))
287 cleanupAndThrow(mesh, "get_var(Points_Color)");
288 if (!nc_var_temp->get(&mesh->Points->Color[0], num_Points))
289 cleanupAndThrow(mesh, "get(Points_Color)");
290 // Points_Nodes
291 int *Points_Nodes = new int[num_Points];
292 if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
293 delete[] Points_Nodes;
294 cleanupAndThrow(mesh, "get_var(Points_Nodes)");
295 }
296 if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
297 delete[] Points_Nodes;
298 cleanupAndThrow(mesh, "get(Points_Nodes)");
299 }
300 // Copy temp array into mesh->Points->Nodes
301 for (int i=0; i<num_Points; i++) {
302 mesh->Points->Id[mesh->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
303 }
304 delete[] Points_Nodes;
305 } // num_Points > 0
306 mesh->Points->updateTagList();
307
308 // get the tags
309 if (num_Tags > 0) {
310 // Temp storage to gather node IDs
311 int *Tags_keys = new int[num_Tags];
312 char name_temp[4096];
313 int i;
314
315 // Tags_keys
316 if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) ) {
317 delete[] Tags_keys;
318 cleanupAndThrow(mesh, "get_var(Tags_keys)");
319 }
320 if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
321 delete[] Tags_keys;
322 cleanupAndThrow(mesh, "get(Tags_keys)");
323 }
324 for (i=0; i<num_Tags; i++) {
325 // Retrieve tag name
326 sprintf(name_temp, "Tags_name_%d", i);
327 if (! (attr=dataFile.get_att(name_temp)) ) {
328 delete[] Tags_keys;
329 stringstream msg;
330 msg << "get_att(" << name_temp << ")";
331 cleanupAndThrow(mesh, msg.str());
332 }
333 boost::scoped_array<char> name(attr->as_string(0));
334 delete attr;
335 mesh->addTagMap(name.get(), Tags_keys[i]);
336 }
337 delete[] Tags_keys;
338 }
339
340 // Nodes_DofDistribution
341 std::vector<index_t> first_DofComponent(mpi_size+1);
342 if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) ) {
343 cleanupAndThrow(mesh, "get_var(Nodes_DofDistribution)");
344 }
345 if (! nc_var_temp->get(&first_DofComponent[0], mpi_size+1) ) {
346 cleanupAndThrow(mesh, "get(Nodes_DofDistribution)");
347 }
348
349 // Nodes_NodeDistribution
350 std::vector<index_t> first_NodeComponent(mpi_size+1);
351 if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) ) {
352 cleanupAndThrow(mesh, "get_var(Nodes_NodeDistribution)");
353 }
354 if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
355 cleanupAndThrow(mesh, "get(Nodes_NodeDistribution)");
356 }
357 mesh->createMappings(first_DofComponent, first_NodeComponent);
358
359 AbstractContinuousDomain* dom(new MeshAdapter(mesh));
360 return dom->getPtr();
361 #else
362 throw DataException("loadMesh: not compiled with NetCDF. Please contact your installation manager.");
363 #endif /* USE_NETCDF */
364 }
365
366 Domain_ptr readMesh(const std::string& fileName, int integrationOrder,
367 int reducedIntegrationOrder, bool optimize)
368 {
369 escript::JMPI mpiInfo = escript::makeInfo(MPI_COMM_WORLD);
370 Mesh* fMesh = Mesh::read(mpiInfo, fileName, optimize);
371 AbstractContinuousDomain* temp = new MeshAdapter(fMesh);
372 return temp->getPtr();
373 }
374
375 Domain_ptr readGmsh(const std::string& fileName, int numDim,
376 int integrationOrder, int reducedIntegrationOrder,
377 bool optimize)
378 {
379 escript::JMPI mpiInfo = escript::makeInfo(MPI_COMM_WORLD);
380 Mesh* fMesh = Mesh::readGmsh(mpiInfo, fileName, numDim, optimize);
381 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
382 return temp->getPtr();
383 }
384
385 Domain_ptr brick(escript::JMPI& mpi_info, double n0, double n1,double n2,
386 int order, double l0, double l1, double l2, int periodic0,
387 int periodic1, int periodic2, int integrationOrder,
388 int reducedIntegrationOrder, int useElementsOnFace,
389 int useFullElementOrder, bool optimize)
390 {
391 int numElements[]={static_cast<int>(n0),static_cast<int>(n1),static_cast<int>(n2)};
392 double length[]={l0,l1,l2};
393
394 // we don't support periodic boundary conditions
395 if (periodic0 || periodic1)
396 throw DudleyException("Dudley does not support periodic boundary conditions.");
397
398 if (integrationOrder>3 || reducedIntegrationOrder>1)
399 throw DudleyException("Dudley does not support the requested integrationOrders.");
400
401 if (useElementsOnFace || useFullElementOrder)
402 throw DudleyException("Dudley does not support useElementsOnFace or useFullElementOrder.");
403
404 if (order > 1)
405 throw DudleyException("Dudley does not support element order greater than 1.");
406
407 Mesh* fMesh = TriangularMesh_Tet4(numElements, length, optimize, mpi_info);
408
409 AbstractContinuousDomain* temp(new MeshAdapter(fMesh));
410 return temp->getPtr();
411 }
412
413 Domain_ptr brick_driver(const boost::python::list& args)
414 {
415 using boost::python::extract;
416 boost::python::object pworld=args[15];
417 escript::JMPI info;
418 if (!pworld.is_none()) {
419 extract<SubWorld_ptr> ex(pworld);
420 if (!ex.check()) {
421 throw DudleyException("Invalid escriptworld parameter.");
422 }
423 info=ex()->getMPI();
424 } else {
425 info=escript::makeInfo(MPI_COMM_WORLD);
426
427 }
428 return brick(info, static_cast<int>(extract<float>(args[0])),
429 static_cast<int>(extract<float>(args[1])),
430 static_cast<int>(extract<float>(args[2])),
431 extract<int>(args[3]), extract<double>(args[4]),
432 extract<double>(args[5]), extract<double>(args[6]),
433 extract<int>(args[7]), extract<int>(args[8]),
434 extract<int>(args[9]), extract<int>(args[10]),
435 extract<int>(args[11]), extract<int>(args[12]),
436 extract<int>(args[13]), extract<int>(args[14])
437 );
438 }
439
440
441 Domain_ptr rectangle_driver(const boost::python::list& args)
442 {
443 using boost::python::extract;
444 boost::python::object pworld=args[12];
445 escript::JMPI info;
446 if (!pworld.is_none()) {
447 extract<SubWorld_ptr> ex(pworld);
448 if (!ex.check()) {
449 throw DudleyException("Invalid escriptworld parameter.");
450 }
451 info=ex()->getMPI();
452 } else {
453 info=escript::makeInfo(MPI_COMM_WORLD);
454 }
455
456 return rectangle(info, static_cast<int>(extract<float>(args[0])),
457 static_cast<int>(extract<float>(args[1])),
458 extract<int>(args[2]), extract<double>(args[3]),
459 extract<double>(args[4]), extract<int>(args[5]),
460 extract<int>(args[6]), extract<int>(args[7]),
461 extract<int>(args[8]), extract<int>(args[9]),
462 extract<int>(args[10]), extract<int>(args[11])
463 );
464 }
465
466
467
468 Domain_ptr rectangle(escript::JMPI& mpi_info, double n0, double n1, int order,
469 double l0, double l1,
470 int periodic0,int periodic1,
471 int integrationOrder,
472 int reducedIntegrationOrder,
473 int useElementsOnFace,
474 int useFullElementOrder,
475 bool optimize)
476 {
477 int numElements[]={static_cast<int>(n0), static_cast<int>(n1)};
478 double length[]={l0,l1};
479
480 if (periodic0 || periodic1) // we don't support periodic boundary conditions
481 throw DudleyException("Dudley does not support periodic boundary conditions.");
482 if (integrationOrder>3 || reducedIntegrationOrder>1)
483 throw DudleyException("Dudley does not support the requested integrationOrders.");
484 if (useElementsOnFace || useFullElementOrder)
485 throw DudleyException("Dudley does not support useElementsOnFace or useFullElementOrder.");
486
487 if (order > 1)
488 throw DudleyException("Dudley does not support element order greater than 1.");
489 Mesh* fMesh = TriangularMesh_Tri3(numElements, length, optimize, mpi_info);
490 //
491 // Convert any dudley errors into a C++ exception
492 MeshAdapter* ma = new MeshAdapter(fMesh);
493 return Domain_ptr(ma);
494 }
495
496
497 } // namespace dudley
498

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26