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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5136 - (show annotations)
Tue Sep 9 07:13:55 2014 UTC (4 years, 5 months ago) by caltinay
File size: 72466 byte(s)
ripley now supports paso solvers again and returns an appropriate matrix type
id. Changed the getSystemMatrixTypeId() method to take a full SolverBuddy
instance and made some other simplifications.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26