/[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 3344 - (show annotations)
Thu Nov 11 23:26:52 2010 UTC (8 years, 2 months ago) by caltinay
File size: 72559 byte(s)
Phew!
-escript, finley, and dudley now uses weipa's saveVTK implementation
-moved tests from finley to weipa accordingly; dudley still to do
-rebaselined all test files
-fixed a few issues in weipa.saveVTK, e.g. saving metadata without schema
-added a deprecation warning to esys.escript.util.saveVTK
-todo: change doco, tests and other places to use weipa.saveVTK


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26