/[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 6009 - (show annotations)
Wed Mar 2 04:13:26 2016 UTC (3 years, 1 month ago) by caltinay
File size: 22884 byte(s)
Much needed sync with trunk...

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(Dudley_Mesh* mesh, string msg)
57 {
58 Dudley_Mesh_free(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 Dudley_Mesh *mesh_p=NULL;
68 const string fName(mpi_info->appendRankToFileName(fileName));
69
70 int *first_DofComponent, *first_NodeComponent;
71
72 // Open NetCDF file for reading
73 NcAtt *attr;
74 NcVar *nc_var_temp;
75 // netCDF error handler
76 NcError err(NcError::silent_nonfatal);
77 // Create the NetCDF file.
78 NcFile dataFile(fName.c_str(), NcFile::ReadOnly);
79 if (!dataFile.is_valid()) {
80 stringstream msg;
81 msg << "loadMesh: Opening NetCDF file '" << fName << "' for reading failed.";
82 throw escript::IOError(msg.str());
83 }
84
85 // Read NetCDF integer attributes
86
87 // index_size was only introduced with 64-bit index support so fall back
88 // to 32 bits if not found.
89 int index_size;
90 try {
91 index_size = ncReadAtt<int>(&dataFile, fName, "index_size");
92 } catch (escript::IOError& e) {
93 index_size = 4;
94 }
95 // technically we could cast if reading 32-bit data on 64-bit escript
96 // but cost-benefit analysis clearly favours this implementation for now
97 if (sizeof(index_t) != index_size) {
98 throw escript::IOError("loadMesh: size of index types at runtime differ from dump file");
99 }
100
101 int mpi_size = ncReadAtt<int>(&dataFile, fName, "mpi_size");
102 int mpi_rank = ncReadAtt<int>(&dataFile, fName, "mpi_rank");
103 int numDim = ncReadAtt<int>(&dataFile, fName, "numDim");
104 dim_t numNodes = ncReadAtt<dim_t>(&dataFile, fName, "numNodes");
105 dim_t num_Elements = ncReadAtt<dim_t>(&dataFile, fName, "num_Elements");
106 dim_t num_FaceElements = ncReadAtt<dim_t>(&dataFile, fName, "num_FaceElements");
107 dim_t num_Points = ncReadAtt<dim_t>(&dataFile, fName, "num_Points");
108 int num_Elements_numNodes = ncReadAtt<int>(&dataFile, fName, "num_Elements_numNodes");
109 int Elements_TypeId = ncReadAtt<int>(&dataFile, fName, "Elements_TypeId");
110 int num_FaceElements_numNodes = ncReadAtt<int>(&dataFile, fName, "num_FaceElements_numNodes");
111 int FaceElements_TypeId = ncReadAtt<int>(&dataFile, fName, "FaceElements_TypeId");
112 int Points_TypeId = ncReadAtt<int>(&dataFile, fName, "Points_TypeId");
113 int num_Tags = ncReadAtt<int>(&dataFile, fName, "num_Tags");
114
115 // Verify size and rank
116 if (mpi_info->size != mpi_size) {
117 stringstream msg;
118 msg << "loadMesh: The NetCDF file '" << fName
119 << "' can only be read on " << mpi_size
120 << " CPUs. Currently running: " << mpi_info->size;
121 throw DudleyException(msg.str());
122 }
123 if (mpi_info->rank != mpi_rank) {
124 stringstream msg;
125 msg << "loadMesh: The NetCDF file '" << fName
126 << "' should be read on CPU #" << mpi_rank
127 << " and NOT on #" << mpi_info->rank;
128 throw DudleyException(msg.str());
129 }
130
131 // Read mesh name
132 if (! (attr=dataFile.get_att("Name")) ) {
133 stringstream msg;
134 msg << "loadMesh: Error retrieving mesh name from NetCDF file '"
135 << fName << "'";
136 throw escript::IOError(msg.str());
137 }
138 boost::scoped_array<char> name(attr->as_string(0));
139 delete attr;
140
141 // allocate mesh
142 mesh_p = Dudley_Mesh_alloc(name.get(), numDim, mpi_info);
143 // read nodes
144 Dudley_NodeFile_allocTable(mesh_p->Nodes, 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
181 Dudley_NodeFile_setTagsInUse(mesh_p->Nodes);
182
183 /* read elements */
184 mesh_p->Elements=Dudley_ElementFile_alloc((Dudley_ElementTypeId)Elements_TypeId, mpi_info);
185 Dudley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
186 mesh_p->Elements->minColor=0;
187 mesh_p->Elements->maxColor=num_Elements-1;
188 if (num_Elements>0) {
189 // Elements_Id
190 if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
191 cleanupAndThrow(mesh_p, "get_var(Elements_Id)");
192 if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) )
193 cleanupAndThrow(mesh_p, "get(Elements_Id)");
194 // Elements_Tag
195 if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
196 cleanupAndThrow(mesh_p, "get_var(Elements_Tag)");
197 if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) )
198 cleanupAndThrow(mesh_p, "get(Elements_Tag)");
199 // Elements_Owner
200 if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
201 cleanupAndThrow(mesh_p, "get_var(Elements_Owner)");
202 if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) )
203 cleanupAndThrow(mesh_p, "get(Elements_Owner)");
204 // Elements_Color
205 if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
206 cleanupAndThrow(mesh_p, "get_var(Elements_Color)");
207 if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) )
208 cleanupAndThrow(mesh_p, "get(Elements_Color)");
209 // Elements_Nodes
210 int *Elements_Nodes = new int[num_Elements*num_Elements_numNodes];
211 if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
212 delete[] Elements_Nodes;
213 cleanupAndThrow(mesh_p, "get_var(Elements_Nodes)");
214 }
215 if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
216 delete[] Elements_Nodes;
217 cleanupAndThrow(mesh_p, "get(Elements_Nodes)");
218 }
219
220 // Copy temp array into mesh_p->Elements->Nodes
221 for (int i=0; i<num_Elements; i++) {
222 for (int j=0; j<num_Elements_numNodes; j++) {
223 mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
224 = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
225 }
226 }
227 delete[] Elements_Nodes;
228 } /* num_Elements>0 */
229 Dudley_ElementFile_setTagsInUse(mesh_p->Elements);
230
231 /* get the face elements */
232 mesh_p->FaceElements=Dudley_ElementFile_alloc((Dudley_ElementTypeId)FaceElements_TypeId, mpi_info);
233 Dudley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
234 mesh_p->FaceElements->minColor=0;
235 mesh_p->FaceElements->maxColor=num_FaceElements-1;
236 if (num_FaceElements>0) {
237 // FaceElements_Id
238 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
239 cleanupAndThrow(mesh_p, "get_var(FaceElements_Id)");
240 if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) )
241 cleanupAndThrow(mesh_p, "get(FaceElements_Id)");
242 // FaceElements_Tag
243 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
244 cleanupAndThrow(mesh_p, "get_var(FaceElements_Tag)");
245 if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) )
246 cleanupAndThrow(mesh_p, "get(FaceElements_Tag)");
247 // FaceElements_Owner
248 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
249 cleanupAndThrow(mesh_p, "get_var(FaceElements_Owner)");
250 if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) )
251 cleanupAndThrow(mesh_p, "get(FaceElements_Owner)");
252 // FaceElements_Color
253 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
254 cleanupAndThrow(mesh_p, "get_var(FaceElements_Color)");
255 if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) )
256 cleanupAndThrow(mesh_p, "get(FaceElements_Color)");
257 // FaceElements_Nodes
258 int *FaceElements_Nodes = new int[num_FaceElements*num_FaceElements_numNodes];
259 if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
260 delete[] FaceElements_Nodes;
261 cleanupAndThrow(mesh_p, "get_var(FaceElements_Nodes)");
262 }
263 if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
264 delete[] FaceElements_Nodes;
265 cleanupAndThrow(mesh_p, "get(FaceElements_Nodes)");
266 }
267 // Copy temp array into mesh_p->FaceElements->Nodes
268 for (int i=0; i<num_FaceElements; i++) {
269 for (int j=0; j<num_FaceElements_numNodes; j++) {
270 mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
271 }
272 }
273 delete[] FaceElements_Nodes;
274 } /* num_FaceElements>0 */
275 Dudley_ElementFile_setTagsInUse(mesh_p->FaceElements);
276
277 /* get the Points (nodal elements) */
278 mesh_p->Points=Dudley_ElementFile_alloc((Dudley_ElementTypeId)Points_TypeId, mpi_info);
279 Dudley_ElementFile_allocTable(mesh_p->Points, num_Points);
280 mesh_p->Points->minColor=0;
281 mesh_p->Points->maxColor=num_Points-1;
282 if (num_Points>0) {
283 // Points_Id
284 if (! ( nc_var_temp = dataFile.get_var("Points_Id")))
285 cleanupAndThrow(mesh_p, "get_var(Points_Id)");
286 if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points))
287 cleanupAndThrow(mesh_p, "get(Points_Id)");
288 // Points_Tag
289 if (! ( nc_var_temp = dataFile.get_var("Points_Tag")))
290 cleanupAndThrow(mesh_p, "get_var(Points_Tag)");
291 if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points))
292 cleanupAndThrow(mesh_p, "get(Points_Tag)");
293 // Points_Owner
294 if (! ( nc_var_temp = dataFile.get_var("Points_Owner")))
295 cleanupAndThrow(mesh_p, "get_var(Points_Owner)");
296 if (!nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points))
297 cleanupAndThrow(mesh_p, "get(Points_Owner)");
298 // Points_Color
299 if (! ( nc_var_temp = dataFile.get_var("Points_Color")))
300 cleanupAndThrow(mesh_p, "get_var(Points_Color)");
301 if (!nc_var_temp->get(&mesh_p->Points->Color[0], num_Points))
302 cleanupAndThrow(mesh_p, "get(Points_Color)");
303 // Points_Nodes
304 int *Points_Nodes = new int[num_Points];
305 if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
306 delete[] Points_Nodes;
307 cleanupAndThrow(mesh_p, "get_var(Points_Nodes)");
308 }
309 if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
310 delete[] Points_Nodes;
311 cleanupAndThrow(mesh_p, "get(Points_Nodes)");
312 }
313 // Copy temp array into mesh_p->Points->Nodes
314 for (int i=0; i<num_Points; i++) {
315 mesh_p->Points->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
316 }
317 delete[] Points_Nodes;
318 } /* num_Points>0 */
319 Dudley_ElementFile_setTagsInUse(mesh_p->Points);
320
321 /* get the tags */
322 if (num_Tags>0) {
323 // Temp storage to gather node IDs
324 int *Tags_keys = new int[num_Tags];
325 char name_temp[4096];
326 int i;
327
328 // Tags_keys
329 if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) ) {
330 delete[] Tags_keys;
331 cleanupAndThrow(mesh_p, "get_var(Tags_keys)");
332 }
333 if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
334 delete[] Tags_keys;
335 cleanupAndThrow(mesh_p, "get(Tags_keys)");
336 }
337 for (i=0; i<num_Tags; i++) {
338 // Retrieve tag name
339 sprintf(name_temp, "Tags_name_%d", i);
340 if (! (attr=dataFile.get_att(name_temp)) ) {
341 delete[] Tags_keys;
342 stringstream msg;
343 msg << "get_att(" << name_temp << ")";
344 cleanupAndThrow(mesh_p, msg.str());
345 }
346 boost::scoped_array<char> name(attr->as_string(0));
347 delete attr;
348 Dudley_Mesh_addTagMap(mesh_p, name.get(), Tags_keys[i]);
349 }
350 delete[] Tags_keys;
351 }
352
353 // Nodes_DofDistribution
354 first_DofComponent = new index_t[mpi_size+1];
355 if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) ) {
356 delete[] first_DofComponent;
357 cleanupAndThrow(mesh_p, "get_var(Nodes_DofDistribution)");
358 }
359 if (! nc_var_temp->get(&first_DofComponent[0], mpi_size+1) ) {
360 delete[] first_DofComponent;
361 cleanupAndThrow(mesh_p, "get(Nodes_DofDistribution)");
362 }
363
364 // Nodes_NodeDistribution
365 first_NodeComponent = new index_t[mpi_size+1];
366 if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) ) {
367 delete[] first_DofComponent;
368 delete[] first_NodeComponent;
369 cleanupAndThrow(mesh_p, "get_var(Nodes_NodeDistribution)");
370 }
371 if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
372 delete[] first_DofComponent;
373 delete[] first_NodeComponent;
374 cleanupAndThrow(mesh_p, "get(Nodes_NodeDistribution)");
375 }
376 Dudley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
377 delete[] first_DofComponent;
378 delete[] first_NodeComponent;
379
380 AbstractContinuousDomain* dom=new MeshAdapter(mesh_p);
381 return dom->getPtr();
382 #else
383 throw DataException("loadMesh: not compiled with NetCDF. Please contact your installation manager.");
384 #endif /* USE_NETCDF */
385 }
386
387 Domain_ptr readMesh(const std::string& fileName,
388 int integrationOrder,
389 int reducedIntegrationOrder,
390 int optimize)
391 {
392 //
393 // create a copy of the filename to overcome the non-constness of call
394 // to Dudley_Mesh_read
395 Dudley_Mesh* fMesh=0;
396 // Win32 refactor
397 if( fileName.size() == 0 )
398 {
399 throw DataException("Null file name!");
400 }
401
402 char *fName = new char[fileName.size()+1];
403
404 strcpy(fName,fileName.c_str());
405
406 fMesh=Dudley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, optimize!=0);
407 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
408
409 delete[] fName;
410
411 return temp->getPtr();
412 }
413
414 Domain_ptr readGmsh(const std::string& fileName,
415 int numDim,
416 int integrationOrder,
417 int reducedIntegrationOrder,
418 int optimize)
419 {
420 //
421 // create a copy of the filename to overcome the non-constness of call
422 // to Dudley_Mesh_read
423 Dudley_Mesh* fMesh=0;
424 if (fileName.size() == 0)
425 throw DudleyException("Null file name!");
426
427 char *fName = new char[fileName.size()+1];
428 strcpy(fName,fileName.c_str());
429 fMesh=Dudley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, optimize!=0);
430 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
431 delete[] fName;
432 return temp->getPtr();
433 }
434
435 Domain_ptr brick(escript::JMPI& mpi_info, double n0, double n1,double n2,int order,
436 double l0,double l1,double l2,
437 int periodic0,int periodic1,
438 int periodic2,
439 int integrationOrder,
440 int reducedIntegrationOrder,
441 int useElementsOnFace,
442 int useFullElementOrder,
443 int optimize)
444 {
445 int numElements[]={static_cast<int>(n0),static_cast<int>(n1),static_cast<int>(n2)};
446 double length[]={l0,l1,l2};
447
448 if (periodic0 || periodic1) // we don't support periodic boundary conditions
449 {
450 throw DudleyException("Dudley does not support periodic boundary conditions.");
451 }
452 else if (integrationOrder>3 || reducedIntegrationOrder>1)
453 {
454 throw DudleyException("Dudley does not support the requested integrationOrders.");
455 }
456 else if (useElementsOnFace || useFullElementOrder)
457 {
458 throw DudleyException("Dudley does not support useElementsOnFace or useFullElementOrder.");
459 }
460 if (order>1)
461 {
462 throw DudleyException("Dudley does not support element order greater than 1.");
463 }
464
465 //
466 // linearInterpolation
467 Dudley_Mesh* fMesh=NULL;
468
469 fMesh=Dudley_TriangularMesh_Tet4(numElements, length, integrationOrder,
470 reducedIntegrationOrder, optimize!=0,
471 mpi_info);
472
473 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
474 return temp->getPtr();
475 }
476
477 Domain_ptr brick_driver(const boost::python::list& args)
478 {
479 using boost::python::extract;
480 boost::python::object pworld=args[15];
481 escript::JMPI info;
482 if (!pworld.is_none()) {
483 extract<SubWorld_ptr> ex(pworld);
484 if (!ex.check()) {
485 throw DudleyException("Invalid escriptworld parameter.");
486 }
487 info=ex()->getMPI();
488 } else {
489 info=escript::makeInfo(MPI_COMM_WORLD);
490
491 }
492 return brick(info, static_cast<int>(extract<float>(args[0])),
493 static_cast<int>(extract<float>(args[1])),
494 static_cast<int>(extract<float>(args[2])),
495 extract<int>(args[3]), extract<double>(args[4]),
496 extract<double>(args[5]), extract<double>(args[6]),
497 extract<int>(args[7]), extract<int>(args[8]),
498 extract<int>(args[9]), extract<int>(args[10]),
499 extract<int>(args[11]), extract<int>(args[12]),
500 extract<int>(args[13]), extract<int>(args[14])
501 );
502 }
503
504
505 Domain_ptr rectangle_driver(const boost::python::list& args)
506 {
507 using boost::python::extract;
508 boost::python::object pworld=args[12];
509 escript::JMPI info;
510 if (!pworld.is_none()) {
511 extract<SubWorld_ptr> ex(pworld);
512 if (!ex.check()) {
513 throw DudleyException("Invalid escriptworld parameter.");
514 }
515 info=ex()->getMPI();
516 } else {
517 info=escript::makeInfo(MPI_COMM_WORLD);
518 }
519
520 return rectangle(info, static_cast<int>(extract<float>(args[0])),
521 static_cast<int>(extract<float>(args[1])),
522 extract<int>(args[2]), extract<double>(args[3]),
523 extract<double>(args[4]), extract<int>(args[5]),
524 extract<int>(args[6]), extract<int>(args[7]),
525 extract<int>(args[8]), extract<int>(args[9]),
526 extract<int>(args[10]), extract<int>(args[11])
527 );
528 }
529
530
531
532 Domain_ptr rectangle(escript::JMPI& mpi_info, double n0, double n1, int order,
533 double l0, double l1,
534 int periodic0,int periodic1,
535 int integrationOrder,
536 int reducedIntegrationOrder,
537 int useElementsOnFace,
538 int useFullElementOrder,
539 int optimize)
540 {
541 int numElements[]={static_cast<int>(n0), static_cast<int>(n1)};
542 double length[]={l0,l1};
543
544 if (periodic0 || periodic1) // we don't support periodic boundary conditions
545 {
546 throw DudleyException("Dudley does not support periodic boundary conditions.");
547 }
548 else if (integrationOrder>3 || reducedIntegrationOrder>1)
549 {
550 throw DudleyException("Dudley does not support the requested integrationOrders.");
551 }
552 else if (useElementsOnFace || useFullElementOrder)
553 {
554 throw DudleyException("Dudley does not support useElementsOnFace or useFullElementOrder.");
555 }
556
557 if (order>1)
558 {
559 throw DudleyException("Dudley does not support element order greater than 1.");
560 }
561 Dudley_Mesh* fMesh=Dudley_TriangularMesh_Tri3(numElements, length,
562 integrationOrder, reducedIntegrationOrder, optimize!=0,
563 mpi_info);
564 //
565 // Convert any dudley errors into a C++ exception
566 MeshAdapter* ma=new MeshAdapter(fMesh);
567 return Domain_ptr(ma);
568 }
569
570
571 } // namespace dudley
572

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26