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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3317 - (show annotations)
Thu Oct 28 00:50:41 2010 UTC (8 years, 9 months ago) by caltinay
File size: 72706 byte(s)
Removed bogus mpi.h includes, replaced others by our esysUtils wrapper
and rearranged so that the wrapper always comes before netcdf which fixes
linking problems when disabling mpi on systems with netCDF 4.x.

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
15 #include "MeshAdapter.h"
16 #include "escript/Data.h"
17 #include "escript/DataFactory.h"
18 #ifdef USE_NETCDF
19 #include <netcdfcpp.h>
20 #endif
21 extern "C" {
22 #include "esysUtils/blocktimer.h"
23 }
24
25 using namespace std;
26 using namespace escript;
27
28 namespace dudley {
29
30 //
31 // define the static constants
32 MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
33 const int MeshAdapter::DegreesOfFreedom=DUDLEY_DEGREES_OF_FREEDOM;
34 const int MeshAdapter::ReducedDegreesOfFreedom=DUDLEY_REDUCED_DEGREES_OF_FREEDOM;
35 const int MeshAdapter::Nodes=DUDLEY_NODES;
36 const int MeshAdapter::ReducedNodes=DUDLEY_REDUCED_NODES;
37 const int MeshAdapter::Elements=DUDLEY_ELEMENTS;
38 const int MeshAdapter::ReducedElements=DUDLEY_REDUCED_ELEMENTS;
39 const int MeshAdapter::FaceElements=DUDLEY_FACE_ELEMENTS;
40 const int MeshAdapter::ReducedFaceElements=DUDLEY_REDUCED_FACE_ELEMENTS;
41 const int MeshAdapter::Points=DUDLEY_POINTS;
42
43 MeshAdapter::MeshAdapter(Dudley_Mesh* dudleyMesh)
44 {
45 setFunctionSpaceTypeNames();
46 //
47 // need to use a null_deleter as Dudley_Mesh_free deletes the pointer
48 // for us.
49 m_dudleyMesh.reset(dudleyMesh,null_deleter());
50 }
51
52 //
53 // The copy constructor should just increment the use count
54 MeshAdapter::MeshAdapter(const MeshAdapter& in):
55 m_dudleyMesh(in.m_dudleyMesh)
56 {
57 setFunctionSpaceTypeNames();
58 }
59
60 MeshAdapter::~MeshAdapter()
61 {
62 //
63 // I hope the case for the pointer being zero has been taken care of.
64 // cout << "In MeshAdapter destructor." << endl;
65 if (m_dudleyMesh.unique()) {
66 Dudley_Mesh_free(m_dudleyMesh.get());
67 }
68 }
69
70 int MeshAdapter::getMPISize() const
71 {
72 return m_dudleyMesh.get()->MPIInfo->size;
73 }
74 int MeshAdapter::getMPIRank() const
75 {
76 return m_dudleyMesh.get()->MPIInfo->rank;
77 }
78 void MeshAdapter::MPIBarrier() const
79 {
80 #ifdef ESYS_MPI
81 MPI_Barrier(m_dudleyMesh.get()->MPIInfo->comm);
82 #endif
83 return;
84 }
85 bool MeshAdapter::onMasterProcessor() const
86 {
87 return m_dudleyMesh.get()->MPIInfo->rank == 0;
88 }
89
90
91 #ifdef ESYS_MPI
92 MPI_Comm
93 #else
94 unsigned int
95 #endif
96 MeshAdapter::getMPIComm() const
97 {
98 #ifdef ESYS_MPI
99 return m_dudleyMesh->MPIInfo->comm;
100 #else
101 return 0;
102 #endif
103 }
104
105
106 Dudley_Mesh* MeshAdapter::getDudley_Mesh() const {
107 return m_dudleyMesh.get();
108 }
109
110 void MeshAdapter::write(const string& fileName) const
111 {
112 char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
113 strcpy(fName,fileName.c_str());
114 Dudley_Mesh_write(m_dudleyMesh.get(),fName);
115 checkDudleyError();
116 TMPMEMFREE(fName);
117 }
118
119 void MeshAdapter::Print_Mesh_Info(const bool full) const
120 {
121 Dudley_PrintMesh_Info(m_dudleyMesh.get(), full);
122 }
123
124 void MeshAdapter::dump(const string& fileName) const
125 {
126 #ifdef USE_NETCDF
127 const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
128 NcVar *ids;
129 int *int_ptr;
130 Dudley_Mesh *mesh = m_dudleyMesh.get();
131 Dudley_TagMap* tag_map;
132 int num_Tags = 0;
133 int mpi_size = mesh->MPIInfo->size;
134 int mpi_rank = mesh->MPIInfo->rank;
135 int numDim = mesh->Nodes->numDim;
136 int numNodes = mesh->Nodes->numNodes;
137 int num_Elements = mesh->Elements->numElements;
138 int num_FaceElements = mesh->FaceElements->numElements;
139 int num_Points = mesh->Points->numElements;
140 int num_Elements_numNodes = mesh->Elements->numNodes;
141 int num_FaceElements_numNodes = mesh->FaceElements->numNodes;
142 #ifdef ESYS_MPI
143 MPI_Status status;
144 #endif
145
146 /* Incoming token indicates it's my turn to write */
147 #ifdef ESYS_MPI
148 if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
149 #endif
150
151 char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
152 mpi_size, mpi_rank);
153
154 /* Figure out how much storage is required for tags */
155 tag_map = mesh->TagMap;
156 num_Tags = 0;
157 while (tag_map) {
158 num_Tags++;
159 tag_map=tag_map->next;
160 }
161
162 // NetCDF error handler
163 NcError err(NcError::verbose_nonfatal);
164 // Create the file.
165 NcFile dataFile(newFileName, NcFile::Replace);
166 string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
167 // check if writing was successful
168 if (!dataFile.is_valid())
169 throw DataException(msgPrefix+"Open file for output");
170
171 // Define dimensions (num_Elements and dim_Elements are identical,
172 // dim_Elements only appears if > 0)
173 if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
174 throw DataException(msgPrefix+"add_dim(numNodes)");
175 if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
176 throw DataException(msgPrefix+"add_dim(numDim)");
177 if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )
178 throw DataException(msgPrefix+"add_dim(mpi_size)");
179 if (num_Elements>0)
180 if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
181 throw DataException(msgPrefix+"add_dim(dim_Elements)");
182 if (num_FaceElements>0)
183 if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
184 throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
185 if (num_Points>0)
186 if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
187 throw DataException(msgPrefix+"add_dim(dim_Points)");
188 if (num_Elements>0)
189 if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
190 throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
191 if (num_FaceElements>0)
192 if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
193 throw DataException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
194 if (num_Tags>0)
195 if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
196 throw DataException(msgPrefix+"add_dim(dim_Tags)");
197
198 // Attributes: MPI size, MPI rank, Name, order, reduced_order
199 if (!dataFile.add_att("mpi_size", mpi_size) )
200 throw DataException(msgPrefix+"add_att(mpi_size)");
201 if (!dataFile.add_att("mpi_rank", mpi_rank) )
202 throw DataException(msgPrefix+"add_att(mpi_rank)");
203 if (!dataFile.add_att("Name",mesh->Name) )
204 throw DataException(msgPrefix+"add_att(Name)");
205 if (!dataFile.add_att("numDim",numDim) )
206 throw DataException(msgPrefix+"add_att(order)");
207 if (!dataFile.add_att("order",mesh->integrationOrder) )
208 throw DataException(msgPrefix+"add_att(order)");
209 if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
210 throw DataException(msgPrefix+"add_att(reduced_order)");
211 if (!dataFile.add_att("numNodes",numNodes) )
212 throw DataException(msgPrefix+"add_att(numNodes)");
213 if (!dataFile.add_att("num_Elements",num_Elements) )
214 throw DataException(msgPrefix+"add_att(num_Elements)");
215 if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
216 throw DataException(msgPrefix+"add_att(num_FaceElements)");
217 if (!dataFile.add_att("num_Points",num_Points) )
218 throw DataException(msgPrefix+"add_att(num_Points)");
219 if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
220 throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
221 if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
222 throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
223 if (!dataFile.add_att("Elements_TypeId", mesh->Elements->etype) )
224 throw DataException(msgPrefix+"add_att(Elements_TypeId)");
225 if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->etype) )
226 throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
227 if (!dataFile.add_att("Points_TypeId", mesh->Points->etype) )
228 throw DataException(msgPrefix+"add_att(Points_TypeId)");
229 if (!dataFile.add_att("num_Tags", num_Tags) )
230 throw DataException(msgPrefix+"add_att(num_Tags)");
231
232 // // // // // Nodes // // // // //
233
234 // Nodes nodeDistribution
235 if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
236 throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
237 int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
238 if (! (ids->put(int_ptr, mpi_size+1)) )
239 throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
240
241 // Nodes degreesOfFreedomDistribution
242 if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
243 throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
244 int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
245 if (! (ids->put(int_ptr, mpi_size+1)) )
246 throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
247
248 // Only write nodes if non-empty because NetCDF doesn't like empty arrays
249 // (it treats them as NC_UNLIMITED)
250 if (numNodes>0) {
251
252 // Nodes Id
253 if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
254 throw DataException(msgPrefix+"add_var(Nodes_Id)");
255 int_ptr = &mesh->Nodes->Id[0];
256 if (! (ids->put(int_ptr, numNodes)) )
257 throw DataException(msgPrefix+"put(Nodes_Id)");
258
259 // Nodes Tag
260 if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
261 throw DataException(msgPrefix+"add_var(Nodes_Tag)");
262 int_ptr = &mesh->Nodes->Tag[0];
263 if (! (ids->put(int_ptr, numNodes)) )
264 throw DataException(msgPrefix+"put(Nodes_Tag)");
265
266 // Nodes gDOF
267 if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
268 throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
269 int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
270 if (! (ids->put(int_ptr, numNodes)) )
271 throw DataException(msgPrefix+"put(Nodes_gDOF)");
272
273 // Nodes global node index
274 if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
275 throw DataException(msgPrefix+"add_var(Nodes_gNI)");
276 int_ptr = &mesh->Nodes->globalNodesIndex[0];
277 if (! (ids->put(int_ptr, numNodes)) )
278 throw DataException(msgPrefix+"put(Nodes_gNI)");
279
280 // Nodes grDof
281 if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
282 throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
283 int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
284 if (! (ids->put(int_ptr, numNodes)) )
285 throw DataException(msgPrefix+"put(Nodes_grDfI)");
286
287 // Nodes grNI
288 if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
289 throw DataException(msgPrefix+"add_var(Nodes_grNI)");
290 int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
291 if (! (ids->put(int_ptr, numNodes)) )
292 throw DataException(msgPrefix+"put(Nodes_grNI)");
293
294 // Nodes Coordinates
295 if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
296 throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
297 if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
298 throw DataException(msgPrefix+"put(Nodes_Coordinates)");
299
300 }
301
302 // // // // // Elements // // // // //
303
304 if (num_Elements>0) {
305
306 // Elements_Id
307 if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
308 throw DataException(msgPrefix+"add_var(Elements_Id)");
309 int_ptr = &mesh->Elements->Id[0];
310 if (! (ids->put(int_ptr, num_Elements)) )
311 throw DataException(msgPrefix+"put(Elements_Id)");
312
313 // Elements_Tag
314 if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
315 throw DataException(msgPrefix+"add_var(Elements_Tag)");
316 int_ptr = &mesh->Elements->Tag[0];
317 if (! (ids->put(int_ptr, num_Elements)) )
318 throw DataException(msgPrefix+"put(Elements_Tag)");
319
320 // Elements_Owner
321 if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
322 throw DataException(msgPrefix+"add_var(Elements_Owner)");
323 int_ptr = &mesh->Elements->Owner[0];
324 if (! (ids->put(int_ptr, num_Elements)) )
325 throw DataException(msgPrefix+"put(Elements_Owner)");
326
327 // Elements_Color
328 if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
329 throw DataException(msgPrefix+"add_var(Elements_Color)");
330 int_ptr = &mesh->Elements->Color[0];
331 if (! (ids->put(int_ptr, num_Elements)) )
332 throw DataException(msgPrefix+"put(Elements_Color)");
333
334 // Elements_Nodes
335 if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
336 throw DataException(msgPrefix+"add_var(Elements_Nodes)");
337 if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
338 throw DataException(msgPrefix+"put(Elements_Nodes)");
339
340 }
341
342 // // // // // Face_Elements // // // // //
343
344 if (num_FaceElements>0) {
345
346 // FaceElements_Id
347 if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
348 throw DataException(msgPrefix+"add_var(FaceElements_Id)");
349 int_ptr = &mesh->FaceElements->Id[0];
350 if (! (ids->put(int_ptr, num_FaceElements)) )
351 throw DataException(msgPrefix+"put(FaceElements_Id)");
352
353 // FaceElements_Tag
354 if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
355 throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
356 int_ptr = &mesh->FaceElements->Tag[0];
357 if (! (ids->put(int_ptr, num_FaceElements)) )
358 throw DataException(msgPrefix+"put(FaceElements_Tag)");
359
360 // FaceElements_Owner
361 if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
362 throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
363 int_ptr = &mesh->FaceElements->Owner[0];
364 if (! (ids->put(int_ptr, num_FaceElements)) )
365 throw DataException(msgPrefix+"put(FaceElements_Owner)");
366
367 // FaceElements_Color
368 if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
369 throw DataException(msgPrefix+"add_var(FaceElements_Color)");
370 int_ptr = &mesh->FaceElements->Color[0];
371 if (! (ids->put(int_ptr, num_FaceElements)) )
372 throw DataException(msgPrefix+"put(FaceElements_Color)");
373
374 // FaceElements_Nodes
375 if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
376 throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
377 if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
378 throw DataException(msgPrefix+"put(FaceElements_Nodes)");
379
380 }
381
382 // // // // // Points // // // // //
383
384 if (num_Points>0) {
385
386 fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");
387
388 // Points_Id
389 if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
390 throw DataException(msgPrefix+"add_var(Points_Id)");
391 int_ptr = &mesh->Points->Id[0];
392 if (! (ids->put(int_ptr, num_Points)) )
393 throw DataException(msgPrefix+"put(Points_Id)");
394
395 // Points_Tag
396 if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
397 throw DataException(msgPrefix+"add_var(Points_Tag)");
398 int_ptr = &mesh->Points->Tag[0];
399 if (! (ids->put(int_ptr, num_Points)) )
400 throw DataException(msgPrefix+"put(Points_Tag)");
401
402 // Points_Owner
403 if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
404 throw DataException(msgPrefix+"add_var(Points_Owner)");
405 int_ptr = &mesh->Points->Owner[0];
406 if (! (ids->put(int_ptr, num_Points)) )
407 throw DataException(msgPrefix+"put(Points_Owner)");
408
409 // Points_Color
410 if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
411 throw DataException(msgPrefix+"add_var(Points_Color)");
412 int_ptr = &mesh->Points->Color[0];
413 if (! (ids->put(int_ptr, num_Points)) )
414 throw DataException(msgPrefix+"put(Points_Color)");
415
416 // Points_Nodes
417 // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
418 if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
419 throw DataException(msgPrefix+"add_var(Points_Nodes)");
420 if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
421 throw DataException(msgPrefix+"put(Points_Nodes)");
422
423 }
424
425 // // // // // TagMap // // // // //
426
427 if (num_Tags>0) {
428
429 // Temp storage to gather node IDs
430 int *Tags_keys = TMPMEMALLOC(num_Tags, int);
431 char name_temp[4096];
432
433 /* Copy tag data into temp arrays */
434 tag_map = mesh->TagMap;
435 if (tag_map) {
436 int i = 0;
437 while (tag_map) {
438 Tags_keys[i++] = tag_map->tag_key;
439 tag_map=tag_map->next;
440 }
441 }
442
443 // Tags_keys
444 if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
445 throw DataException(msgPrefix+"add_var(Tags_keys)");
446 int_ptr = &Tags_keys[0];
447 if (! (ids->put(int_ptr, num_Tags)) )
448 throw DataException(msgPrefix+"put(Tags_keys)");
449
450 // Tags_names_*
451 // This is an array of strings, it should be stored as an array but
452 // instead I have hacked in one attribute per string because the NetCDF
453 // manual doesn't tell how to do an array of strings
454 tag_map = mesh->TagMap;
455 if (tag_map) {
456 int i = 0;
457 while (tag_map) {
458 sprintf(name_temp, "Tags_name_%d", i);
459 if (!dataFile.add_att(name_temp, tag_map->name) )
460 throw DataException(msgPrefix+"add_att(Tags_names_XX)");
461 tag_map=tag_map->next;
462 i++;
463 }
464 }
465
466 TMPMEMFREE(Tags_keys);
467 }
468
469 /* Send token to next MPI process so he can take his turn */
470 #ifdef ESYS_MPI
471 if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
472 #endif
473
474 // NetCDF file is closed by destructor of NcFile object
475
476 #else
477 Dudley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
478 #endif /* USE_NETCDF */
479 checkDudleyError();
480 }
481
482 string MeshAdapter::getDescription() const
483 {
484 return "DudleyMesh";
485 }
486
487 string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
488 {
489 FunctionSpaceNamesMapType::iterator loc;
490 loc=m_functionSpaceTypeNames.find(functionSpaceType);
491 if (loc==m_functionSpaceTypeNames.end()) {
492 return "Invalid function space type code.";
493 } else {
494 return loc->second;
495 }
496 }
497
498 bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
499 {
500 FunctionSpaceNamesMapType::iterator loc;
501 loc=m_functionSpaceTypeNames.find(functionSpaceType);
502 return (loc!=m_functionSpaceTypeNames.end());
503 }
504
505 void MeshAdapter::setFunctionSpaceTypeNames()
506 {
507 m_functionSpaceTypeNames.insert
508 (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Dudley_DegreesOfFreedom"));
509 m_functionSpaceTypeNames.insert
510 (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Dudley_ReducedDegreesOfFreedom"));
511 m_functionSpaceTypeNames.insert
512 (FunctionSpaceNamesMapType::value_type(Nodes,"Dudley_Nodes"));
513 m_functionSpaceTypeNames.insert
514 (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Dudley_Reduced_Nodes"));
515 m_functionSpaceTypeNames.insert
516 (FunctionSpaceNamesMapType::value_type(Elements,"Dudley_Elements"));
517 m_functionSpaceTypeNames.insert
518 (FunctionSpaceNamesMapType::value_type(ReducedElements,"Dudley_Reduced_Elements"));
519 m_functionSpaceTypeNames.insert
520 (FunctionSpaceNamesMapType::value_type(FaceElements,"Dudley_Face_Elements"));
521 m_functionSpaceTypeNames.insert
522 (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Dudley_Reduced_Face_Elements"));
523 m_functionSpaceTypeNames.insert
524 (FunctionSpaceNamesMapType::value_type(Points,"Dudley_Points"));
525 }
526
527 int MeshAdapter::getContinuousFunctionCode() const
528 {
529 return Nodes;
530 }
531 int MeshAdapter::getReducedContinuousFunctionCode() const
532 {
533 return ReducedNodes;
534 }
535
536 int MeshAdapter::getFunctionCode() const
537 {
538 return Elements;
539 }
540 int MeshAdapter::getReducedFunctionCode() const
541 {
542 return ReducedElements;
543 }
544
545 int MeshAdapter::getFunctionOnBoundaryCode() const
546 {
547 return FaceElements;
548 }
549 int MeshAdapter::getReducedFunctionOnBoundaryCode() const
550 {
551 return ReducedFaceElements;
552 }
553
554 int MeshAdapter::getFunctionOnContactZeroCode() const
555 {
556 throw DudleyAdapterException("Dudley does not support contact elements.");
557 }
558
559 int MeshAdapter::getReducedFunctionOnContactZeroCode() const
560 {
561 throw DudleyAdapterException("Dudley does not support contact elements.");
562 }
563
564 int MeshAdapter::getFunctionOnContactOneCode() const
565 {
566 throw DudleyAdapterException("Dudley does not support contact elements.");
567 }
568
569 int MeshAdapter::getReducedFunctionOnContactOneCode() const
570 {
571 throw DudleyAdapterException("Dudley does not support contact elements.");
572 }
573
574 int MeshAdapter::getSolutionCode() const
575 {
576 return DegreesOfFreedom;
577 }
578
579 int MeshAdapter::getReducedSolutionCode() const
580 {
581 return ReducedDegreesOfFreedom;
582 }
583
584 int MeshAdapter::getDiracDeltaFunctionCode() const
585 {
586 return Points;
587 }
588
589 //
590 // return the spatial dimension of the Mesh:
591 //
592 int MeshAdapter::getDim() const
593 {
594 Dudley_Mesh* mesh=m_dudleyMesh.get();
595 int numDim=Dudley_Mesh_getDim(mesh);
596 checkDudleyError();
597 return numDim;
598 }
599
600 //
601 // Return the number of data points summed across all MPI processes
602 //
603 int MeshAdapter::getNumDataPointsGlobal() const
604 {
605 return Dudley_NodeFile_getGlobalNumNodes(m_dudleyMesh.get()->Nodes);
606 }
607
608 //
609 // return the number of data points per sample and the number of samples
610 // needed to represent data on a parts of the mesh.
611 //
612 pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
613 {
614 int numDataPointsPerSample=0;
615 int numSamples=0;
616 Dudley_Mesh* mesh=m_dudleyMesh.get();
617 switch (functionSpaceCode) {
618 case(Nodes):
619 numDataPointsPerSample=1;
620 numSamples=Dudley_NodeFile_getNumNodes(mesh->Nodes);
621 break;
622 case(ReducedNodes):
623 numDataPointsPerSample=1;
624 numSamples=Dudley_NodeFile_getNumReducedNodes(mesh->Nodes);
625 break;
626 case(Elements):
627 if (mesh->Elements!=NULL) {
628 numSamples=mesh->Elements->numElements;
629 numDataPointsPerSample=mesh->Elements->numLocalDim+1/*referenceElementSet->referenceElement->BasisFunctions->numQuadNodes*/;
630 }
631 break;
632 case(ReducedElements):
633 if (mesh->Elements!=NULL) {
634 numSamples=mesh->Elements->numElements;
635 numDataPointsPerSample=(mesh->Elements->numLocalDim==0)?0:1;
636 }
637 break;
638 case(FaceElements):
639 if (mesh->FaceElements!=NULL) {
640 numDataPointsPerSample=mesh->FaceElements->numLocalDim+1/*referenceElementSet->referenceElement->BasisFunctions->numQuadNodes*/;
641 numSamples=mesh->FaceElements->numElements;
642 }
643 break;
644 case(ReducedFaceElements):
645 if (mesh->FaceElements!=NULL) {
646 numDataPointsPerSample=(mesh->FaceElements->numLocalDim==0)?0:1/*referenceElementSet->referenceElementReducedQuadrature->BasisFunctions->numQuadNodes*/;
647 numSamples=mesh->FaceElements->numElements;
648 }
649 break;
650 case(Points):
651 if (mesh->Points!=NULL) {
652 numDataPointsPerSample=1;
653 numSamples=mesh->Points->numElements;
654 }
655 break;
656 case(DegreesOfFreedom):
657 if (mesh->Nodes!=NULL) {
658 numDataPointsPerSample=1;
659 numSamples=Dudley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
660 }
661 break;
662 case(ReducedDegreesOfFreedom):
663 if (mesh->Nodes!=NULL) {
664 numDataPointsPerSample=1;
665 numSamples=Dudley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
666 }
667 break;
668 default:
669 stringstream temp;
670 temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
671 throw DudleyAdapterException(temp.str());
672 break;
673 }
674 return pair<int,int>(numDataPointsPerSample,numSamples);
675 }
676
677 //
678 // adds linear PDE of second order into a given stiffness matrix and right hand side:
679 //
680 void MeshAdapter::addPDEToSystem(
681 AbstractSystemMatrix& mat, escript::Data& rhs,
682 const escript::Data& A, const escript::Data& B, const escript::Data& C,const escript::Data& D,const escript::Data& X,const escript::Data& Y,
683 const escript::Data& d, const escript::Data& y,
684 const escript::Data& d_contact, const escript::Data& y_contact) const
685 {
686 if (!d_contact.isEmpty())
687 {
688 throw DudleyAdapterException("Dudley does not support d_contact");
689 }
690 if (!y_contact.isEmpty())
691 {
692 throw DudleyAdapterException("Dudley does not support y_contact");
693 }
694 SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
695 if (smat==0)
696 {
697 throw DudleyAdapterException("Dudley only accepts its own system matrices");
698 }
699 escriptDataC _rhs=rhs.getDataC();
700 escriptDataC _A =A.getDataC();
701 escriptDataC _B=B.getDataC();
702 escriptDataC _C=C.getDataC();
703 escriptDataC _D=D.getDataC();
704 escriptDataC _X=X.getDataC();
705 escriptDataC _Y=Y.getDataC();
706 escriptDataC _d=d.getDataC();
707 escriptDataC _y=y.getDataC();
708
709 Dudley_Mesh* mesh=m_dudleyMesh.get();
710
711 Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
712 checkDudleyError();
713
714
715 Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
716 checkDudleyError();
717
718
719 checkDudleyError();
720 }
721
722 void MeshAdapter::addPDEToLumpedSystem(
723 escript::Data& mat,
724 const escript::Data& D,
725 const escript::Data& d) const
726 {
727 escriptDataC _mat=mat.getDataC();
728 escriptDataC _D=D.getDataC();
729 escriptDataC _d=d.getDataC();
730
731 Dudley_Mesh* mesh=m_dudleyMesh.get();
732
733 Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
734 checkDudleyError();
735
736 Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
737 checkDudleyError();
738
739 }
740
741
742 //
743 // adds linear PDE of second order into the right hand side only
744 //
745 void MeshAdapter::addPDEToRHS( escript::Data& rhs, const escript::Data& X,const escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
746 {
747 if (!y_contact.isEmpty())
748 {
749 throw DudleyAdapterException("Dudley does not support y_contact");
750 }
751 Dudley_Mesh* mesh=m_dudleyMesh.get();
752
753 escriptDataC _rhs=rhs.getDataC();
754 escriptDataC _X=X.getDataC();
755 escriptDataC _Y=Y.getDataC();
756 escriptDataC _y=y.getDataC();
757
758 Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
759 checkDudleyError();
760
761 Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
762 checkDudleyError();
763
764 checkDudleyError();
765 }
766 //
767 // adds PDE of second order into a transport problem
768 //
769 void MeshAdapter::addPDEToTransportProblem(
770 AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
771 const escript::Data& A, const escript::Data& B, const escript::Data& C,
772 const escript::Data& D,const escript::Data& X,const escript::Data& Y,
773 const escript::Data& d, const escript::Data& y,
774 const escript::Data& d_contact,const escript::Data& y_contact) const
775 {
776 if (!d_contact.isEmpty())
777 {
778 throw DudleyAdapterException("Dudley does not support d_contact");
779 }
780 if (!y_contact.isEmpty())
781 {
782 throw DudleyAdapterException("Dudley does not support y_contact");
783 }
784 TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
785 if (tpa==0)
786 {
787 throw DudleyAdapterException("Dudley only accepts its own Transportproblems");
788 }
789 DataTypes::ShapeType shape;
790 source.expand();
791 escriptDataC _source=source.getDataC();
792 escriptDataC _M=M.getDataC();
793 escriptDataC _A=A.getDataC();
794 escriptDataC _B=B.getDataC();
795 escriptDataC _C=C.getDataC();
796 escriptDataC _D=D.getDataC();
797 escriptDataC _X=X.getDataC();
798 escriptDataC _Y=Y.getDataC();
799 escriptDataC _d=d.getDataC();
800 escriptDataC _y=y.getDataC();
801 escriptDataC _d_contact=d_contact.getDataC();
802 escriptDataC _y_contact=y_contact.getDataC();
803
804 Dudley_Mesh* mesh=m_dudleyMesh.get();
805 Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
806
807 Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
808 checkDudleyError();
809
810 Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
811 checkDudleyError();
812
813 Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
814 checkDudleyError();
815
816 checkDudleyError();
817 }
818
819 //
820 // interpolates data between different function spaces:
821 //
822 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
823 {
824 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
825 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
826 if (inDomain!=*this)
827 throw DudleyAdapterException("Error - Illegal domain of interpolant.");
828 if (targetDomain!=*this)
829 throw DudleyAdapterException("Error - Illegal domain of interpolation target.");
830
831 Dudley_Mesh* mesh=m_dudleyMesh.get();
832 escriptDataC _target=target.getDataC();
833 escriptDataC _in=in.getDataC();
834 switch(in.getFunctionSpace().getTypeCode()) {
835 case(Nodes):
836 switch(target.getFunctionSpace().getTypeCode()) {
837 case(Nodes):
838 case(ReducedNodes):
839 case(DegreesOfFreedom):
840 case(ReducedDegreesOfFreedom):
841 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
842 break;
843 case(Elements):
844 case(ReducedElements):
845 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
846 break;
847 case(FaceElements):
848 case(ReducedFaceElements):
849 Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
850 break;
851 case(Points):
852 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
853 break;
854 default:
855 stringstream temp;
856 temp << "Error - Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
857 throw DudleyAdapterException(temp.str());
858 break;
859 }
860 break;
861 case(ReducedNodes):
862 switch(target.getFunctionSpace().getTypeCode()) {
863 case(Nodes):
864 case(ReducedNodes):
865 case(DegreesOfFreedom):
866 case(ReducedDegreesOfFreedom):
867 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
868 break;
869 case(Elements):
870 case(ReducedElements):
871 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
872 break;
873 case(FaceElements):
874 case(ReducedFaceElements):
875 Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
876 break;
877 case(Points):
878 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
879 break;
880 default:
881 stringstream temp;
882 temp << "Error - Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
883 throw DudleyAdapterException(temp.str());
884 break;
885 }
886 break;
887 case(Elements):
888 if (target.getFunctionSpace().getTypeCode()==Elements) {
889 Dudley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
890 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
891 Dudley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
892 } else {
893 throw DudleyAdapterException("Error - No interpolation with data on elements possible.");
894 }
895 break;
896 case(ReducedElements):
897 if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
898 Dudley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
899 } else {
900 throw DudleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
901 }
902 break;
903 case(FaceElements):
904 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
905 Dudley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
906 } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
907 Dudley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
908 } else {
909 throw DudleyAdapterException("Error - No interpolation with data on face elements possible.");
910 }
911 break;
912 case(ReducedFaceElements):
913 if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
914 Dudley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
915 } else {
916 throw DudleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
917 }
918 break;
919 case(Points):
920 if (target.getFunctionSpace().getTypeCode()==Points) {
921 Dudley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
922 } else {
923 throw DudleyAdapterException("Error - No interpolation with data on points possible.");
924 }
925 break;
926 case(DegreesOfFreedom):
927 switch(target.getFunctionSpace().getTypeCode()) {
928 case(ReducedDegreesOfFreedom):
929 case(DegreesOfFreedom):
930 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
931 break;
932
933 case(Nodes):
934 case(ReducedNodes):
935 if (getMPISize()>1) {
936 escript::Data temp=escript::Data(in);
937 temp.expand();
938 escriptDataC _in2 = temp.getDataC();
939 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
940 } else {
941 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
942 }
943 break;
944 case(Elements):
945 case(ReducedElements):
946 if (getMPISize()>1) {
947 escript::Data temp=escript::Data( in, continuousFunction(*this) );
948 escriptDataC _in2 = temp.getDataC();
949 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
950 } else {
951 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
952 }
953 break;
954 case(FaceElements):
955 case(ReducedFaceElements):
956 if (getMPISize()>1) {
957 escript::Data temp=escript::Data( in, continuousFunction(*this) );
958 escriptDataC _in2 = temp.getDataC();
959 Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
960
961 } else {
962 Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
963 }
964 break;
965 case(Points):
966 if (getMPISize()>1) {
967 escript::Data temp=escript::Data( in, continuousFunction(*this) );
968 escriptDataC _in2 = temp.getDataC();
969 } else {
970 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
971 }
972 break;
973 default:
974 stringstream temp;
975 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
976 throw DudleyAdapterException(temp.str());
977 break;
978 }
979 break;
980 case(ReducedDegreesOfFreedom):
981 switch(target.getFunctionSpace().getTypeCode()) {
982 case(Nodes):
983 throw DudleyAdapterException("Error - Dudley does not support interpolation from reduced degrees of freedom to mesh nodes.");
984 break;
985 case(ReducedNodes):
986 if (getMPISize()>1) {
987 escript::Data temp=escript::Data(in);
988 temp.expand();
989 escriptDataC _in2 = temp.getDataC();
990 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
991 } else {
992 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
993 }
994 break;
995 case(DegreesOfFreedom):
996 throw DudleyAdapterException("Error - Dudley does not support interpolation from reduced degrees of freedom to degrees of freedom");
997 break;
998 case(ReducedDegreesOfFreedom):
999 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1000 break;
1001 case(Elements):
1002 case(ReducedElements):
1003 if (getMPISize()>1) {
1004 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
1005 escriptDataC _in2 = temp.getDataC();
1006 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1007 } else {
1008 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1009 }
1010 break;
1011 case(FaceElements):
1012 case(ReducedFaceElements):
1013 if (getMPISize()>1) {
1014 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
1015 escriptDataC _in2 = temp.getDataC();
1016 Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1017 } else {
1018 Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1019 }
1020 break;
1021 case(Points):
1022 if (getMPISize()>1) {
1023 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
1024 escriptDataC _in2 = temp.getDataC();
1025 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1026 } else {
1027 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1028 }
1029 break;
1030 default:
1031 stringstream temp;
1032 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1033 throw DudleyAdapterException(temp.str());
1034 break;
1035 }
1036 break;
1037 default:
1038 stringstream temp;
1039 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1040 throw DudleyAdapterException(temp.str());
1041 break;
1042 }
1043 checkDudleyError();
1044 }
1045
1046 //
1047 // copies the locations of sample points into x:
1048 //
1049 void MeshAdapter::setToX(escript::Data& arg) const
1050 {
1051 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1052 if (argDomain!=*this)
1053 throw DudleyAdapterException("Error - Illegal domain of data point locations");
1054 Dudley_Mesh* mesh=m_dudleyMesh.get();
1055 // in case of values node coordinates we can do the job directly:
1056 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1057 escriptDataC _arg=arg.getDataC();
1058 Dudley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1059 } else {
1060 escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1061 escriptDataC _tmp_data=tmp_data.getDataC();
1062 Dudley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1063 // this is then interpolated onto arg:
1064 interpolateOnDomain(arg,tmp_data);
1065 }
1066 checkDudleyError();
1067 }
1068
1069 //
1070 // return the normal vectors at the location of data points as a Data object:
1071 //
1072 void MeshAdapter::setToNormal(escript::Data& normal) const
1073 {
1074 /* const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1075 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1076 if (normalDomain!=*this)
1077 throw DudleyAdapterException("Error - Illegal domain of normal locations");
1078 Dudley_Mesh* mesh=m_dudleyMesh.get();
1079 escriptDataC _normal=normal.getDataC();
1080 switch(normal.getFunctionSpace().getTypeCode()) {
1081 case(Nodes):
1082 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for nodes");
1083 break;
1084 case(ReducedNodes):
1085 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for reduced nodes");
1086 break;
1087 case(Elements):
1088 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for elements");
1089 break;
1090 case(ReducedElements):
1091 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for elements with reduced integration order");
1092 break;
1093 case (FaceElements):
1094 Dudley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1095 break;
1096 case (ReducedFaceElements):
1097 Dudley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1098 break;
1099 case(Points):
1100 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for point elements");
1101 break;
1102 case(DegreesOfFreedom):
1103 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for degrees of freedom.");
1104 break;
1105 case(ReducedDegreesOfFreedom):
1106 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for reduced degrees of freedom.");
1107 break;
1108 default:
1109 stringstream temp;
1110 temp << "Error - Normal Vectors: Dudley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1111 throw DudleyAdapterException(temp.str());
1112 break;
1113 }
1114 checkDudleyError();
1115 }
1116
1117 //
1118 // interpolates data to other domain:
1119 //
1120 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1121 {
1122 const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1123 const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1124 if (targetDomain!=this)
1125 throw DudleyAdapterException("Error - Illegal domain of interpolation target");
1126
1127 throw DudleyAdapterException("Error - Dudley does not allow interpolation across domains yet.");
1128 }
1129
1130 //
1131 // calculates the integral of a function defined of arg:
1132 //
1133 void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1134 {
1135 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1136 if (argDomain!=*this)
1137 throw DudleyAdapterException("Error - Illegal domain of integration kernel");
1138
1139 double blocktimer_start = blocktimer_time();
1140 Dudley_Mesh* mesh=m_dudleyMesh.get();
1141 escriptDataC _temp;
1142 escript::Data temp;
1143 escriptDataC _arg=arg.getDataC();
1144 switch(arg.getFunctionSpace().getTypeCode()) {
1145 case(Nodes):
1146 temp=escript::Data( arg, escript::function(*this) );
1147 _temp=temp.getDataC();
1148 Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1149 break;
1150 case(ReducedNodes):
1151 temp=escript::Data( arg, escript::function(*this) );
1152 _temp=temp.getDataC();
1153 Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1154 break;
1155 case(Elements):
1156 Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1157 break;
1158 case(ReducedElements):
1159 Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1160 break;
1161 case(FaceElements):
1162 Dudley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1163 break;
1164 case(ReducedFaceElements):
1165 Dudley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1166 break;
1167 case(Points):
1168 throw DudleyAdapterException("Error - Integral of data on points is not supported.");
1169 break;
1170 case(DegreesOfFreedom):
1171 temp=escript::Data( arg, escript::function(*this) );
1172 _temp=temp.getDataC();
1173 Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1174 break;
1175 case(ReducedDegreesOfFreedom):
1176 temp=escript::Data( arg, escript::function(*this) );
1177 _temp=temp.getDataC();
1178 Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1179 break;
1180 default:
1181 stringstream temp;
1182 temp << "Error - Integrals: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1183 throw DudleyAdapterException(temp.str());
1184 break;
1185 }
1186 checkDudleyError();
1187 blocktimer_increment("integrate()", blocktimer_start);
1188 }
1189
1190 //
1191 // calculates the gradient of arg:
1192 //
1193 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1194 {
1195 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1196 if (argDomain!=*this)
1197 throw DudleyAdapterException("Error - Illegal domain of gradient argument");
1198 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1199 if (gradDomain!=*this)
1200 throw DudleyAdapterException("Error - Illegal domain of gradient");
1201
1202 Dudley_Mesh* mesh=m_dudleyMesh.get();
1203 escriptDataC _grad=grad.getDataC();
1204 escriptDataC nodeDataC;
1205 escript::Data temp;
1206 if (getMPISize()>1) {
1207 if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1208 temp=escript::Data( arg, continuousFunction(*this) );
1209 nodeDataC = temp.getDataC();
1210 } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1211 temp=escript::Data( arg, reducedContinuousFunction(*this) );
1212 nodeDataC = temp.getDataC();
1213 } else {
1214 nodeDataC = arg.getDataC();
1215 }
1216 } else {
1217 nodeDataC = arg.getDataC();
1218 }
1219 switch(grad.getFunctionSpace().getTypeCode()) {
1220 case(Nodes):
1221 throw DudleyAdapterException("Error - Gradient at nodes is not supported.");
1222 break;
1223 case(ReducedNodes):
1224 throw DudleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1225 break;
1226 case(Elements):
1227 Dudley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1228 break;
1229 case(ReducedElements):
1230 Dudley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1231 break;
1232 case(FaceElements):
1233 Dudley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1234 break;
1235 case(ReducedFaceElements):
1236 Dudley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1237 break;
1238 case(Points):
1239 throw DudleyAdapterException("Error - Gradient at points is not supported.");
1240 break;
1241 case(DegreesOfFreedom):
1242 throw DudleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1243 break;
1244 case(ReducedDegreesOfFreedom):
1245 throw DudleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1246 break;
1247 default:
1248 stringstream temp;
1249 temp << "Error - Gradient: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1250 throw DudleyAdapterException(temp.str());
1251 break;
1252 }
1253 checkDudleyError();
1254 }
1255
1256 //
1257 // returns the size of elements:
1258 //
1259 void MeshAdapter::setToSize(escript::Data& size) const
1260 {
1261 Dudley_Mesh* mesh=m_dudleyMesh.get();
1262 escriptDataC tmp=size.getDataC();
1263 switch(size.getFunctionSpace().getTypeCode()) {
1264 case(Nodes):
1265 throw DudleyAdapterException("Error - Size of nodes is not supported.");
1266 break;
1267 case(ReducedNodes):
1268 throw DudleyAdapterException("Error - Size of reduced nodes is not supported.");
1269 break;
1270 case(Elements):
1271 Dudley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1272 break;
1273 case(ReducedElements):
1274 Dudley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1275 break;
1276 case(FaceElements):
1277 Dudley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1278 break;
1279 case(ReducedFaceElements):
1280 Dudley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1281 break;
1282 case(Points):
1283 throw DudleyAdapterException("Error - Size of point elements is not supported.");
1284 break;
1285 case(DegreesOfFreedom):
1286 throw DudleyAdapterException("Error - Size of degrees of freedom is not supported.");
1287 break;
1288 case(ReducedDegreesOfFreedom):
1289 throw DudleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1290 break;
1291 default:
1292 stringstream temp;
1293 temp << "Error - Element size: Dudley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1294 throw DudleyAdapterException(temp.str());
1295 break;
1296 }
1297 checkDudleyError();
1298 }
1299
1300 //
1301 // sets the location of nodes
1302 //
1303 void MeshAdapter::setNewX(const escript::Data& new_x)
1304 {
1305 Dudley_Mesh* mesh=m_dudleyMesh.get();
1306 escriptDataC tmp;
1307 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1308 if (newDomain!=*this)
1309 throw DudleyAdapterException("Error - Illegal domain of new point locations");
1310 if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1311 tmp = new_x.getDataC();
1312 Dudley_Mesh_setCoordinates(mesh,&tmp);
1313 } else {
1314 escript::Data new_x_inter=escript::Data( new_x, continuousFunction(*this) );
1315 tmp = new_x_inter.getDataC();
1316 Dudley_Mesh_setCoordinates(mesh,&tmp);
1317 }
1318 checkDudleyError();
1319 }
1320
1321 //
1322 // Helper for the save* methods. Extracts optional data variable names and the
1323 // corresponding pointers from python dictionary. Caller must free arrays.
1324 //
1325 void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1326 {
1327 numData = boost::python::extract<int>(arg.attr("__len__")());
1328 /* win32 refactor */
1329 names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1330 data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1331 dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1332
1333 boost::python::list keys=arg.keys();
1334 for (int i=0; i<numData; ++i) {
1335 string n=boost::python::extract<string>(keys[i]);
1336 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1337 if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1338 throw DudleyAdapterException("Error: Data must be defined on same Domain");
1339 data[i] = d.getDataC();
1340 dataPtr[i] = &(data[i]);
1341 names[i] = TMPMEMALLOC(n.length()+1, char);
1342 strcpy(names[i], n.c_str());
1343 }
1344 }
1345
1346 //
1347 // saves mesh and optionally data arrays in openDX format
1348 //
1349 void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1350 {
1351 int num_data;
1352 char **names;
1353 escriptDataC *data;
1354 escriptDataC **ptr_data;
1355
1356 extractArgsFromDict(arg, num_data, names, data, ptr_data);
1357 Dudley_Mesh_saveDX(filename.c_str(), m_dudleyMesh.get(), num_data, names, ptr_data);
1358 checkDudleyError();
1359
1360 /* win32 refactor */
1361 TMPMEMFREE(data);
1362 TMPMEMFREE(ptr_data);
1363 for(int i=0; i<num_data; i++)
1364 {
1365 TMPMEMFREE(names[i]);
1366 }
1367 TMPMEMFREE(names);
1368
1369 return;
1370 }
1371
1372 //
1373 // saves mesh and optionally data arrays in VTK format
1374 //
1375 void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg, const string& metadata, const string& metadata_schema) const
1376 {
1377 int num_data;
1378 char **names;
1379 escriptDataC *data;
1380 escriptDataC **ptr_data;
1381
1382 extractArgsFromDict(arg, num_data, names, data, ptr_data);
1383 Dudley_Mesh_saveVTK(filename.c_str(), m_dudleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());
1384 checkDudleyError();
1385
1386 /* win32 refactor */
1387 TMPMEMFREE(data);
1388 TMPMEMFREE(ptr_data);
1389 for(int i=0; i<num_data; i++)
1390 {
1391 TMPMEMFREE(names[i]);
1392 }
1393 TMPMEMFREE(names);
1394 }
1395
1396 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1397 {
1398 #ifdef ESYS_MPI
1399 index_t myFirstNode=0, myLastNode=0, k=0;
1400 index_t* globalNodeIndex=0;
1401 Dudley_Mesh* mesh_p=m_dudleyMesh.get();
1402 if (fs_code == DUDLEY_REDUCED_NODES)
1403 {
1404 myFirstNode = Dudley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1405 myLastNode = Dudley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1406 globalNodeIndex = Dudley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1407 }
1408 else
1409 {
1410 myFirstNode = Dudley_NodeFile_getFirstNode(mesh_p->Nodes);
1411 myLastNode = Dudley_NodeFile_getLastNode(mesh_p->Nodes);
1412 globalNodeIndex = Dudley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1413 }
1414 k=globalNodeIndex[id];
1415 return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1416 #endif
1417 return true;
1418 }
1419
1420
1421
1422 //
1423 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1424 //
1425 ASM_ptr MeshAdapter::newSystemMatrix(
1426 const int row_blocksize,
1427 const escript::FunctionSpace& row_functionspace,
1428 const int column_blocksize,
1429 const escript::FunctionSpace& column_functionspace,
1430 const int type) const
1431 {
1432 int reduceRowOrder=0;
1433 int reduceColOrder=0;
1434 // is the domain right?
1435 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1436 if (row_domain!=*this)
1437 throw DudleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1438 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1439 if (col_domain!=*this)
1440 throw DudleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1441 // is the function space type right
1442 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1443 reduceRowOrder=0;
1444 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1445 reduceRowOrder=1;
1446 } else {
1447 throw DudleyAdapterException("Error - illegal function space type for system matrix rows.");
1448 }
1449 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1450 reduceColOrder=0;
1451 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1452 reduceColOrder=1;
1453 } else {
1454 throw DudleyAdapterException("Error - illegal function space type for system matrix columns.");
1455 }
1456 // generate matrix:
1457
1458 Paso_SystemMatrixPattern* fsystemMatrixPattern=Dudley_getPattern(getDudley_Mesh(),reduceRowOrder,reduceColOrder);
1459 checkDudleyError();
1460 Paso_SystemMatrix* fsystemMatrix;
1461 int trilinos = 0;
1462 if (trilinos) {
1463 #ifdef TRILINOS
1464 /* Allocation Epetra_VrbMatrix here */
1465 #endif
1466 }
1467 else {
1468 fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1469 }
1470 checkPasoError();
1471 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1472 SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace, column_blocksize,column_functionspace);
1473 return ASM_ptr(sma);
1474 }
1475
1476 //
1477 // creates a TransportProblemAdapter
1478 //
1479 ATP_ptr MeshAdapter::newTransportProblem(
1480 const bool useBackwardEuler,
1481 const int blocksize,
1482 const escript::FunctionSpace& functionspace,
1483 const int type) const
1484 {
1485 int reduceOrder=0;
1486 // is the domain right?
1487 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1488 if (domain!=*this)
1489 throw DudleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1490 // is the function space type right
1491 if (functionspace.getTypeCode()==DegreesOfFreedom) {
1492 reduceOrder=0;
1493 } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1494 reduceOrder=1;
1495 } else {
1496 throw DudleyAdapterException("Error - illegal function space type for system matrix rows.");
1497 }
1498 // generate matrix:
1499
1500 Paso_SystemMatrixPattern* fsystemMatrixPattern=Dudley_getPattern(getDudley_Mesh(),reduceOrder,reduceOrder);
1501 checkDudleyError();
1502 Paso_TransportProblem* transportProblem;
1503 transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1504 checkPasoError();
1505 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1506 AbstractTransportProblem* atp=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1507 return ATP_ptr(atp);
1508 }
1509
1510 //
1511 // vtkObject MeshAdapter::createVtkObject() const
1512 // TODO:
1513 //
1514
1515 //
1516 // returns true if data at the atom_type is considered as being cell centered:
1517 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1518 {
1519 switch(functionSpaceCode) {
1520 case(Nodes):
1521 case(DegreesOfFreedom):
1522 case(ReducedDegreesOfFreedom):
1523 return false;
1524 break;
1525 case(Elements):
1526 case(FaceElements):
1527 case(Points):
1528 case(ReducedElements):
1529 case(ReducedFaceElements):
1530 return true;
1531 break;
1532 default:
1533 stringstream temp;
1534 temp << "Error - Cell: Dudley does not know anything about function space type " << functionSpaceCode;
1535 throw DudleyAdapterException(temp.str());
1536 break;
1537 }
1538 checkDudleyError();
1539 return false;
1540 }
1541
1542 bool
1543 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1544 {
1545 /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1546 class 1: DOF <-> Nodes
1547 class 2: ReducedDOF <-> ReducedNodes
1548 class 3: Points
1549 class 4: Elements
1550 class 5: ReducedElements
1551 class 6: FaceElements
1552 class 7: ReducedFaceElements
1553 class 8: ContactElementZero <-> ContactElementOne
1554 class 9: ReducedContactElementZero <-> ReducedContactElementOne
1555
1556 There is also a set of lines. Interpolation is possible down a line but not between lines.
1557 class 1 and 2 belong to all lines so aren't considered.
1558 line 0: class 3
1559 line 1: class 4,5
1560 line 2: class 6,7
1561 line 3: class 8,9
1562
1563 For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1564 eg hasnodes is true if we have at least one instance of Nodes.
1565 */
1566 if (fs.empty())
1567 {
1568 return false;
1569 }
1570 vector<int> hasclass(10);
1571 vector<int> hasline(4);
1572 bool hasnodes=false;
1573 bool hasrednodes=false;
1574 for (int i=0;i<fs.size();++i)
1575 {
1576 switch(fs[i])
1577 {
1578 case(Nodes): hasnodes=true; // no break is deliberate
1579 case(DegreesOfFreedom):
1580 hasclass[1]=1;
1581 break;
1582 case(ReducedNodes): hasrednodes=true; // no break is deliberate
1583 case(ReducedDegreesOfFreedom):
1584 hasclass[2]=1;
1585 break;
1586 case(Points):
1587 hasline[0]=1;
1588 hasclass[3]=1;
1589 break;
1590 case(Elements):
1591 hasclass[4]=1;
1592 hasline[1]=1;
1593 break;
1594 case(ReducedElements):
1595 hasclass[5]=1;
1596 hasline[1]=1;
1597 break;
1598 case(FaceElements):
1599 hasclass[6]=1;
1600 hasline[2]=1;
1601 break;
1602 case(ReducedFaceElements):
1603 hasclass[7]=1;
1604 hasline[2]=1;
1605 break;
1606 default:
1607 return false;
1608 }
1609 }
1610 int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1611 // fail if we have more than one leaf group
1612
1613 if (totlines>1)
1614 {
1615 return false; // there are at least two branches we can't interpolate between
1616 }
1617 else if (totlines==1)
1618 {
1619 if (hasline[0]==1) // we have points
1620 {
1621 resultcode=Points;
1622 }
1623 else if (hasline[1]==1)
1624 {
1625 if (hasclass[5]==1)
1626 {
1627 resultcode=ReducedElements;
1628 }
1629 else
1630 {
1631 resultcode=Elements;
1632 }
1633 }
1634 else if (hasline[2]==1)
1635 {
1636 if (hasclass[7]==1)
1637 {
1638 resultcode=ReducedFaceElements;
1639 }
1640 else
1641 {
1642 resultcode=FaceElements;
1643 }
1644 }
1645 else // so we must be in line3
1646 {
1647
1648 throw DudleyAdapterException("Programmer Error - choosing between contact elements - we should never get here.");
1649
1650 }
1651 }
1652 else // totlines==0
1653 {
1654 if (hasclass[2]==1)
1655 {
1656 // something from class 2
1657 resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1658 }
1659 else
1660 { // something from class 1
1661 resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1662 }
1663 }
1664 return true;
1665 }
1666
1667 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1668 {
1669 switch(functionSpaceType_source) {
1670 case(Nodes):
1671 switch(functionSpaceType_target) {
1672 case(Nodes):
1673 case(ReducedNodes):
1674 case(ReducedDegreesOfFreedom):
1675 case(DegreesOfFreedom):
1676 case(Elements):
1677 case(ReducedElements):
1678 case(FaceElements):
1679 case(ReducedFaceElements):
1680 case(Points):
1681 return true;
1682 default:
1683 stringstream temp;
1684 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1685 throw DudleyAdapterException(temp.str());
1686 }
1687 break;
1688 case(ReducedNodes):
1689 switch(functionSpaceType_target) {
1690 case(ReducedNodes):
1691 case(ReducedDegreesOfFreedom):
1692 case(Elements):
1693 case(ReducedElements):
1694 case(FaceElements):
1695 case(ReducedFaceElements):
1696 case(Points):
1697 return true;
1698 case(Nodes):
1699 case(DegreesOfFreedom):
1700 return false;
1701 default:
1702 stringstream temp;
1703 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1704 throw DudleyAdapterException(temp.str());
1705 }
1706 break;
1707 case(Elements):
1708 if (functionSpaceType_target==Elements) {
1709 return true;
1710 } else if (functionSpaceType_target==ReducedElements) {
1711 return true;
1712 } else {
1713 return false;
1714 }
1715 case(ReducedElements):
1716 if (functionSpaceType_target==ReducedElements) {
1717 return true;
1718 } else {
1719 return false;
1720 }
1721 case(FaceElements):
1722 if (functionSpaceType_target==FaceElements) {
1723 return true;
1724 } else if (functionSpaceType_target==ReducedFaceElements) {
1725 return true;
1726 } else {
1727 return false;
1728 }
1729 case(ReducedFaceElements):
1730 if (functionSpaceType_target==ReducedFaceElements) {
1731 return true;
1732 } else {
1733 return false;
1734 }
1735 case(Points):
1736 if (functionSpaceType_target==Points) {
1737 return true;
1738 } else {
1739 return false;
1740 }
1741 case(DegreesOfFreedom):
1742 switch(functionSpaceType_target) {
1743 case(ReducedDegreesOfFreedom):
1744 case(DegreesOfFreedom):
1745 case(Nodes):
1746 case(ReducedNodes):
1747 case(Elements):
1748 case(ReducedElements):
1749 case(Points):
1750 case(FaceElements):
1751 case(ReducedFaceElements):
1752 return true;
1753 default:
1754 stringstream temp;
1755 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1756 throw DudleyAdapterException(temp.str());
1757 }
1758 break;
1759 case(ReducedDegreesOfFreedom):
1760 switch(functionSpaceType_target) {
1761 case(ReducedDegreesOfFreedom):
1762 case(ReducedNodes):
1763 case(Elements):
1764 case(ReducedElements):
1765 case(FaceElements):
1766 case(ReducedFaceElements):
1767 case(Points):
1768 return true;
1769 case(Nodes):
1770 case(DegreesOfFreedom):
1771 return false;
1772 default:
1773 stringstream temp;
1774 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1775 throw DudleyAdapterException(temp.str());
1776 }
1777 break;
1778 default:
1779 stringstream temp;
1780 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_source;
1781 throw DudleyAdapterException(temp.str());
1782 break;
1783 }
1784 checkDudleyError();
1785 return false;
1786 }
1787
1788 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1789 {
1790 return false;
1791 }
1792
1793 bool MeshAdapter::operator==(const AbstractDomain& other) const
1794 {
1795 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1796 if (temp!=0) {
1797 return (m_dudleyMesh==temp->m_dudleyMesh);
1798 } else {
1799 return false;
1800 }
1801 }
1802
1803 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1804 {
1805 return !(operator==(other));
1806 }
1807
1808 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1809 {
1810 Dudley_Mesh* mesh=m_dudleyMesh.get();
1811 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1812 checkPasoError();
1813 return out;
1814 }
1815 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1816 {
1817 Dudley_Mesh* mesh=m_dudleyMesh.get();
1818 int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1819 checkPasoError();
1820 return out;
1821 }
1822
1823 escript::Data MeshAdapter::getX() const
1824 {
1825 return continuousFunction(*this).getX();
1826 }
1827
1828 escript::Data MeshAdapter::getNormal() const
1829 {
1830 return functionOnBoundary(*this).getNormal();
1831 }
1832
1833 escript::Data MeshAdapter::getSize() const
1834 {
1835 return escript::function(*this).getSize();
1836 }
1837
1838 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1839 {
1840 int *out = NULL;
1841 Dudley_Mesh* mesh=m_dudleyMesh.get();
1842 switch (functionSpaceType) {
1843 case(Nodes):
1844 out=mesh->Nodes->Id;
1845 break;
1846 case(ReducedNodes):
1847 out=mesh->Nodes->reducedNodesId;
1848 break;
1849 case(Elements):
1850 out=mesh->Elements->Id;
1851 break;
1852 case(ReducedElements):
1853 out=mesh->Elements->Id;
1854 break;
1855 case(FaceElements):
1856 out=mesh->FaceElements->Id;
1857 break;
1858 case(ReducedFaceElements):
1859 out=mesh->FaceElements->Id;
1860 break;
1861 case(Points):
1862 out=mesh->Points->Id;
1863 break;
1864 case(DegreesOfFreedom):
1865 out=mesh->Nodes->degreesOfFreedomId;
1866 break;
1867 case(ReducedDegreesOfFreedom):
1868 out=mesh->Nodes->reducedDegreesOfFreedomId;
1869 break;
1870 default:
1871 stringstream temp;
1872 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1873 throw DudleyAdapterException(temp.str());
1874 break;
1875 }
1876 return out;
1877 }
1878 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1879 {
1880 int out=0;
1881 Dudley_Mesh* mesh=m_dudleyMesh.get();
1882 switch (functionSpaceType) {
1883 case(Nodes):
1884 out=mesh->Nodes->Tag[sampleNo];
1885 break;
1886 case(ReducedNodes):
1887 throw DudleyAdapterException(" Error - ReducedNodes does not support tags.");
1888 break;
1889 case(Elements):
1890 out=mesh->Elements->Tag[sampleNo];
1891 break;
1892 case(ReducedElements):
1893 out=mesh->Elements->Tag[sampleNo];
1894 break;
1895 case(FaceElements):
1896 out=mesh->FaceElements->Tag[sampleNo];
1897 break;
1898 case(ReducedFaceElements):
1899 out=mesh->FaceElements->Tag[sampleNo];
1900 break;
1901 case(Points):
1902 out=mesh->Points->Tag[sampleNo];
1903 break;
1904 case(DegreesOfFreedom):
1905 throw DudleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1906 break;
1907 case(ReducedDegreesOfFreedom):
1908 throw DudleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1909 break;
1910 default:
1911 stringstream temp;
1912 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1913 throw DudleyAdapterException(temp.str());
1914 break;
1915 }
1916 return out;
1917 }
1918
1919
1920 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1921 {
1922 Dudley_Mesh* mesh=m_dudleyMesh.get();
1923 escriptDataC tmp=mask.getDataC();
1924 switch(functionSpaceType) {
1925 case(Nodes):
1926 Dudley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1927 break;
1928 case(ReducedNodes):
1929 throw DudleyAdapterException("Error - ReducedNodes does not support tags");
1930 break;
1931 case(DegreesOfFreedom):
1932 throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
1933 break;
1934 case(ReducedDegreesOfFreedom):
1935 throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1936 break;
1937 case(Elements):
1938 Dudley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1939 break;
1940 case(ReducedElements):
1941 Dudley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1942 break;
1943 case(FaceElements):
1944 Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1945 break;
1946 case(ReducedFaceElements):
1947 Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1948 break;
1949 case(Points):
1950 Dudley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1951 break;
1952 default:
1953 stringstream temp;
1954 temp << "Error - Dudley does not know anything about function space type " << functionSpaceType;
1955 throw DudleyAdapterException(temp.str());
1956 }
1957 checkDudleyError();
1958 return;
1959 }
1960
1961 void MeshAdapter::setTagMap(const string& name, int tag)
1962 {
1963 Dudley_Mesh* mesh=m_dudleyMesh.get();
1964 Dudley_Mesh_addTagMap(mesh, name.c_str(),tag);
1965 checkPasoError();
1966 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1967 }
1968
1969 int MeshAdapter::getTag(const string& name) const
1970 {
1971 Dudley_Mesh* mesh=m_dudleyMesh.get();
1972 int tag=0;
1973 tag=Dudley_Mesh_getTag(mesh, name.c_str());
1974 checkPasoError();
1975 // throwStandardException("MeshAdapter::getTag is not implemented.");
1976 return tag;
1977 }
1978
1979 bool MeshAdapter::isValidTagName(const string& name) const
1980 {
1981 Dudley_Mesh* mesh=m_dudleyMesh.get();
1982 return Dudley_Mesh_isValidTagName(mesh,name.c_str());
1983 }
1984
1985 string MeshAdapter::showTagNames() const
1986 {
1987 stringstream temp;
1988 Dudley_Mesh* mesh=m_dudleyMesh.get();
1989 Dudley_TagMap* tag_map=mesh->TagMap;
1990 while (tag_map) {
1991 temp << tag_map->name;
1992 tag_map=tag_map->next;
1993 if (tag_map) temp << ", ";
1994 }
1995 return temp.str();
1996 }
1997
1998 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
1999 {
2000 Dudley_Mesh* mesh=m_dudleyMesh.get();
2001 dim_t numTags=0;
2002 switch(functionSpaceCode) {
2003 case(Nodes):
2004 numTags=mesh->Nodes->numTagsInUse;
2005 break;
2006 case(ReducedNodes):
2007 throw DudleyAdapterException("Error - ReducedNodes does not support tags");
2008 break;
2009 case(DegreesOfFreedom):
2010 throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
2011 break;
2012 case(ReducedDegreesOfFreedom):
2013 throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2014 break;
2015 case(Elements):
2016 case(ReducedElements):
2017 numTags=mesh->Elements->numTagsInUse;
2018 break;
2019 case(FaceElements):
2020 case(ReducedFaceElements):
2021 numTags=mesh->FaceElements->numTagsInUse;
2022 break;
2023 case(Points):
2024 numTags=mesh->Points->numTagsInUse;
2025 break;
2026 default:
2027 stringstream temp;
2028 temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2029 throw DudleyAdapterException(temp.str());
2030 }
2031 return numTags;
2032 }
2033
2034 const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2035 {
2036 Dudley_Mesh* mesh=m_dudleyMesh.get();
2037 index_t* tags=NULL;
2038 switch(functionSpaceCode) {
2039 case(Nodes):
2040 tags=mesh->Nodes->tagsInUse;
2041 break;
2042 case(ReducedNodes):
2043 throw DudleyAdapterException("Error - ReducedNodes does not support tags");
2044 break;
2045 case(DegreesOfFreedom):
2046 throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
2047 break;
2048 case(ReducedDegreesOfFreedom):
2049 throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2050 break;
2051 case(Elements):
2052 case(ReducedElements):
2053 tags=mesh->Elements->tagsInUse;
2054 break;
2055 case(FaceElements):
2056 case(ReducedFaceElements):
2057 tags=mesh->FaceElements->tagsInUse;
2058 break;
2059 case(Points):
2060 tags=mesh->Points->tagsInUse;
2061 break;
2062 default:
2063 stringstream temp;
2064 temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2065 throw DudleyAdapterException(temp.str());
2066 }
2067 return tags;
2068 }
2069
2070
2071 bool MeshAdapter::canTag(int functionSpaceCode) const
2072 {
2073 switch(functionSpaceCode) {
2074 case(Nodes):
2075 case(Elements):
2076 case(ReducedElements):
2077 case(FaceElements):
2078 case(ReducedFaceElements):
2079 case(Points):
2080 return true;
2081 case(ReducedNodes):
2082 case(DegreesOfFreedom):
2083 case(ReducedDegreesOfFreedom):
2084 return false;
2085 default:
2086 return false;
2087 }
2088 }
2089
2090 AbstractDomain::StatusType MeshAdapter::getStatus() const
2091 {
2092 Dudley_Mesh* mesh=m_dudleyMesh.get();
2093 return Dudley_Mesh_getStatus(mesh);
2094 }
2095
2096 int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2097 {
2098
2099 Dudley_Mesh* mesh=m_dudleyMesh.get();
2100 int order =-1;
2101 switch(functionSpaceCode) {
2102 case(Nodes):
2103 case(DegreesOfFreedom):
2104 order=mesh->approximationOrder;
2105 break;
2106 case(ReducedNodes):
2107 case(ReducedDegreesOfFreedom):
2108 order=mesh->reducedApproximationOrder;
2109 break;
2110 case(Elements):
2111 case(FaceElements):
2112 case(Points):
2113 order=mesh->integrationOrder;
2114 break;
2115 case(ReducedElements):
2116 case(ReducedFaceElements):
2117 order=mesh->reducedIntegrationOrder;
2118 break;
2119 default:
2120 stringstream temp;
2121 temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2122 throw DudleyAdapterException(temp.str());
2123 }
2124 return order;
2125 }
2126
2127
2128 bool MeshAdapter::supportsContactElements() const
2129 {
2130 return false;
2131 }
2132
2133 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26