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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6009 - (show annotations)
Wed Mar 2 04:13:26 2016 UTC (3 years, 1 month ago) by caltinay
File size: 69620 byte(s)
Much needed sync with trunk...

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include "MeshAdapter.h"
18
19 #include <escript/Data.h>
20 #include <escript/DataFactory.h>
21 #include <escript/Random.h>
22 #include <escript/SolverOptions.h>
23
24 #include <paso/SystemMatrix.h>
25 #include <paso/Transport.h>
26
27 #ifdef USE_NETCDF
28 #include <netcdfcpp.h>
29 #endif
30 #include <boost/python/tuple.hpp>
31
32 using namespace std;
33 using namespace escript;
34 namespace bp = boost::python;
35
36 namespace dudley {
37
38 //
39 // define the static constants
40 MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
41 const int MeshAdapter::DegreesOfFreedom=DUDLEY_DEGREES_OF_FREEDOM;
42 const int MeshAdapter::ReducedDegreesOfFreedom=DUDLEY_REDUCED_DEGREES_OF_FREEDOM;
43 const int MeshAdapter::Nodes=DUDLEY_NODES;
44 const int MeshAdapter::ReducedNodes=DUDLEY_REDUCED_NODES;
45 const int MeshAdapter::Elements=DUDLEY_ELEMENTS;
46 const int MeshAdapter::ReducedElements=DUDLEY_REDUCED_ELEMENTS;
47 const int MeshAdapter::FaceElements=DUDLEY_FACE_ELEMENTS;
48 const int MeshAdapter::ReducedFaceElements=DUDLEY_REDUCED_FACE_ELEMENTS;
49 const int MeshAdapter::Points=DUDLEY_POINTS;
50
51 MeshAdapter::MeshAdapter(Dudley_Mesh* dudleyMesh)
52 {
53 setFunctionSpaceTypeNames();
54 //
55 // need to use a null_deleter as Dudley_Mesh_free deletes the pointer
56 // for us.
57 m_dudleyMesh.reset(dudleyMesh,null_deleter());
58 }
59
60 //
61 // The copy constructor should just increment the use count
62 MeshAdapter::MeshAdapter(const MeshAdapter& in):
63 m_dudleyMesh(in.m_dudleyMesh)
64 {
65 setFunctionSpaceTypeNames();
66 }
67
68 MeshAdapter::~MeshAdapter()
69 {
70 //
71 // I hope the case for the pointer being zero has been taken care of.
72 // cout << "In MeshAdapter destructor." << endl;
73 if (m_dudleyMesh.unique()) {
74 Dudley_Mesh_free(m_dudleyMesh.get());
75 }
76 }
77
78 escript::JMPI MeshAdapter::getMPI() const
79 {
80 return m_dudleyMesh.get()->MPIInfo;
81 }
82
83 int MeshAdapter::getMPISize() const
84 {
85 return m_dudleyMesh.get()->MPIInfo->size;
86 }
87
88 int MeshAdapter::getMPIRank() const
89 {
90 return m_dudleyMesh.get()->MPIInfo->rank;
91 }
92
93 void MeshAdapter::MPIBarrier() const
94 {
95 #ifdef ESYS_MPI
96 MPI_Barrier(m_dudleyMesh.get()->MPIInfo->comm);
97 #endif
98 }
99
100 bool MeshAdapter::onMasterProcessor() const
101 {
102 return m_dudleyMesh.get()->MPIInfo->rank == 0;
103 }
104
105 MPI_Comm MeshAdapter::getMPIComm() const
106 {
107 return m_dudleyMesh->MPIInfo->comm;
108 }
109
110 Dudley_Mesh* MeshAdapter::getDudley_Mesh() const
111 {
112 return m_dudleyMesh.get();
113 }
114
115 void MeshAdapter::write(const string& fileName) const
116 {
117 char *fName = (fileName.size()+1>0) ? new char[fileName.size()+1] : (char*)NULL;
118 strcpy(fName,fileName.c_str());
119 Dudley_Mesh_write(m_dudleyMesh.get(),fName);
120 delete[] fName;
121 }
122
123 void MeshAdapter::Print_Mesh_Info(const bool full) const
124 {
125 Dudley_PrintMesh_Info(m_dudleyMesh.get(), full);
126 }
127
128 void MeshAdapter::dump(const string& fileName) const
129 {
130 #ifdef USE_NETCDF
131 const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
132 NcVar *ids;
133 int *int_ptr;
134 Dudley_Mesh *mesh = m_dudleyMesh.get();
135 Dudley_TagMap* tag_map;
136 int num_Tags = 0;
137 int mpi_size = mesh->MPIInfo->size;
138 int mpi_rank = mesh->MPIInfo->rank;
139 int numDim = mesh->Nodes->numDim;
140 int numNodes = mesh->Nodes->numNodes;
141 int num_Elements = mesh->Elements->numElements;
142 int num_FaceElements = mesh->FaceElements->numElements;
143 int num_Points = mesh->Points->numElements;
144 int num_Elements_numNodes = mesh->Elements->numNodes;
145 int num_FaceElements_numNodes = mesh->FaceElements->numNodes;
146 #ifdef ESYS_MPI
147 MPI_Status status;
148 #endif
149
150 /* Incoming token indicates it's my turn to write */
151 #ifdef ESYS_MPI
152 if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
153 #endif
154
155 string newFileName(mesh->MPIInfo->appendRankToFileName(fileName));
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.c_str(), 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 = new int[num_Tags];
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 delete[] 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 throw DudleyException("MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
481 #endif /* USE_NETCDF */
482 }
483
484 string MeshAdapter::getDescription() const
485 {
486 return "DudleyMesh";
487 }
488
489 string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
490 {
491 FunctionSpaceNamesMapType::iterator loc;
492 loc=m_functionSpaceTypeNames.find(functionSpaceType);
493 if (loc==m_functionSpaceTypeNames.end()) {
494 return "Invalid function space type code.";
495 } else {
496 return loc->second;
497 }
498 }
499
500 bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
501 {
502 FunctionSpaceNamesMapType::iterator loc;
503 loc=m_functionSpaceTypeNames.find(functionSpaceType);
504 return (loc!=m_functionSpaceTypeNames.end());
505 }
506
507 void MeshAdapter::setFunctionSpaceTypeNames()
508 {
509 m_functionSpaceTypeNames.insert
510 (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Dudley_DegreesOfFreedom [Solution(domain)]"));
511 m_functionSpaceTypeNames.insert
512 (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Dudley_ReducedDegreesOfFreedom [ReducedSolution(domain)]"));
513 m_functionSpaceTypeNames.insert
514 (FunctionSpaceNamesMapType::value_type(Nodes,"Dudley_Nodes [ContinuousFunction(domain)]"));
515 m_functionSpaceTypeNames.insert
516 (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Dudley_Reduced_Nodes [ReducedContinuousFunction(domain)]"));
517 m_functionSpaceTypeNames.insert
518 (FunctionSpaceNamesMapType::value_type(Elements,"Dudley_Elements [Function(domain)]"));
519 m_functionSpaceTypeNames.insert
520 (FunctionSpaceNamesMapType::value_type(ReducedElements,"Dudley_Reduced_Elements [ReducedFunction(domain)]"));
521 m_functionSpaceTypeNames.insert
522 (FunctionSpaceNamesMapType::value_type(FaceElements,"Dudley_Face_Elements [FunctionOnBoundary(domain)]"));
523 m_functionSpaceTypeNames.insert
524 (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Dudley_Reduced_Face_Elements [ReducedFunctionOnBoundary(domain)]"));
525 m_functionSpaceTypeNames.insert
526 (FunctionSpaceNamesMapType::value_type(Points,"Dudley_Points [DiracDeltaFunctions(domain)]"));
527 }
528
529 int MeshAdapter::getContinuousFunctionCode() const
530 {
531 return Nodes;
532 }
533 int MeshAdapter::getReducedContinuousFunctionCode() const
534 {
535 return ReducedNodes;
536 }
537
538 int MeshAdapter::getFunctionCode() const
539 {
540 return Elements;
541 }
542 int MeshAdapter::getReducedFunctionCode() const
543 {
544 return ReducedElements;
545 }
546
547 int MeshAdapter::getFunctionOnBoundaryCode() const
548 {
549 return FaceElements;
550 }
551 int MeshAdapter::getReducedFunctionOnBoundaryCode() const
552 {
553 return ReducedFaceElements;
554 }
555
556 int MeshAdapter::getFunctionOnContactZeroCode() const
557 {
558 throw DudleyException("Dudley does not support contact elements.");
559 }
560
561 int MeshAdapter::getReducedFunctionOnContactZeroCode() const
562 {
563 throw DudleyException("Dudley does not support contact elements.");
564 }
565
566 int MeshAdapter::getFunctionOnContactOneCode() const
567 {
568 throw DudleyException("Dudley does not support contact elements.");
569 }
570
571 int MeshAdapter::getReducedFunctionOnContactOneCode() const
572 {
573 throw DudleyException("Dudley does not support contact elements.");
574 }
575
576 int MeshAdapter::getSolutionCode() const
577 {
578 return DegreesOfFreedom;
579 }
580
581 int MeshAdapter::getReducedSolutionCode() const
582 {
583 return ReducedDegreesOfFreedom;
584 }
585
586 int MeshAdapter::getDiracDeltaFunctionsCode() const
587 {
588 return Points;
589 }
590
591 //
592 // return the spatial dimension of the Mesh:
593 //
594 int MeshAdapter::getDim() const
595 {
596 Dudley_Mesh* mesh=m_dudleyMesh.get();
597 int numDim=Dudley_Mesh_getDim(mesh);
598 return numDim;
599 }
600
601 //
602 // Return the number of data points summed across all MPI processes
603 //
604 int MeshAdapter::getNumDataPointsGlobal() const
605 {
606 return Dudley_NodeFile_getGlobalNumNodes(m_dudleyMesh.get()->Nodes);
607 }
608
609 //
610 // return the number of data points per sample and the number of samples
611 // needed to represent data on a parts of the mesh.
612 //
613 pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
614 {
615 int numDataPointsPerSample=0;
616 int numSamples=0;
617 Dudley_Mesh* mesh=m_dudleyMesh.get();
618 switch (functionSpaceCode) {
619 case(Nodes):
620 numDataPointsPerSample=1;
621 numSamples=Dudley_NodeFile_getNumNodes(mesh->Nodes);
622 break;
623 case(ReducedNodes):
624 numDataPointsPerSample=1;
625 numSamples=Dudley_NodeFile_getNumReducedNodes(mesh->Nodes);
626 break;
627 case(Elements):
628 if (mesh->Elements!=NULL) {
629 numSamples=mesh->Elements->numElements;
630 numDataPointsPerSample=mesh->Elements->numLocalDim+1/*referenceElementSet->referenceElement->BasisFunctions->numQuadNodes*/;
631 }
632 break;
633 case(ReducedElements):
634 if (mesh->Elements!=NULL) {
635 numSamples=mesh->Elements->numElements;
636 numDataPointsPerSample=(mesh->Elements->numLocalDim==0)?0:1;
637 }
638 break;
639 case(FaceElements):
640 if (mesh->FaceElements!=NULL) {
641 numDataPointsPerSample=mesh->FaceElements->numLocalDim+1/*referenceElementSet->referenceElement->BasisFunctions->numQuadNodes*/;
642 numSamples=mesh->FaceElements->numElements;
643 }
644 break;
645 case(ReducedFaceElements):
646 if (mesh->FaceElements!=NULL) {
647 numDataPointsPerSample=(mesh->FaceElements->numLocalDim==0)?0:1/*referenceElementSet->referenceElementReducedQuadrature->BasisFunctions->numQuadNodes*/;
648 numSamples=mesh->FaceElements->numElements;
649 }
650 break;
651 case(Points):
652 if (mesh->Points!=NULL) {
653 numDataPointsPerSample=1;
654 numSamples=mesh->Points->numElements;
655 }
656 break;
657 case(DegreesOfFreedom):
658 if (mesh->Nodes!=NULL) {
659 numDataPointsPerSample=1;
660 numSamples=Dudley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
661 }
662 break;
663 case(ReducedDegreesOfFreedom):
664 if (mesh->Nodes!=NULL) {
665 numDataPointsPerSample=1;
666 numSamples=Dudley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
667 }
668 break;
669 default:
670 stringstream temp;
671 temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
672 throw DudleyException(temp.str());
673 break;
674 }
675 return pair<int,int>(numDataPointsPerSample,numSamples);
676 }
677
678 //
679 // adds linear PDE of second order into a given stiffness matrix and right hand side:
680 //
681 void MeshAdapter::addPDEToSystem(AbstractSystemMatrix& mat, escript::Data& rhs,
682 const escript::Data& A, const escript::Data& B,
683 const escript::Data& C, const escript::Data& D,
684 const escript::Data& X, const escript::Data& Y,
685 const escript::Data& d, const escript::Data& y,
686 const escript::Data& d_contact,
687 const escript::Data& y_contact,
688 const escript::Data& d_dirac,
689 const escript::Data& y_dirac) const
690 {
691 if (!d_contact.isEmpty() || !y_contact.isEmpty())
692 throw DudleyException("Dudley does not support contact elements");
693
694 paso::SystemMatrix* smat = dynamic_cast<paso::SystemMatrix*>(&mat);
695 if (smat) {
696 paso::SystemMatrix_ptr S(smat->shared_from_this());
697 Dudley_Mesh* mesh = m_dudleyMesh.get();
698
699 Assemble_PDE(mesh->Nodes, mesh->Elements, S, rhs, A, B, C, D, X, Y);
700
701 Assemble_PDE(mesh->Nodes, mesh->FaceElements, S, rhs, escript::Data(),
702 escript::Data(), escript::Data(), d, escript::Data(), y);
703
704 Assemble_PDE(mesh->Nodes, mesh->Points, S, rhs, escript::Data(),
705 escript::Data(), escript::Data(), d_dirac,
706 escript::Data(), y_dirac);
707 return;
708 }
709 throw DudleyException("Dudley only accepts Paso system matrices");
710 }
711
712 void MeshAdapter::addPDEToLumpedSystem(escript::Data& mat,
713 const escript::Data& D,
714 const escript::Data& d,
715 const escript::Data& d_dirac,
716 bool useHRZ) const
717 {
718 Dudley_Mesh* mesh=m_dudleyMesh.get();
719
720 Assemble_LumpedSystem(mesh->Nodes,mesh->Elements, mat, D, useHRZ);
721 Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements, mat, d, useHRZ);
722 Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements, mat, d_dirac, useHRZ);
723 }
724
725
726 //
727 // adds linear PDE of second order into the right hand side only
728 //
729 void MeshAdapter::addPDEToRHS(escript::Data& rhs, const escript::Data& X,const escript::Data& Y, const escript::Data& y, const escript::Data& y_contact, const escript::Data& y_dirac) const
730 {
731 if (!y_contact.isEmpty())
732 {
733 throw DudleyException("Dudley does not support y_contact");
734 }
735 Dudley_Mesh* mesh=m_dudleyMesh.get();
736
737 Assemble_PDE(mesh->Nodes, mesh->Elements, escript::ASM_ptr(), rhs,
738 escript::Data(), escript::Data(), escript::Data(),
739 escript::Data(), X, Y);
740
741 Assemble_PDE(mesh->Nodes, mesh->FaceElements, escript::ASM_ptr(), rhs,
742 escript::Data(), escript::Data(), escript::Data(),
743 escript::Data(), escript::Data(), y);
744
745 Assemble_PDE(mesh->Nodes, mesh->Points, escript::ASM_ptr(), rhs,
746 escript::Data(), escript::Data(), escript::Data(),
747 escript::Data(), escript::Data(), y_dirac);
748 }
749 //
750 // adds PDE of second order into a transport problem
751 //
752 void MeshAdapter::addPDEToTransportProblem(
753 AbstractTransportProblem& tp, escript::Data& source,
754 const escript::Data& M, const escript::Data& A, const escript::Data& B,
755 const escript::Data& C, const escript::Data& D,const escript::Data& X,
756 const escript::Data& Y, const escript::Data& d, const escript::Data& y,
757 const escript::Data& d_contact,const escript::Data& y_contact,
758 const escript::Data& d_dirac, const escript::Data& y_dirac) const
759 {
760 if (!d_contact.isEmpty())
761 {
762 throw DudleyException("Dudley does not support d_contact");
763 }
764 if (!y_contact.isEmpty())
765 {
766 throw DudleyException("Dudley does not support y_contact");
767 }
768 paso::TransportProblem* ptp = dynamic_cast<paso::TransportProblem*>(&tp);
769 if (!ptp)
770 {
771 throw DudleyException("Dudley only accepts Paso transport problems");
772 }
773 DataTypes::ShapeType shape;
774 source.expand();
775
776 Dudley_Mesh* mesh=m_dudleyMesh.get();
777
778 Assemble_PDE(mesh->Nodes, mesh->Elements, ptp->borrowMassMatrix(), source,
779 escript::Data(), escript::Data(), escript::Data(), M,
780 escript::Data(), escript::Data());
781
782 Assemble_PDE(mesh->Nodes, mesh->Elements, ptp->borrowTransportMatrix(),
783 source, A, B, C, D, X, Y);
784
785 Assemble_PDE(mesh->Nodes, mesh->FaceElements, ptp->borrowTransportMatrix(),
786 source, escript::Data(), escript::Data(), escript::Data(), d,
787 escript::Data(), y);
788
789 Assemble_PDE(mesh->Nodes, mesh->Points, ptp->borrowTransportMatrix(),
790 source, escript::Data(), escript::Data(), escript::Data(),
791 d_dirac, escript::Data(), y_dirac);
792 }
793
794 //
795 // interpolates data between different function spaces:
796 //
797 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
798 {
799 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
800 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
801 if (inDomain!=*this)
802 throw DudleyException("Illegal domain of interpolant.");
803 if (targetDomain!=*this)
804 throw DudleyException("Illegal domain of interpolation target.");
805
806 Dudley_Mesh* mesh=m_dudleyMesh.get();
807 switch(in.getFunctionSpace().getTypeCode()) {
808 case(Nodes):
809 switch(target.getFunctionSpace().getTypeCode()) {
810 case(Nodes):
811 case(ReducedNodes):
812 case(DegreesOfFreedom):
813 case(ReducedDegreesOfFreedom):
814 Assemble_CopyNodalData(mesh->Nodes,&target,&in);
815 break;
816 case(Elements):
817 case(ReducedElements):
818 Assemble_interpolate(mesh->Nodes,mesh->Elements,&in,&target);
819 break;
820 case(FaceElements):
821 case(ReducedFaceElements):
822 Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&in,&target);
823 break;
824 case(Points):
825 Assemble_interpolate(mesh->Nodes,mesh->Points,&in,&target);
826 break;
827 default:
828 stringstream temp;
829 temp << "Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
830 throw DudleyException(temp.str());
831 break;
832 }
833 break;
834 case(ReducedNodes):
835 switch(target.getFunctionSpace().getTypeCode()) {
836 case(Nodes):
837 case(ReducedNodes):
838 case(DegreesOfFreedom):
839 case(ReducedDegreesOfFreedom):
840 Assemble_CopyNodalData(mesh->Nodes,&target,&in);
841 break;
842 case(Elements):
843 case(ReducedElements):
844 Assemble_interpolate(mesh->Nodes,mesh->Elements,&in,&target);
845 break;
846 case(FaceElements):
847 case(ReducedFaceElements):
848 Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&in,&target);
849 break;
850 case(Points):
851 Assemble_interpolate(mesh->Nodes,mesh->Points,&in,&target);
852 break;
853 default:
854 stringstream temp;
855 temp << "Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
856 throw DudleyException(temp.str());
857 break;
858 }
859 break;
860 case(Elements):
861 if (target.getFunctionSpace().getTypeCode()==Elements) {
862 Assemble_CopyElementData(mesh->Elements,&target,&in);
863 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
864 Assemble_AverageElementData(mesh->Elements,&target,&in);
865 } else {
866 throw DudleyException("No interpolation with data on elements possible.");
867 }
868 break;
869 case(ReducedElements):
870 if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
871 Assemble_CopyElementData(mesh->Elements,&target,&in);
872 } else {
873 throw DudleyException("No interpolation with data on elements with reduced integration order possible.");
874 }
875 break;
876 case(FaceElements):
877 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
878 Assemble_CopyElementData(mesh->FaceElements,&target,&in);
879 } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
880 Assemble_AverageElementData(mesh->FaceElements,&target,&in);
881 } else {
882 throw DudleyException("No interpolation with data on face elements possible.");
883 }
884 break;
885 case(ReducedFaceElements):
886 if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
887 Assemble_CopyElementData(mesh->FaceElements,&target,&in);
888 } else {
889 throw DudleyException("No interpolation with data on face elements with reduced integration order possible.");
890 }
891 break;
892 case(Points):
893 if (target.getFunctionSpace().getTypeCode()==Points) {
894 Assemble_CopyElementData(mesh->Points,&target,&in);
895 } else {
896 throw DudleyException("No interpolation with data on points possible.");
897 }
898 break;
899 case(DegreesOfFreedom):
900 switch(target.getFunctionSpace().getTypeCode()) {
901 case(ReducedDegreesOfFreedom):
902 case(DegreesOfFreedom):
903 Assemble_CopyNodalData(mesh->Nodes,&target,&in);
904 break;
905
906 case(Nodes):
907 case(ReducedNodes):
908 if (getMPISize()>1) {
909 escript::Data temp=escript::Data(in);
910 temp.expand();
911 Assemble_CopyNodalData(mesh->Nodes,&target,&temp);
912 } else {
913 Assemble_CopyNodalData(mesh->Nodes,&target,&in);
914 }
915 break;
916 case(Elements):
917 case(ReducedElements):
918 if (getMPISize()>1) {
919 escript::Data temp=escript::Data( in, continuousFunction(*this) );
920 Assemble_interpolate(mesh->Nodes,mesh->Elements,&temp,&target);
921 } else {
922 Assemble_interpolate(mesh->Nodes,mesh->Elements,&in,&target);
923 }
924 break;
925 case(FaceElements):
926 case(ReducedFaceElements):
927 if (getMPISize()>1) {
928 escript::Data temp=escript::Data( in, continuousFunction(*this) );
929 Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&temp,&target);
930
931 } else {
932 Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&in,&target);
933 }
934 break;
935 case(Points):
936 if (getMPISize()>1) {
937 //escript::Data temp=escript::Data( in, continuousFunction(*this) );
938 //escriptDataC _in2 = temp.getDataC();
939 } else {
940 Assemble_interpolate(mesh->Nodes,mesh->Points,&in,&target);
941 }
942 break;
943 default:
944 stringstream temp;
945 temp << "Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
946 throw DudleyException(temp.str());
947 break;
948 }
949 break;
950 case(ReducedDegreesOfFreedom):
951 switch(target.getFunctionSpace().getTypeCode()) {
952 case(Nodes):
953 throw DudleyException("Dudley does not support interpolation from reduced degrees of freedom to mesh nodes.");
954 break;
955 case(ReducedNodes):
956 if (getMPISize()>1) {
957 escript::Data temp=escript::Data(in);
958 temp.expand();
959 Assemble_CopyNodalData(mesh->Nodes,&target,&temp);
960 } else {
961 Assemble_CopyNodalData(mesh->Nodes,&target,&in);
962 }
963 break;
964 case(DegreesOfFreedom):
965 throw DudleyException("Dudley does not support interpolation from reduced degrees of freedom to degrees of freedom");
966 break;
967 case(ReducedDegreesOfFreedom):
968 Assemble_CopyNodalData(mesh->Nodes,&target,&in);
969 break;
970 case(Elements):
971 case(ReducedElements):
972 if (getMPISize()>1) {
973 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
974 Assemble_interpolate(mesh->Nodes,mesh->Elements,&temp,&target);
975 } else {
976 Assemble_interpolate(mesh->Nodes,mesh->Elements,&in,&target);
977 }
978 break;
979 case(FaceElements):
980 case(ReducedFaceElements):
981 if (getMPISize()>1) {
982 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
983 Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&temp,&target);
984 } else {
985 Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&in,&target);
986 }
987 break;
988 case(Points):
989 if (getMPISize()>1) {
990 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
991 Assemble_interpolate(mesh->Nodes,mesh->Points,&temp,&target);
992 } else {
993 Assemble_interpolate(mesh->Nodes,mesh->Points,&in,&target);
994 }
995 break;
996 default:
997 stringstream temp;
998 temp << "Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
999 throw DudleyException(temp.str());
1000 break;
1001 }
1002 break;
1003 default:
1004 stringstream temp;
1005 temp << "Interpolation On Domain: Dudley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1006 throw DudleyException(temp.str());
1007 break;
1008 }
1009 }
1010
1011 //
1012 // copies the locations of sample points into x:
1013 //
1014 void MeshAdapter::setToX(escript::Data& arg) const
1015 {
1016 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1017 if (argDomain!=*this)
1018 throw DudleyException("Illegal domain of data point locations");
1019 Dudley_Mesh* mesh=m_dudleyMesh.get();
1020 // in case of values node coordinates we can do the job directly:
1021 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1022 Assemble_NodeCoordinates(mesh->Nodes,&arg);
1023 } else {
1024 escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1025 Assemble_NodeCoordinates(mesh->Nodes,&tmp_data);
1026 // this is then interpolated onto arg:
1027 interpolateOnDomain(arg,tmp_data);
1028 }
1029 }
1030
1031 //
1032 // return the normal vectors at the location of data points as a Data object:
1033 //
1034 void MeshAdapter::setToNormal(escript::Data& normal) const
1035 {
1036 /* const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1037 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1038 if (normalDomain!=*this)
1039 throw DudleyException("Illegal domain of normal locations");
1040 Dudley_Mesh* mesh=m_dudleyMesh.get();
1041 switch(normal.getFunctionSpace().getTypeCode()) {
1042 case(Nodes):
1043 throw DudleyException("Dudley does not support surface normal vectors for nodes");
1044 break;
1045 case(ReducedNodes):
1046 throw DudleyException("Dudley does not support surface normal vectors for reduced nodes");
1047 break;
1048 case(Elements):
1049 throw DudleyException("Dudley does not support surface normal vectors for elements");
1050 break;
1051 case(ReducedElements):
1052 throw DudleyException("Dudley does not support surface normal vectors for elements with reduced integration order");
1053 break;
1054 case (FaceElements):
1055 Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&normal);
1056 break;
1057 case (ReducedFaceElements):
1058 Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&normal);
1059 break;
1060 case(Points):
1061 throw DudleyException("Dudley does not support surface normal vectors for point elements");
1062 break;
1063 case(DegreesOfFreedom):
1064 throw DudleyException("Dudley does not support surface normal vectors for degrees of freedom.");
1065 break;
1066 case(ReducedDegreesOfFreedom):
1067 throw DudleyException("Dudley does not support surface normal vectors for reduced degrees of freedom.");
1068 break;
1069 default:
1070 stringstream temp;
1071 temp << "Normal Vectors: Dudley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1072 throw DudleyException(temp.str());
1073 break;
1074 }
1075 }
1076
1077 //
1078 // interpolates data to other domain:
1079 //
1080 void MeshAdapter::interpolateAcross(escript::Data& target,const escript::Data& source) const
1081 {
1082 const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1083 const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1084 if (targetDomain!=this)
1085 throw DudleyException("Illegal domain of interpolation target");
1086
1087 throw DudleyException("Dudley does not allow interpolation across domains yet.");
1088 }
1089
1090 //
1091 // calculates the integral of a function defined of arg:
1092 //
1093 void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1094 {
1095 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1096 if (argDomain!=*this)
1097 throw DudleyException("Illegal domain of integration kernel");
1098
1099 Dudley_Mesh* mesh=m_dudleyMesh.get();
1100 escript::Data temp;
1101 switch(arg.getFunctionSpace().getTypeCode()) {
1102 case(Nodes):
1103 temp=escript::Data( arg, escript::function(*this) );
1104 Assemble_integrate(mesh->Nodes,mesh->Elements,&temp,&integrals[0]);
1105 break;
1106 case(ReducedNodes):
1107 temp=escript::Data( arg, escript::function(*this) );
1108 Assemble_integrate(mesh->Nodes,mesh->Elements,&temp,&integrals[0]);
1109 break;
1110 case(Elements):
1111 Assemble_integrate(mesh->Nodes,mesh->Elements,&arg,&integrals[0]);
1112 break;
1113 case(ReducedElements):
1114 Assemble_integrate(mesh->Nodes,mesh->Elements,&arg,&integrals[0]);
1115 break;
1116 case(FaceElements):
1117 Assemble_integrate(mesh->Nodes,mesh->FaceElements,&arg,&integrals[0]);
1118 break;
1119 case(ReducedFaceElements):
1120 Assemble_integrate(mesh->Nodes,mesh->FaceElements,&arg,&integrals[0]);
1121 break;
1122 case(Points):
1123 throw DudleyException("Integral of data on points is not supported.");
1124 break;
1125 case(DegreesOfFreedom):
1126 temp=escript::Data( arg, escript::function(*this) );
1127 Assemble_integrate(mesh->Nodes,mesh->Elements,&temp,&integrals[0]);
1128 break;
1129 case(ReducedDegreesOfFreedom):
1130 temp=escript::Data( arg, escript::function(*this) );
1131 Assemble_integrate(mesh->Nodes,mesh->Elements,&temp,&integrals[0]);
1132 break;
1133 default:
1134 stringstream temp;
1135 temp << "Integrals: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1136 throw DudleyException(temp.str());
1137 break;
1138 }
1139 }
1140
1141 //
1142 // calculates the gradient of arg:
1143 //
1144 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1145 {
1146 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1147 if (argDomain!=*this)
1148 throw DudleyException("Illegal domain of gradient argument");
1149 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1150 if (gradDomain!=*this)
1151 throw DudleyException("Illegal domain of gradient");
1152
1153 Dudley_Mesh* mesh=m_dudleyMesh.get();
1154 const escript::Data* nodeData=0;
1155 escript::Data temp;
1156 if (getMPISize()>1) {
1157 if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1158 temp=escript::Data( arg, continuousFunction(*this) );
1159 nodeData = &temp;
1160 } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1161 temp=escript::Data( arg, reducedContinuousFunction(*this) );
1162 nodeData = &temp;
1163 } else {
1164 nodeData = &arg;
1165 }
1166 } else {
1167 nodeData = &arg;
1168 }
1169 switch(grad.getFunctionSpace().getTypeCode()) {
1170 case(Nodes):
1171 throw DudleyException("Gradient at nodes is not supported.");
1172 break;
1173 case(ReducedNodes):
1174 throw DudleyException("Gradient at reduced nodes is not supported.");
1175 break;
1176 case(Elements):
1177 Assemble_gradient(mesh->Nodes,mesh->Elements,&grad, nodeData);
1178 break;
1179 case(ReducedElements):
1180 Assemble_gradient(mesh->Nodes,mesh->Elements,&grad, nodeData);
1181 break;
1182 case(FaceElements):
1183 Assemble_gradient(mesh->Nodes,mesh->FaceElements,&grad, nodeData);
1184 break;
1185 case(ReducedFaceElements):
1186 Assemble_gradient(mesh->Nodes,mesh->FaceElements,&grad, nodeData);
1187 break;
1188 case(Points):
1189 throw DudleyException("Gradient at points is not supported.");
1190 break;
1191 case(DegreesOfFreedom):
1192 throw DudleyException("Gradient at degrees of freedom is not supported.");
1193 break;
1194 case(ReducedDegreesOfFreedom):
1195 throw DudleyException("Gradient at reduced degrees of freedom is not supported.");
1196 break;
1197 default:
1198 stringstream temp;
1199 temp << "Gradient: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1200 throw DudleyException(temp.str());
1201 break;
1202 }
1203 }
1204
1205 //
1206 // returns the size of elements:
1207 //
1208 void MeshAdapter::setToSize(escript::Data& size) const
1209 {
1210 Dudley_Mesh* mesh=m_dudleyMesh.get();
1211 switch(size.getFunctionSpace().getTypeCode()) {
1212 case(Nodes):
1213 throw DudleyException("Size of nodes is not supported.");
1214 break;
1215 case(ReducedNodes):
1216 throw DudleyException("Size of reduced nodes is not supported.");
1217 break;
1218 case(Elements):
1219 Assemble_getSize(mesh->Nodes,mesh->Elements,&size);
1220 break;
1221 case(ReducedElements):
1222 Assemble_getSize(mesh->Nodes,mesh->Elements,&size);
1223 break;
1224 case(FaceElements):
1225 Assemble_getSize(mesh->Nodes,mesh->FaceElements,&size);
1226 break;
1227 case(ReducedFaceElements):
1228 Assemble_getSize(mesh->Nodes,mesh->FaceElements,&size);
1229 break;
1230 case(Points):
1231 throw DudleyException("Size of point elements is not supported.");
1232 break;
1233 case(DegreesOfFreedom):
1234 throw DudleyException("Size of degrees of freedom is not supported.");
1235 break;
1236 case(ReducedDegreesOfFreedom):
1237 throw DudleyException("Size of reduced degrees of freedom is not supported.");
1238 break;
1239 default:
1240 stringstream temp;
1241 temp << "Element size: Dudley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1242 throw DudleyException(temp.str());
1243 break;
1244 }
1245 }
1246
1247 //
1248 // sets the location of nodes
1249 //
1250 void MeshAdapter::setNewX(const escript::Data& new_x)
1251 {
1252 Dudley_Mesh* mesh=m_dudleyMesh.get();
1253 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1254 if (newDomain!=*this)
1255 throw DudleyException("Illegal domain of new point locations");
1256 if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1257 Dudley_Mesh_setCoordinates(mesh,&new_x);
1258 } else {
1259 throw DudleyException("As of version escript3.3 - SetNewX only accepts ContinuousFunction arguments please interpolate.");
1260
1261 }
1262 }
1263
1264 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1265 {
1266 if (getMPISize()>1) {
1267 #ifdef ESYS_MPI
1268 index_t myFirstNode=0, myLastNode=0, k=0;
1269 index_t* globalNodeIndex=0;
1270 Dudley_Mesh* mesh_p=m_dudleyMesh.get();
1271 if (fs_code == DUDLEY_REDUCED_NODES)
1272 {
1273 myFirstNode = Dudley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1274 myLastNode = Dudley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1275 globalNodeIndex = Dudley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1276 }
1277 else if (fs_code == DUDLEY_NODES)
1278 {
1279 myFirstNode = Dudley_NodeFile_getFirstNode(mesh_p->Nodes);
1280 myLastNode = Dudley_NodeFile_getLastNode(mesh_p->Nodes);
1281 globalNodeIndex = Dudley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1282 }
1283 else
1284 {
1285 throw DudleyException("unsupported function space type for ownSample()");
1286 }
1287 k=globalNodeIndex[id];
1288 return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1289 #endif
1290 }
1291 return true;
1292 }
1293
1294
1295
1296 //
1297 // creates a stiffness matrix an initializes it with zeros
1298 //
1299 ASM_ptr MeshAdapter::newSystemMatrix(const int row_blocksize,
1300 const escript::FunctionSpace& row_functionspace,
1301 const int column_blocksize,
1302 const escript::FunctionSpace& column_functionspace,
1303 const int type) const
1304 {
1305 int reduceRowOrder=0;
1306 int reduceColOrder=0;
1307 // is the domain right?
1308 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1309 if (row_domain!=*this)
1310 throw DudleyException("domain of row function space does not match the domain of matrix generator.");
1311 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1312 if (col_domain!=*this)
1313 throw DudleyException("domain of column function space does not match the domain of matrix generator.");
1314 // is the function space type right
1315 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1316 reduceRowOrder=0;
1317 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1318 reduceRowOrder=1;
1319 } else {
1320 throw DudleyException("illegal function space type for system matrix rows.");
1321 }
1322 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1323 reduceColOrder=0;
1324 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1325 reduceColOrder=1;
1326 } else {
1327 throw DudleyException("illegal function space type for system matrix columns.");
1328 }
1329 // generate matrix:
1330
1331 paso::SystemMatrixPattern_ptr fsystemMatrixPattern(Dudley_getPattern(getDudley_Mesh(),reduceRowOrder,reduceColOrder));
1332 paso::SystemMatrix_ptr sm;
1333 int trilinos = 0;
1334 if (trilinos) {
1335 #ifdef TRILINOS
1336 /* Allocation Epetra_VrbMatrix here */
1337 #endif
1338 }
1339 else {
1340 sm.reset(new paso::SystemMatrix(type, fsystemMatrixPattern,
1341 row_blocksize, column_blocksize, false, row_functionspace,
1342 column_functionspace));
1343 }
1344 return sm;
1345 }
1346
1347 //
1348 // creates a TransportProblem
1349 //
1350 ATP_ptr MeshAdapter::newTransportProblem(int blocksize,
1351 const escript::FunctionSpace& fs,
1352 int type) const
1353 {
1354 int reduceOrder=0;
1355 // is the domain right?
1356 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(fs.getDomain()));
1357 if (domain!=*this)
1358 throw DudleyException("domain of function space does not match the domain of transport problem generator.");
1359 // is the function space type right
1360 if (fs.getTypeCode()==DegreesOfFreedom) {
1361 reduceOrder=0;
1362 } else if (fs.getTypeCode()==ReducedDegreesOfFreedom) {
1363 reduceOrder=1;
1364 } else {
1365 throw DudleyException("illegal function space type for system matrix rows.");
1366 }
1367 // generate matrix:
1368
1369 paso::SystemMatrixPattern_ptr fsystemMatrixPattern(Dudley_getPattern(
1370 getDudley_Mesh(),reduceOrder,reduceOrder));
1371 paso::TransportProblem_ptr transportProblem(new paso::TransportProblem(
1372 fsystemMatrixPattern, blocksize,
1373 fs));
1374 return transportProblem;
1375 }
1376
1377 //
1378 // vtkObject MeshAdapter::createVtkObject() const
1379 // TODO:
1380 //
1381
1382 //
1383 // returns true if data at the atom_type is considered as being cell centered:
1384 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1385 {
1386 switch(functionSpaceCode) {
1387 case(Nodes):
1388 case(DegreesOfFreedom):
1389 case(ReducedDegreesOfFreedom):
1390 return false;
1391 break;
1392 case(Elements):
1393 case(FaceElements):
1394 case(Points):
1395 case(ReducedElements):
1396 case(ReducedFaceElements):
1397 return true;
1398 break;
1399 default:
1400 stringstream temp;
1401 temp << "Cell: Dudley does not know anything about function space type " << functionSpaceCode;
1402 throw DudleyException(temp.str());
1403 break;
1404 }
1405 return false;
1406 }
1407
1408 bool
1409 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1410 {
1411 /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1412 class 1: DOF <-> Nodes
1413 class 2: ReducedDOF <-> ReducedNodes
1414 class 3: Points
1415 class 4: Elements
1416 class 5: ReducedElements
1417 class 6: FaceElements
1418 class 7: ReducedFaceElements
1419 class 8: ContactElementZero <-> ContactElementOne
1420 class 9: ReducedContactElementZero <-> ReducedContactElementOne
1421
1422 There is also a set of lines. Interpolation is possible down a line but not between lines.
1423 class 1 and 2 belong to all lines so aren't considered.
1424 line 0: class 3
1425 line 1: class 4,5
1426 line 2: class 6,7
1427 line 3: class 8,9
1428
1429 For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1430 eg hasnodes is true if we have at least one instance of Nodes.
1431 */
1432 if (fs.empty())
1433 {
1434 return false;
1435 }
1436 vector<int> hasclass(10);
1437 vector<int> hasline(4);
1438 bool hasnodes=false;
1439 bool hasrednodes=false;
1440 for (int i=0;i<fs.size();++i)
1441 {
1442 switch(fs[i])
1443 {
1444 case(Nodes): hasnodes=true; // no break is deliberate
1445 case(DegreesOfFreedom):
1446 hasclass[1]=1;
1447 break;
1448 case(ReducedNodes): hasrednodes=true; // no break is deliberate
1449 case(ReducedDegreesOfFreedom):
1450 hasclass[2]=1;
1451 break;
1452 case(Points):
1453 hasline[0]=1;
1454 hasclass[3]=1;
1455 break;
1456 case(Elements):
1457 hasclass[4]=1;
1458 hasline[1]=1;
1459 break;
1460 case(ReducedElements):
1461 hasclass[5]=1;
1462 hasline[1]=1;
1463 break;
1464 case(FaceElements):
1465 hasclass[6]=1;
1466 hasline[2]=1;
1467 break;
1468 case(ReducedFaceElements):
1469 hasclass[7]=1;
1470 hasline[2]=1;
1471 break;
1472 default:
1473 return false;
1474 }
1475 }
1476 int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1477 // fail if we have more than one leaf group
1478
1479 if (totlines>1)
1480 {
1481 return false; // there are at least two branches we can't interpolate between
1482 }
1483 else if (totlines==1)
1484 {
1485 if (hasline[0]==1) // we have points
1486 {
1487 resultcode=Points;
1488 }
1489 else if (hasline[1]==1)
1490 {
1491 if (hasclass[5]==1)
1492 {
1493 resultcode=ReducedElements;
1494 }
1495 else
1496 {
1497 resultcode=Elements;
1498 }
1499 }
1500 else if (hasline[2]==1)
1501 {
1502 if (hasclass[7]==1)
1503 {
1504 resultcode=ReducedFaceElements;
1505 }
1506 else
1507 {
1508 resultcode=FaceElements;
1509 }
1510 }
1511 else // so we must be in line3
1512 {
1513
1514 throw DudleyException("choosing between contact elements - we should never get here.");
1515
1516 }
1517 }
1518 else // totlines==0
1519 {
1520 if (hasclass[2]==1)
1521 {
1522 // something from class 2
1523 resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1524 }
1525 else
1526 { // something from class 1
1527 resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1528 }
1529 }
1530 return true;
1531 }
1532
1533 signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1534 {
1535 if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1536 {
1537 return 1;
1538 }
1539 else if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1540 {
1541 return -1;
1542 }
1543 return 0;
1544 }
1545
1546
1547
1548 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1549 {
1550 switch(functionSpaceType_source) {
1551 case(Nodes):
1552 switch(functionSpaceType_target) {
1553 case(Nodes):
1554 case(ReducedNodes):
1555 case(ReducedDegreesOfFreedom):
1556 case(DegreesOfFreedom):
1557 case(Elements):
1558 case(ReducedElements):
1559 case(FaceElements):
1560 case(ReducedFaceElements):
1561 case(Points):
1562 return true;
1563 default:
1564 stringstream temp;
1565 temp << "Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1566 throw DudleyException(temp.str());
1567 }
1568 break;
1569 case(ReducedNodes):
1570 switch(functionSpaceType_target) {
1571 case(ReducedNodes):
1572 case(ReducedDegreesOfFreedom):
1573 case(Elements):
1574 case(ReducedElements):
1575 case(FaceElements):
1576 case(ReducedFaceElements):
1577 case(Points):
1578 return true;
1579 case(Nodes):
1580 case(DegreesOfFreedom):
1581 return false;
1582 default:
1583 stringstream temp;
1584 temp << "Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1585 throw DudleyException(temp.str());
1586 }
1587 break;
1588 case(Elements):
1589 if (functionSpaceType_target==Elements) {
1590 return true;
1591 } else if (functionSpaceType_target==ReducedElements) {
1592 return true;
1593 } else {
1594 return false;
1595 }
1596 case(ReducedElements):
1597 if (functionSpaceType_target==ReducedElements) {
1598 return true;
1599 } else {
1600 return false;
1601 }
1602 case(FaceElements):
1603 if (functionSpaceType_target==FaceElements) {
1604 return true;
1605 } else if (functionSpaceType_target==ReducedFaceElements) {
1606 return true;
1607 } else {
1608 return false;
1609 }
1610 case(ReducedFaceElements):
1611 if (functionSpaceType_target==ReducedFaceElements) {
1612 return true;
1613 } else {
1614 return false;
1615 }
1616 case(Points):
1617 if (functionSpaceType_target==Points) {
1618 return true;
1619 } else {
1620 return false;
1621 }
1622 case(DegreesOfFreedom):
1623 switch(functionSpaceType_target) {
1624 case(ReducedDegreesOfFreedom):
1625 case(DegreesOfFreedom):
1626 case(Nodes):
1627 case(ReducedNodes):
1628 case(Elements):
1629 case(ReducedElements):
1630 case(Points):
1631 case(FaceElements):
1632 case(ReducedFaceElements):
1633 return true;
1634 default:
1635 stringstream temp;
1636 temp << "Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1637 throw DudleyException(temp.str());
1638 }
1639 break;
1640 case(ReducedDegreesOfFreedom):
1641 switch(functionSpaceType_target) {
1642 case(ReducedDegreesOfFreedom):
1643 case(ReducedNodes):
1644 case(Elements):
1645 case(ReducedElements):
1646 case(FaceElements):
1647 case(ReducedFaceElements):
1648 case(Points):
1649 return true;
1650 case(Nodes):
1651 case(DegreesOfFreedom):
1652 return false;
1653 default:
1654 stringstream temp;
1655 temp << "Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1656 throw DudleyException(temp.str());
1657 }
1658 break;
1659 default:
1660 stringstream temp;
1661 temp << "Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_source;
1662 throw DudleyException(temp.str());
1663 break;
1664 }
1665 return false;
1666 }
1667
1668 bool MeshAdapter::probeInterpolationAcross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1669 {
1670 return false;
1671 }
1672
1673 bool MeshAdapter::operator==(const AbstractDomain& other) const
1674 {
1675 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1676 if (temp!=0) {
1677 return (m_dudleyMesh==temp->m_dudleyMesh);
1678 } else {
1679 return false;
1680 }
1681 }
1682
1683 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1684 {
1685 return !(operator==(other));
1686 }
1687
1688 int MeshAdapter::getSystemMatrixTypeId(const bp::object& options) const
1689 {
1690 const escript::SolverBuddy& sb = bp::extract<escript::SolverBuddy>(options);
1691
1692 return paso::SystemMatrix::getSystemMatrixTypeId(sb.getSolverMethod(),
1693 sb.getPreconditioner(), sb.getPackage(), sb.isSymmetric(),
1694 m_dudleyMesh->MPIInfo);
1695 }
1696
1697 int MeshAdapter::getTransportTypeId(int solver, int preconditioner,
1698 int package, bool symmetry) const
1699 {
1700 Dudley_Mesh* mesh=m_dudleyMesh.get();
1701 return paso::TransportProblem::getTypeId(solver, preconditioner, package,
1702 symmetry, mesh->MPIInfo);
1703 }
1704
1705 escript::Data MeshAdapter::getX() const
1706 {
1707 return continuousFunction(*this).getX();
1708 }
1709
1710 escript::Data MeshAdapter::getNormal() const
1711 {
1712 return functionOnBoundary(*this).getNormal();
1713 }
1714
1715 escript::Data MeshAdapter::getSize() const
1716 {
1717 return escript::function(*this).getSize();
1718 }
1719
1720 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1721 {
1722 int *out = NULL;
1723 Dudley_Mesh* mesh=m_dudleyMesh.get();
1724 switch (functionSpaceType) {
1725 case(Nodes):
1726 out=mesh->Nodes->Id;
1727 break;
1728 case(ReducedNodes):
1729 out=mesh->Nodes->reducedNodesId;
1730 break;
1731 case(Elements):
1732 out=mesh->Elements->Id;
1733 break;
1734 case(ReducedElements):
1735 out=mesh->Elements->Id;
1736 break;
1737 case(FaceElements):
1738 out=mesh->FaceElements->Id;
1739 break;
1740 case(ReducedFaceElements):
1741 out=mesh->FaceElements->Id;
1742 break;
1743 case(Points):
1744 out=mesh->Points->Id;
1745 break;
1746 case(DegreesOfFreedom):
1747 out=mesh->Nodes->degreesOfFreedomId;
1748 break;
1749 case(ReducedDegreesOfFreedom):
1750 out=mesh->Nodes->reducedDegreesOfFreedomId;
1751 break;
1752 default:
1753 stringstream temp;
1754 temp << "Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1755 throw DudleyException(temp.str());
1756 break;
1757 }
1758 return out;
1759 }
1760 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1761 {
1762 int out=0;
1763 Dudley_Mesh* mesh=m_dudleyMesh.get();
1764 switch (functionSpaceType) {
1765 case(Nodes):
1766 out=mesh->Nodes->Tag[sampleNo];
1767 break;
1768 case(ReducedNodes):
1769 throw DudleyException("ReducedNodes does not support tags.");
1770 break;
1771 case(Elements):
1772 out=mesh->Elements->Tag[sampleNo];
1773 break;
1774 case(ReducedElements):
1775 out=mesh->Elements->Tag[sampleNo];
1776 break;
1777 case(FaceElements):
1778 out=mesh->FaceElements->Tag[sampleNo];
1779 break;
1780 case(ReducedFaceElements):
1781 out=mesh->FaceElements->Tag[sampleNo];
1782 break;
1783 case(Points):
1784 out=mesh->Points->Tag[sampleNo];
1785 break;
1786 case(DegreesOfFreedom):
1787 throw DudleyException("DegreesOfFreedom does not support tags.");
1788 break;
1789 case(ReducedDegreesOfFreedom):
1790 throw DudleyException("ReducedDegreesOfFreedom does not support tags.");
1791 break;
1792 default:
1793 stringstream temp;
1794 temp << "Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1795 throw DudleyException(temp.str());
1796 break;
1797 }
1798 return out;
1799 }
1800
1801
1802 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1803 {
1804 Dudley_Mesh* mesh=m_dudleyMesh.get();
1805 switch(functionSpaceType) {
1806 case(Nodes):
1807 Dudley_NodeFile_setTags(mesh->Nodes,newTag,&mask);
1808 break;
1809 case(ReducedNodes):
1810 throw DudleyException("ReducedNodes does not support tags");
1811 break;
1812 case(DegreesOfFreedom):
1813 throw DudleyException("DegreesOfFreedom does not support tags");
1814 break;
1815 case(ReducedDegreesOfFreedom):
1816 throw DudleyException("ReducedDegreesOfFreedom does not support tags");
1817 break;
1818 case(Elements):
1819 Dudley_ElementFile_setTags(mesh->Elements,newTag,&mask);
1820 break;
1821 case(ReducedElements):
1822 Dudley_ElementFile_setTags(mesh->Elements,newTag,&mask);
1823 break;
1824 case(FaceElements):
1825 Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&mask);
1826 break;
1827 case(ReducedFaceElements):
1828 Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&mask);
1829 break;
1830 case(Points):
1831 Dudley_ElementFile_setTags(mesh->Points,newTag,&mask);
1832 break;
1833 default:
1834 stringstream temp;
1835 temp << "Dudley does not know anything about function space type " << functionSpaceType;
1836 throw DudleyException(temp.str());
1837 }
1838 }
1839
1840 void MeshAdapter::setTagMap(const string& name, int tag)
1841 {
1842 Dudley_Mesh* mesh=m_dudleyMesh.get();
1843 Dudley_Mesh_addTagMap(mesh, name.c_str(),tag);
1844 }
1845
1846 int MeshAdapter::getTag(const string& name) const
1847 {
1848 Dudley_Mesh* mesh=m_dudleyMesh.get();
1849 int tag=0;
1850 tag=Dudley_Mesh_getTag(mesh, name.c_str());
1851 return tag;
1852 }
1853
1854 bool MeshAdapter::isValidTagName(const string& name) const
1855 {
1856 Dudley_Mesh* mesh=m_dudleyMesh.get();
1857 return Dudley_Mesh_isValidTagName(mesh,name.c_str());
1858 }
1859
1860 string MeshAdapter::showTagNames() const
1861 {
1862 stringstream temp;
1863 Dudley_Mesh* mesh=m_dudleyMesh.get();
1864 Dudley_TagMap* tag_map=mesh->TagMap;
1865 while (tag_map) {
1866 temp << tag_map->name;
1867 tag_map=tag_map->next;
1868 if (tag_map) temp << ", ";
1869 }
1870 return temp.str();
1871 }
1872
1873 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
1874 {
1875 Dudley_Mesh* mesh=m_dudleyMesh.get();
1876 dim_t numTags=0;
1877 switch(functionSpaceCode) {
1878 case(Nodes):
1879 numTags=mesh->Nodes->numTagsInUse;
1880 break;
1881 case(ReducedNodes):
1882 throw DudleyException("ReducedNodes does not support tags");
1883 break;
1884 case(DegreesOfFreedom):
1885 throw DudleyException("DegreesOfFreedom does not support tags");
1886 break;
1887 case(ReducedDegreesOfFreedom):
1888 throw DudleyException("ReducedDegreesOfFreedom does not support tags");
1889 break;
1890 case(Elements):
1891 case(ReducedElements):
1892 numTags=mesh->Elements->numTagsInUse;
1893 break;
1894 case(FaceElements):
1895 case(ReducedFaceElements):
1896 numTags=mesh->FaceElements->numTagsInUse;
1897 break;
1898 case(Points):
1899 numTags=mesh->Points->numTagsInUse;
1900 break;
1901 default:
1902 stringstream temp;
1903 temp << "Dudley does not know anything about function space type " << functionSpaceCode;
1904 throw DudleyException(temp.str());
1905 }
1906 return numTags;
1907 }
1908
1909 const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
1910 {
1911 Dudley_Mesh* mesh=m_dudleyMesh.get();
1912 index_t* tags=NULL;
1913 switch(functionSpaceCode) {
1914 case(Nodes):
1915 tags=mesh->Nodes->tagsInUse;
1916 break;
1917 case(ReducedNodes):
1918 throw DudleyException("ReducedNodes does not support tags");
1919 break;
1920 case(DegreesOfFreedom):
1921 throw DudleyException("DegreesOfFreedom does not support tags");
1922 break;
1923 case(ReducedDegreesOfFreedom):
1924 throw DudleyException("ReducedDegreesOfFreedom does not support tags");
1925 break;
1926 case(Elements):
1927 case(ReducedElements):
1928 tags=mesh->Elements->tagsInUse;
1929 break;
1930 case(FaceElements):
1931 case(ReducedFaceElements):
1932 tags=mesh->FaceElements->tagsInUse;
1933 break;
1934 case(Points):
1935 tags=mesh->Points->tagsInUse;
1936 break;
1937 default:
1938 stringstream temp;
1939 temp << "Dudley does not know anything about function space type " << functionSpaceCode;
1940 throw DudleyException(temp.str());
1941 }
1942 return tags;
1943 }
1944
1945
1946 bool MeshAdapter::canTag(int functionSpaceCode) const
1947 {
1948 switch(functionSpaceCode) {
1949 case(Nodes):
1950 case(Elements):
1951 case(ReducedElements):
1952 case(FaceElements):
1953 case(ReducedFaceElements):
1954 case(Points):
1955 return true;
1956 case(ReducedNodes):
1957 case(DegreesOfFreedom):
1958 case(ReducedDegreesOfFreedom):
1959 return false;
1960 default:
1961 return false;
1962 }
1963 }
1964
1965 AbstractDomain::StatusType MeshAdapter::getStatus() const
1966 {
1967 Dudley_Mesh* mesh=m_dudleyMesh.get();
1968 return Dudley_Mesh_getStatus(mesh);
1969 }
1970
1971 int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
1972 {
1973
1974 Dudley_Mesh* mesh=m_dudleyMesh.get();
1975 int order =-1;
1976 switch(functionSpaceCode) {
1977 case(Nodes):
1978 case(DegreesOfFreedom):
1979 order=mesh->approximationOrder;
1980 break;
1981 case(ReducedNodes):
1982 case(ReducedDegreesOfFreedom):
1983 order=mesh->reducedApproximationOrder;
1984 break;
1985 case(Elements):
1986 case(FaceElements):
1987 case(Points):
1988 order=mesh->integrationOrder;
1989 break;
1990 case(ReducedElements):
1991 case(ReducedFaceElements):
1992 order=mesh->reducedIntegrationOrder;
1993 break;
1994 default:
1995 stringstream temp;
1996 temp << "Dudley does not know anything about function space type " << functionSpaceCode;
1997 throw DudleyException(temp.str());
1998 }
1999 return order;
2000 }
2001
2002
2003 bool MeshAdapter::supportsContactElements() const
2004 {
2005 return false;
2006 }
2007
2008 escript::Data MeshAdapter::randomFill(const escript::DataTypes::ShapeType& shape,
2009 const escript::FunctionSpace& what, long seed,
2010 const bp::tuple& filter) const
2011 {
2012 Data towipe(0, shape, what, true);
2013 // since we just made this object, no sharing is possible and we don't need to check for
2014 // exlusive write
2015 escript::DataTypes::RealVectorType& dv=towipe.getExpandedVectorReference();
2016 const size_t dvsize=dv.size();
2017 escript::randomFillArray(seed, &(dv[0]), dvsize);
2018 return towipe;
2019 }
2020
2021 } // end of namespace
2022

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26