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

Contents of /trunk/dudley/src/CPPAdapter/MeshAdapterFactory.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3892 - (show annotations)
Tue Apr 10 08:57:23 2012 UTC (6 years, 11 months ago) by jfenwick
File size: 24479 byte(s)
Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 tests.



1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14 #include "MeshAdapterFactory.h"
15 #include "DudleyError.h"
16 extern "C" {
17 #include "esysUtils/blocktimer.h"
18 #include "dudley/Dudley.h"
19 #include "dudley/Mesh.h"
20 #include "dudley/TriangularMesh.h"
21 #ifdef ESYS_MPI
22 #include "esysUtils/Esys_MPI.h"
23 #endif
24 }
25
26 #ifdef USE_NETCDF
27 #include <netcdfcpp.h>
28 #endif
29
30 #include <boost/python/extract.hpp>
31
32 #include <sstream>
33
34 using namespace std;
35 using namespace escript;
36
37 namespace dudley {
38
39 #ifdef USE_NETCDF
40 // A convenience method to retrieve an integer attribute from a NetCDF file
41 int NetCDF_Get_Int_Attribute(NcFile *dataFile, char *fName, char *attr_name) {
42 NcAtt *attr;
43 char error_msg[LenErrorMsg_MAX];
44 if (! (attr=dataFile->get_att(attr_name)) ) {
45 sprintf(error_msg,"loadMesh: Error retrieving integer attribute '%s' from NetCDF file '%s'", attr_name, fName);
46 throw DataException(error_msg);
47 }
48 int temp = attr->as_int(0);
49 delete attr;
50 return(temp);
51 }
52 #endif
53
54 inline void cleanupAndThrow(Dudley_Mesh* mesh, Esys_MPIInfo* info, string msg)
55 {
56 Dudley_Mesh_free(mesh);
57 Esys_MPIInfo_free(info);
58 string msgPrefix("loadMesh: NetCDF operation failed - ");
59 throw DataException(msgPrefix+msg);
60 }
61
62 // AbstractContinuousDomain* loadMesh(const std::string& fileName)
63 Domain_ptr loadMesh(const std::string& fileName)
64 {
65 #ifdef USE_NETCDF
66 Esys_MPIInfo *mpi_info = Esys_MPIInfo_alloc( MPI_COMM_WORLD );
67 Dudley_Mesh *mesh_p=NULL;
68 char error_msg[LenErrorMsg_MAX];
69
70 char *fName = Esys_MPI_appendRankToFileName(fileName.c_str(),
71 mpi_info->size,
72 mpi_info->rank);
73
74 double blocktimer_start = blocktimer_time();
75 Dudley_resetError();
76 int *first_DofComponent, *first_NodeComponent;
77
78 // Open NetCDF file for reading
79 NcAtt *attr;
80 NcVar *nc_var_temp;
81 // netCDF error handler
82 NcError err(NcError::silent_nonfatal);
83 // Create the NetCDF file.
84 NcFile dataFile(fName, NcFile::ReadOnly);
85 if (!dataFile.is_valid()) {
86 sprintf(error_msg,"loadMesh: Opening NetCDF file '%s' for reading failed.", fName);
87 Dudley_setError(IO_ERROR,error_msg);
88 Esys_MPIInfo_free( mpi_info );
89 throw DataException(error_msg);
90 }
91
92 // Read NetCDF integer attributes
93 int mpi_size = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_size");
94 int mpi_rank = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"mpi_rank");
95 int numDim = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numDim");
96 int numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"numNodes");
97 int num_Elements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements");
98 int num_FaceElements = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements");
99 int num_Points = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Points");
100 int num_Elements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Elements_numNodes");
101 int Elements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Elements_TypeId");
102 int num_FaceElements_numNodes = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_FaceElements_numNodes");
103 int FaceElements_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"FaceElements_TypeId");
104 int Points_TypeId = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"Points_TypeId");
105 int num_Tags = NetCDF_Get_Int_Attribute(&dataFile, fName, (char *)"num_Tags");
106
107 // Verify size and rank
108 if (mpi_info->size != mpi_size) {
109 sprintf(error_msg, "loadMesh: The NetCDF file '%s' can only be read on %d CPUs instead of %d", fName, mpi_size, mpi_info->size);
110 throw DataException(error_msg);
111 }
112 if (mpi_info->rank != mpi_rank) {
113 sprintf(error_msg, "loadMesh: The NetCDF file '%s' should be read on CPU #%d instead of %d", fName, mpi_rank, mpi_info->rank);
114 throw DataException(error_msg);
115 }
116
117 // Read mesh name
118 if (! (attr=dataFile.get_att("Name")) ) {
119 sprintf(error_msg,"loadMesh: Error retrieving mesh name from NetCDF file '%s'", fName);
120 throw DataException(error_msg);
121 }
122 char *name = attr->as_string(0);
123 delete attr;
124
125 TMPMEMFREE(fName);
126
127 /* allocate mesh */
128 mesh_p = Dudley_Mesh_alloc(name,numDim,mpi_info);
129 if (Dudley_noError()) {
130
131 /* read nodes */
132 Dudley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
133 // Nodes_Id
134 if (! ( nc_var_temp = dataFile.get_var("Nodes_Id")) )
135 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_Id)");
136 if (! nc_var_temp->get(&mesh_p->Nodes->Id[0], numNodes) )
137 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_Id)");
138 // Nodes_Tag
139 if (! ( nc_var_temp = dataFile.get_var("Nodes_Tag")) )
140 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_Tag)");
141 if (! nc_var_temp->get(&mesh_p->Nodes->Tag[0], numNodes) )
142 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_Tag)");
143 // Nodes_gDOF
144 if (! ( nc_var_temp = dataFile.get_var("Nodes_gDOF")) )
145 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_gDOF)");
146 if (! nc_var_temp->get(&mesh_p->Nodes->globalDegreesOfFreedom[0], numNodes) )
147 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_gDOF)");
148 // Nodes_gNI
149 if (! ( nc_var_temp = dataFile.get_var("Nodes_gNI")) )
150 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_gNI)");
151 if (! nc_var_temp->get(&mesh_p->Nodes->globalNodesIndex[0], numNodes) )
152 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_gNI)");
153 // Nodes_grDfI
154 if (! ( nc_var_temp = dataFile.get_var("Nodes_grDfI")) )
155 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_grDfI)");
156 if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedDOFIndex[0], numNodes) )
157 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_grDfI)");
158 // Nodes_grNI
159 if (! ( nc_var_temp = dataFile.get_var("Nodes_grNI")) )
160 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_grNI)");
161 if (! nc_var_temp->get(&mesh_p->Nodes->globalReducedNodesIndex[0], numNodes) )
162 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_grNI)");
163 // Nodes_Coordinates
164 if (!(nc_var_temp = dataFile.get_var("Nodes_Coordinates")))
165 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_Coordinates)");
166 if (! nc_var_temp->get(&(mesh_p->Nodes->Coordinates[0]), numNodes, numDim) )
167 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_Coordinates)");
168
169 /* read elements */
170 if (Dudley_noError()) {
171 mesh_p->Elements=Dudley_ElementFile_alloc((Dudley_ElementTypeId)Elements_TypeId, mpi_info);
172 if (Dudley_noError())
173 Dudley_ElementFile_allocTable(mesh_p->Elements, num_Elements);
174 if (Dudley_noError()) {
175 mesh_p->Elements->minColor=0;
176 mesh_p->Elements->maxColor=num_Elements-1;
177 if (num_Elements>0) {
178 // Elements_Id
179 if (! ( nc_var_temp = dataFile.get_var("Elements_Id")) )
180 cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Id)");
181 if (! nc_var_temp->get(&mesh_p->Elements->Id[0], num_Elements) )
182 cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Id)");
183 // Elements_Tag
184 if (! ( nc_var_temp = dataFile.get_var("Elements_Tag")) )
185 cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Tag)");
186 if (! nc_var_temp->get(&mesh_p->Elements->Tag[0], num_Elements) )
187 cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Tag)");
188 // Elements_Owner
189 if (! ( nc_var_temp = dataFile.get_var("Elements_Owner")) )
190 cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Owner)");
191 if (! nc_var_temp->get(&mesh_p->Elements->Owner[0], num_Elements) )
192 cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Owner)");
193 // Elements_Color
194 if (! ( nc_var_temp = dataFile.get_var("Elements_Color")) )
195 cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Color)");
196 if (! nc_var_temp->get(&mesh_p->Elements->Color[0], num_Elements) )
197 cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Color)");
198 // Elements_Nodes
199 int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
200 if (!(nc_var_temp = dataFile.get_var("Elements_Nodes"))) {
201 TMPMEMFREE(Elements_Nodes);
202 cleanupAndThrow(mesh_p, mpi_info, "get_var(Elements_Nodes)");
203 }
204 if (! nc_var_temp->get(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes) ) {
205 TMPMEMFREE(Elements_Nodes);
206 cleanupAndThrow(mesh_p, mpi_info, "get(Elements_Nodes)");
207 }
208
209 // Copy temp array into mesh_p->Elements->Nodes
210 for (int i=0; i<num_Elements; i++) {
211 for (int j=0; j<num_Elements_numNodes; j++) {
212 mesh_p->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]
213 = Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)];
214 }
215 }
216 TMPMEMFREE(Elements_Nodes);
217 } /* num_Elements>0 */
218 }
219 }
220
221 /* get the face elements */
222 if (Dudley_noError()) {
223 mesh_p->FaceElements=Dudley_ElementFile_alloc((Dudley_ElementTypeId)FaceElements_TypeId, mpi_info);
224 if (Dudley_noError())
225 Dudley_ElementFile_allocTable(mesh_p->FaceElements, num_FaceElements);
226 if (Dudley_noError()) {
227 mesh_p->FaceElements->minColor=0;
228 mesh_p->FaceElements->maxColor=num_FaceElements-1;
229 if (num_FaceElements>0) {
230 // FaceElements_Id
231 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Id")) )
232 cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Id)");
233 if (! nc_var_temp->get(&mesh_p->FaceElements->Id[0], num_FaceElements) )
234 cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Id)");
235 // FaceElements_Tag
236 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Tag")) )
237 cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Tag)");
238 if (! nc_var_temp->get(&mesh_p->FaceElements->Tag[0], num_FaceElements) )
239 cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Tag)");
240 // FaceElements_Owner
241 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Owner")) )
242 cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Owner)");
243 if (! nc_var_temp->get(&mesh_p->FaceElements->Owner[0], num_FaceElements) )
244 cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Owner)");
245 // FaceElements_Color
246 if (! ( nc_var_temp = dataFile.get_var("FaceElements_Color")) )
247 cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Color)");
248 if (! nc_var_temp->get(&mesh_p->FaceElements->Color[0], num_FaceElements) )
249 cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Color)");
250 // FaceElements_Nodes
251 int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
252 if (!(nc_var_temp = dataFile.get_var("FaceElements_Nodes"))) {
253 TMPMEMFREE(FaceElements_Nodes);
254 cleanupAndThrow(mesh_p, mpi_info, "get_var(FaceElements_Nodes)");
255 }
256 if (! nc_var_temp->get(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes) ) {
257 TMPMEMFREE(FaceElements_Nodes);
258 cleanupAndThrow(mesh_p, mpi_info, "get(FaceElements_Nodes)");
259 }
260 // Copy temp array into mesh_p->FaceElements->Nodes
261 for (int i=0; i<num_FaceElements; i++) {
262 for (int j=0; j<num_FaceElements_numNodes; j++) {
263 mesh_p->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)];
264 }
265 }
266 TMPMEMFREE(FaceElements_Nodes);
267 } /* num_FaceElements>0 */
268 }
269 }
270
271 /* get the Points (nodal elements) */
272 if (Dudley_noError()) {
273 mesh_p->Points=Dudley_ElementFile_alloc((Dudley_ElementTypeId)Points_TypeId, mpi_info);
274 if (Dudley_noError())
275 Dudley_ElementFile_allocTable(mesh_p->Points, num_Points);
276 if (Dudley_noError()) {
277 mesh_p->Points->minColor=0;
278 mesh_p->Points->maxColor=num_Points-1;
279 if (num_Points>0) {
280 // Points_Id
281 if (! ( nc_var_temp = dataFile.get_var("Points_Id")))
282 cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Id)");
283 if (! nc_var_temp->get(&mesh_p->Points->Id[0], num_Points))
284 cleanupAndThrow(mesh_p, mpi_info, "get(Points_Id)");
285 // Points_Tag
286 if (! ( nc_var_temp = dataFile.get_var("Points_Tag")))
287 cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Tag)");
288 if (! nc_var_temp->get(&mesh_p->Points->Tag[0], num_Points))
289 cleanupAndThrow(mesh_p, mpi_info, "get(Points_Tag)");
290 // Points_Owner
291 if (! ( nc_var_temp = dataFile.get_var("Points_Owner")))
292 cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Owner)");
293 if (!nc_var_temp->get(&mesh_p->Points->Owner[0], num_Points))
294 cleanupAndThrow(mesh_p, mpi_info, "get(Points_Owner)");
295 // Points_Color
296 if (! ( nc_var_temp = dataFile.get_var("Points_Color")))
297 cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Color)");
298 if (!nc_var_temp->get(&mesh_p->Points->Color[0], num_Points))
299 cleanupAndThrow(mesh_p, mpi_info, "get(Points_Color)");
300 // Points_Nodes
301 int *Points_Nodes = TMPMEMALLOC(num_Points,int);
302 if (!(nc_var_temp = dataFile.get_var("Points_Nodes"))) {
303 TMPMEMFREE(Points_Nodes);
304 cleanupAndThrow(mesh_p, mpi_info, "get_var(Points_Nodes)");
305 }
306 if (! nc_var_temp->get(&(Points_Nodes[0]), num_Points) ) {
307 TMPMEMFREE(Points_Nodes);
308 cleanupAndThrow(mesh_p, mpi_info, "get(Points_Nodes)");
309 }
310 // Copy temp array into mesh_p->Points->Nodes
311 for (int i=0; i<num_Points; i++) {
312 mesh_p->Nodes->Id[mesh_p->Points->Nodes[INDEX2(0,i,1)]] = Points_Nodes[i];
313 }
314 TMPMEMFREE(Points_Nodes);
315 } /* num_Points>0 */
316 }
317 }
318
319 /* get the tags */
320 if (Dudley_noError()) {
321 if (num_Tags>0) {
322 // Temp storage to gather node IDs
323 int *Tags_keys = TMPMEMALLOC(num_Tags, int);
324 char name_temp[4096];
325 int i;
326
327 // Tags_keys
328 if (! ( nc_var_temp = dataFile.get_var("Tags_keys")) ) {
329 TMPMEMFREE(Tags_keys);
330 cleanupAndThrow(mesh_p, mpi_info, "get_var(Tags_keys)");
331 }
332 if (! nc_var_temp->get(&Tags_keys[0], num_Tags) ) {
333 TMPMEMFREE(Tags_keys);
334 cleanupAndThrow(mesh_p, mpi_info, "get(Tags_keys)");
335 }
336 for (i=0; i<num_Tags; i++) {
337 // Retrieve tag name
338 sprintf(name_temp, "Tags_name_%d", i);
339 if (! (attr=dataFile.get_att(name_temp)) ) {
340 TMPMEMFREE(Tags_keys);
341 sprintf(error_msg,"get_att(%s)", name_temp);
342 cleanupAndThrow(mesh_p, mpi_info, error_msg);
343 }
344 char *name = attr->as_string(0);
345 delete attr;
346 Dudley_Mesh_addTagMap(mesh_p, name, Tags_keys[i]);
347 }
348 TMPMEMFREE(Tags_keys);
349 }
350 }
351
352 if (Dudley_noError()) {
353 // Nodes_DofDistribution
354 first_DofComponent = TMPMEMALLOC(mpi_size+1,index_t);
355 if (! ( nc_var_temp = dataFile.get_var("Nodes_DofDistribution")) ) {
356 TMPMEMFREE(first_DofComponent);
357 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_DofDistribution)");
358 }
359 if (! nc_var_temp->get(&first_DofComponent[0], mpi_size+1) ) {
360 TMPMEMFREE(first_DofComponent);
361 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_DofDistribution)");
362 }
363
364 // Nodes_NodeDistribution
365 first_NodeComponent = TMPMEMALLOC(mpi_size+1,index_t);
366 if (! ( nc_var_temp = dataFile.get_var("Nodes_NodeDistribution")) ) {
367 TMPMEMFREE(first_DofComponent);
368 TMPMEMFREE(first_NodeComponent);
369 cleanupAndThrow(mesh_p, mpi_info, "get_var(Nodes_NodeDistribution)");
370 }
371 if (! nc_var_temp->get(&first_NodeComponent[0], mpi_size+1) ) {
372 TMPMEMFREE(first_DofComponent);
373 TMPMEMFREE(first_NodeComponent);
374 cleanupAndThrow(mesh_p, mpi_info, "get(Nodes_NodeDistribution)");
375 }
376 Dudley_Mesh_createMappings(mesh_p, first_DofComponent, first_NodeComponent);
377 TMPMEMFREE(first_DofComponent);
378 TMPMEMFREE(first_NodeComponent);
379 }
380
381 } /* Dudley_noError() after Dudley_Mesh_alloc() */
382
383 checkDudleyError();
384 AbstractContinuousDomain* dom=new MeshAdapter(mesh_p);
385
386 if (! Dudley_noError()) {
387 Dudley_Mesh_free(mesh_p);
388 }
389
390 blocktimer_increment("LoadMesh()", blocktimer_start);
391 return dom->getPtr();
392 #else
393 throw DataException("loadMesh: not compiled with NetCDF. Please contact your installation manager.");
394 #endif /* USE_NETCDF */
395 }
396
397 Domain_ptr readMesh(const std::string& fileName,
398 int integrationOrder,
399 int reducedIntegrationOrder,
400 int optimize)
401 {
402 //
403 // create a copy of the filename to overcome the non-constness of call
404 // to Dudley_Mesh_read
405 Dudley_Mesh* fMesh=0;
406 // Win32 refactor
407 if( fileName.size() == 0 )
408 {
409 throw DataException("Null file name!");
410 }
411
412 char *fName = TMPMEMALLOC(fileName.size()+1,char);
413
414 strcpy(fName,fileName.c_str());
415 double blocktimer_start = blocktimer_time();
416
417 fMesh=Dudley_Mesh_read(fName,integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
418 checkDudleyError();
419 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
420
421 /* win32 refactor */
422 TMPMEMFREE(fName);
423
424 blocktimer_increment("ReadMesh()", blocktimer_start);
425 return temp->getPtr();
426 }
427
428 Domain_ptr readGmsh(const std::string& fileName,
429 int numDim,
430 int integrationOrder,
431 int reducedIntegrationOrder,
432 int optimize,
433 int useMacroElements)
434 {
435 //
436 // create a copy of the filename to overcome the non-constness of call
437 // to Dudley_Mesh_read
438 Dudley_Mesh* fMesh=0;
439 // Win32 refactor
440 if( fileName.size() == 0 )
441 {
442 throw DataException("Null file name!");
443 }
444
445 char *fName = TMPMEMALLOC(fileName.size()+1,char);
446
447 strcpy(fName,fileName.c_str());
448 double blocktimer_start = blocktimer_time();
449
450 fMesh=Dudley_Mesh_readGmsh(fName, numDim, integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE), (useMacroElements ? TRUE : FALSE));
451 checkDudleyError();
452 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
453
454 /* win32 refactor */
455 TMPMEMFREE(fName);
456
457 blocktimer_increment("ReadGmsh()", blocktimer_start);
458 return temp->getPtr();
459 }
460
461 Domain_ptr brick(double n0, double n1,double n2,int order,
462 double l0,double l1,double l2,
463 int periodic0,int periodic1,
464 int periodic2,
465 int integrationOrder,
466 int reducedIntegrationOrder,
467 int useElementsOnFace,
468 int useFullElementOrder,
469 int optimize)
470 {
471 int numElements[]={static_cast<int>(n0),static_cast<int>(n1),static_cast<int>(n2)};
472 double length[]={l0,l1,l2};
473
474 if (periodic0 || periodic1) // we don't support periodic boundary conditions
475 {
476 throw DudleyAdapterException("Dudley does not support periodic boundary conditions.");
477 }
478 else if (integrationOrder>3 || reducedIntegrationOrder>1)
479 {
480 throw DudleyAdapterException("Dudley does not support the requested integrationOrders.");
481 }
482 else if (useElementsOnFace || useFullElementOrder)
483 {
484 throw DudleyAdapterException("Dudley does not support useElementsOnFace or useFullElementOrder.");
485 }
486 if (order>1)
487 {
488 throw DudleyAdapterException("Dudley does not support element order greater than 1.");
489 }
490
491 //
492 // linearInterpolation
493 Dudley_Mesh* fMesh=NULL;
494
495 fMesh=Dudley_TriangularMesh_Tet4(numElements, length, integrationOrder,
496 reducedIntegrationOrder, (optimize ? TRUE : FALSE));
497
498 //
499 // Convert any dudley errors into a C++ exception
500 checkDudleyError();
501 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
502 return temp->getPtr();
503 }
504
505 Domain_ptr rectangle(double n0, double n1, int order,
506 double l0, double l1,
507 int periodic0,int periodic1,
508 int integrationOrder,
509 int reducedIntegrationOrder,
510 int useElementsOnFace,
511 int useFullElementOrder,
512 int optimize)
513 {
514 int numElements[]={static_cast<int>(n0), static_cast<int>(n1)};
515 double length[]={l0,l1};
516
517 if (periodic0 || periodic1) // we don't support periodic boundary conditions
518 {
519 throw DudleyAdapterException("Dudley does not support periodic boundary conditions.");
520 }
521 else if (integrationOrder>3 || reducedIntegrationOrder>1)
522 {
523 throw DudleyAdapterException("Dudley does not support the requested integrationOrders.");
524 }
525 else if (useElementsOnFace || useFullElementOrder)
526 {
527 throw DudleyAdapterException("Dudley does not support useElementsOnFace or useFullElementOrder.");
528 }
529
530 if (order>1)
531 {
532 throw DudleyAdapterException("Dudley does not support element order greater than 1.");
533 }
534 Dudley_Mesh* fMesh=Dudley_TriangularMesh_Tri3(numElements, length,
535 integrationOrder, reducedIntegrationOrder, (optimize ? TRUE : FALSE));
536 //
537 // Convert any dudley errors into a C++ exception
538 checkDudleyError();
539 AbstractContinuousDomain* temp=new MeshAdapter(fMesh);
540
541 return temp->getPtr();
542 }
543
544 // end of namespace
545
546 }
547

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26