/[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 6079 - (show annotations)
Mon Mar 21 12:22:38 2016 UTC (2 years, 11 months ago) by caltinay
File size: 59145 byte(s)
Big commit - making dudley much more like finley to make it more
managable. Fixed quite a few issues that had been fixed in finley.
Disposed of all ReducedNode/ReducedDOF entities that dudley never supported.
Compiles and passes tests.

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 #include <dudley/Assemble.h>
19 #include <dudley/DudleyException.h>
20
21 #include <escript/Data.h>
22 #include <escript/DataFactory.h>
23 #include <escript/Random.h>
24 #include <escript/SolverOptions.h>
25
26 #include <paso/SystemMatrix.h>
27 #include <paso/Transport.h>
28
29 #ifdef USE_NETCDF
30 #include <netcdfcpp.h>
31 #endif
32
33 #ifdef USE_TRILINOS
34 #include <trilinoswrap/TrilinosMatrixAdapter.h>
35
36 using esys_trilinos::TrilinosMatrixAdapter;
37 using esys_trilinos::const_TrilinosGraph_ptr;
38 #endif
39
40 using namespace std;
41 namespace bp = boost::python;
42 using escript::ValueError;
43
44 namespace dudley {
45
46 // define the static constants
47 MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
48 // dudley only supports single approximation order
49 const int MeshAdapter::DegreesOfFreedom=DUDLEY_DEGREES_OF_FREEDOM;
50 const int MeshAdapter::Nodes=DUDLEY_NODES;
51 const int MeshAdapter::Elements=DUDLEY_ELEMENTS;
52 const int MeshAdapter::ReducedElements=DUDLEY_REDUCED_ELEMENTS;
53 const int MeshAdapter::FaceElements=DUDLEY_FACE_ELEMENTS;
54 const int MeshAdapter::ReducedFaceElements=DUDLEY_REDUCED_FACE_ELEMENTS;
55 const int MeshAdapter::Points=DUDLEY_POINTS;
56
57 MeshAdapter::MeshAdapter(Mesh* dudleyMesh) :
58 m_dudleyMesh(dudleyMesh)
59 {
60 setFunctionSpaceTypeNames();
61 }
62
63 // The copy constructor should just increment the use count
64 MeshAdapter::MeshAdapter(const MeshAdapter& in) :
65 m_dudleyMesh(in.m_dudleyMesh)
66 {
67 setFunctionSpaceTypeNames();
68 }
69
70 MeshAdapter::~MeshAdapter()
71 {
72 }
73
74 escript::JMPI MeshAdapter::getMPI() const
75 {
76 return m_dudleyMesh->MPIInfo;
77 }
78
79 int MeshAdapter::getMPISize() const
80 {
81 return getMPI()->size;
82 }
83
84 int MeshAdapter::getMPIRank() const
85 {
86 return getMPI()->rank;
87 }
88
89 void MeshAdapter::MPIBarrier() const
90 {
91 #ifdef ESYS_MPI
92 MPI_Barrier(getMPIComm());
93 #endif
94 }
95
96 bool MeshAdapter::onMasterProcessor() const
97 {
98 return getMPIRank() == 0;
99 }
100
101 MPI_Comm MeshAdapter::getMPIComm() const
102 {
103 return getMPI()->comm;
104 }
105
106 Mesh* MeshAdapter::getMesh() const
107 {
108 return m_dudleyMesh.get();
109 }
110
111 void MeshAdapter::write(const string& fileName) const
112 {
113 m_dudleyMesh->write(fileName);
114 }
115
116 void MeshAdapter::Print_Mesh_Info(bool full) const
117 {
118 m_dudleyMesh->printInfo(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 index_t* index_ptr;
127 #ifdef ESYS_INDEXTYPE_LONG
128 NcType ncIdxType = ncLong;
129 #else
130 NcType ncIdxType = ncInt;
131 #endif
132 Mesh* mesh = m_dudleyMesh.get();
133 int num_Tags = 0;
134 int mpi_size = getMPISize();
135 int mpi_rank = getMPIRank();
136 int numDim = mesh->Nodes->numDim;
137 dim_t numNodes = mesh->Nodes->getNumNodes();
138 dim_t num_Elements = mesh->Elements->numElements;
139 dim_t num_FaceElements = mesh->FaceElements->numElements;
140 dim_t num_Points = mesh->Points->numElements;
141 int num_Elements_numNodes = mesh->Elements->numNodes;
142 int num_FaceElements_numNodes = mesh->FaceElements->numNodes;
143 #ifdef ESYS_MPI
144 MPI_Status status;
145 #endif
146
147 // Incoming token indicates it's my turn to write
148 #ifdef ESYS_MPI
149 if (mpi_rank > 0)
150 MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, getMPIComm(), &status);
151 #endif
152
153 const string newFileName(mesh->MPIInfo->appendRankToFileName(fileName));
154
155 // Figure out how much storage is required for tags
156 num_Tags = mesh->tagMap.size();
157
158 // NetCDF error handler
159 NcError err(NcError::verbose_nonfatal);
160 // Create the file
161 NcFile dataFile(newFileName.c_str(), NcFile::Replace);
162 string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
163 // check if writing was successful
164 if (!dataFile.is_valid())
165 throw DudleyException(msgPrefix + "Open file for output");
166
167 // Define dimensions (num_Elements and dim_Elements are identical,
168 // dim_Elements only appears if > 0)
169 if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
170 throw DudleyException(msgPrefix+"add_dim(numNodes)");
171 if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
172 throw DudleyException(msgPrefix+"add_dim(numDim)");
173 if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )
174 throw DudleyException(msgPrefix+"add_dim(mpi_size)");
175 if (num_Elements > 0)
176 if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
177 throw DudleyException(msgPrefix+"add_dim(dim_Elements)");
178 if (num_FaceElements > 0)
179 if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
180 throw DudleyException(msgPrefix+"add_dim(dim_FaceElements)");
181 if (num_Points > 0)
182 if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
183 throw DudleyException(msgPrefix+"add_dim(dim_Points)");
184 if (num_Elements > 0)
185 if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
186 throw DudleyException(msgPrefix+"add_dim(dim_Elements_Nodes)");
187 if (num_FaceElements > 0)
188 if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
189 throw DudleyException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
190 if (num_Tags > 0)
191 if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
192 throw DudleyException(msgPrefix+"add_dim(dim_Tags)");
193
194 // Attributes: MPI size, MPI rank, Name, order, reduced_order
195 if (!dataFile.add_att("index_size", (int)sizeof(index_t)))
196 throw DudleyException(msgPrefix+"add_att(index_size)");
197 if (!dataFile.add_att("mpi_size", mpi_size) )
198 throw DudleyException(msgPrefix+"add_att(mpi_size)");
199 if (!dataFile.add_att("mpi_rank", mpi_rank) )
200 throw DudleyException(msgPrefix+"add_att(mpi_rank)");
201 if (!dataFile.add_att("Name",mesh->m_name.c_str()) )
202 throw DudleyException(msgPrefix+"add_att(Name)");
203 if (!dataFile.add_att("numDim",numDim) )
204 throw DudleyException(msgPrefix+"add_att(order)");
205 if (!dataFile.add_att("order",mesh->integrationOrder) )
206 throw DudleyException(msgPrefix+"add_att(order)");
207 if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
208 throw DudleyException(msgPrefix+"add_att(reduced_order)");
209 if (!dataFile.add_att("numNodes",numNodes) )
210 throw DudleyException(msgPrefix+"add_att(numNodes)");
211 if (!dataFile.add_att("num_Elements",num_Elements) )
212 throw DudleyException(msgPrefix+"add_att(num_Elements)");
213 if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
214 throw DudleyException(msgPrefix+"add_att(num_FaceElements)");
215 if (!dataFile.add_att("num_Points",num_Points) )
216 throw DudleyException(msgPrefix+"add_att(num_Points)");
217 if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
218 throw DudleyException(msgPrefix+"add_att(num_Elements_numNodes)");
219 if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
220 throw DudleyException(msgPrefix+"add_att(num_FaceElements_numNodes)");
221 if (!dataFile.add_att("Elements_TypeId", mesh->Elements->etype) )
222 throw DudleyException(msgPrefix+"add_att(Elements_TypeId)");
223 if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->etype) )
224 throw DudleyException(msgPrefix+"add_att(FaceElements_TypeId)");
225 if (!dataFile.add_att("Points_TypeId", mesh->Points->etype) )
226 throw DudleyException(msgPrefix+"add_att(Points_TypeId)");
227 if (!dataFile.add_att("num_Tags", num_Tags) )
228 throw DudleyException(msgPrefix+"add_att(num_Tags)");
229
230 // // // // // Nodes // // // // //
231
232 // Nodes nodeDistribution
233 if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncIdxType, ncdims[2])) )
234 throw DudleyException(msgPrefix+"add_var(Nodes_NodeDistribution)");
235 index_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
236 if (! (ids->put(index_ptr, mpi_size+1)) )
237 throw DudleyException(msgPrefix+"put(Nodes_NodeDistribution)");
238
239 // Nodes degreesOfFreedomDistribution
240 if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncIdxType, ncdims[2])) )
241 throw DudleyException(msgPrefix+"add_var(Nodes_DofDistribution)");
242 index_ptr = &mesh->Nodes->dofDistribution->first_component[0];
243 if (! (ids->put(index_ptr, mpi_size+1)) )
244 throw DudleyException(msgPrefix+"put(Nodes_DofDistribution)");
245
246 // Only write nodes if non-empty because NetCDF doesn't like empty arrays
247 // (it treats them as NC_UNLIMITED)
248 if (numNodes > 0) {
249 // Nodes Id
250 if (! ( ids = dataFile.add_var("Nodes_Id", ncIdxType, ncdims[0])) )
251 throw DudleyException(msgPrefix+"add_var(Nodes_Id)");
252 if (! (ids->put(mesh->Nodes->Id, numNodes)) )
253 throw DudleyException(msgPrefix+"put(Nodes_Id)");
254
255 // Nodes Tag
256 if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
257 throw DudleyException(msgPrefix+"add_var(Nodes_Tag)");
258 if (! (ids->put(mesh->Nodes->Tag, numNodes)) )
259 throw DudleyException(msgPrefix+"put(Nodes_Tag)");
260
261 // Nodes gDOF
262 if (! ( ids = dataFile.add_var("Nodes_gDOF", ncIdxType, ncdims[0])) )
263 throw DudleyException(msgPrefix+"add_var(Nodes_gDOF)");
264 if (! (ids->put(mesh->Nodes->globalDegreesOfFreedom, numNodes)) )
265 throw DudleyException(msgPrefix+"put(Nodes_gDOF)");
266
267 // Nodes global node index
268 if (! ( ids = dataFile.add_var("Nodes_gNI", ncIdxType, ncdims[0])) )
269 throw DudleyException(msgPrefix+"add_var(Nodes_gNI)");
270 if (! (ids->put(mesh->Nodes->globalNodesIndex, numNodes)) )
271 throw DudleyException(msgPrefix+"put(Nodes_gNI)");
272
273 // Nodes Coordinates
274 if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
275 throw DudleyException(msgPrefix+"add_var(Nodes_Coordinates)");
276 if (! (ids->put(mesh->Nodes->Coordinates, numNodes, numDim)) )
277 throw DudleyException(msgPrefix+"put(Nodes_Coordinates)");
278 }
279
280 // // // // // Elements // // // // //
281 if (num_Elements > 0) {
282 // Elements_Id
283 if (! ( ids = dataFile.add_var("Elements_Id", ncIdxType, ncdims[3])) )
284 throw DudleyException(msgPrefix+"add_var(Elements_Id)");
285 if (! (ids->put(mesh->Elements->Id, num_Elements)) )
286 throw DudleyException(msgPrefix+"put(Elements_Id)");
287
288 // Elements_Tag
289 if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
290 throw DudleyException(msgPrefix+"add_var(Elements_Tag)");
291 if (! (ids->put(mesh->Elements->Tag, num_Elements)) )
292 throw DudleyException(msgPrefix+"put(Elements_Tag)");
293
294 // Elements_Owner
295 if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
296 throw DudleyException(msgPrefix+"add_var(Elements_Owner)");
297 if (! (ids->put(mesh->Elements->Owner, num_Elements)) )
298 throw DudleyException(msgPrefix+"put(Elements_Owner)");
299
300 // Elements_Color
301 if (! ( ids = dataFile.add_var("Elements_Color", ncIdxType, ncdims[3])) )
302 throw DudleyException(msgPrefix+"add_var(Elements_Color)");
303 if (! (ids->put(mesh->Elements->Color, num_Elements)) )
304 throw DudleyException(msgPrefix+"put(Elements_Color)");
305
306 // Elements_Nodes
307 if (! ( ids = dataFile.add_var("Elements_Nodes", ncIdxType, ncdims[3], ncdims[7]) ) )
308 throw DudleyException(msgPrefix+"add_var(Elements_Nodes)");
309 if (! (ids->put(mesh->Elements->Nodes, num_Elements, num_Elements_numNodes)) )
310 throw DudleyException(msgPrefix+"put(Elements_Nodes)");
311 }
312
313 // // // // // Face_Elements // // // // //
314 if (num_FaceElements > 0) {
315 // FaceElements_Id
316 if (!(ids = dataFile.add_var("FaceElements_Id", ncIdxType, ncdims[4])))
317 throw DudleyException(msgPrefix+"add_var(FaceElements_Id)");
318 if (!(ids->put(mesh->FaceElements->Id, num_FaceElements)))
319 throw DudleyException(msgPrefix+"put(FaceElements_Id)");
320
321 // FaceElements_Tag
322 if (!(ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])))
323 throw DudleyException(msgPrefix+"add_var(FaceElements_Tag)");
324 if (!(ids->put(mesh->FaceElements->Tag, num_FaceElements)))
325 throw DudleyException(msgPrefix+"put(FaceElements_Tag)");
326
327 // FaceElements_Owner
328 if (!(ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])))
329 throw DudleyException(msgPrefix+"add_var(FaceElements_Owner)");
330 if (!(ids->put(mesh->FaceElements->Owner, num_FaceElements)))
331 throw DudleyException(msgPrefix+"put(FaceElements_Owner)");
332
333 // FaceElements_Color
334 if (!(ids = dataFile.add_var("FaceElements_Color", ncIdxType, ncdims[4])))
335 throw DudleyException(msgPrefix+"add_var(FaceElements_Color)");
336 if (!(ids->put(mesh->FaceElements->Color, num_FaceElements)))
337 throw DudleyException(msgPrefix+"put(FaceElements_Color)");
338
339 // FaceElements_Nodes
340 if (!(ids = dataFile.add_var("FaceElements_Nodes", ncIdxType, ncdims[4], ncdims[8])))
341 throw DudleyException(msgPrefix+"add_var(FaceElements_Nodes)");
342 if (!(ids->put(mesh->FaceElements->Nodes, num_FaceElements, num_FaceElements_numNodes)))
343 throw DudleyException(msgPrefix+"put(FaceElements_Nodes)");
344 }
345
346 // // // // // Points // // // // //
347 if (num_Points > 0) {
348 // Points_Id
349 if (!(ids = dataFile.add_var("Points_Id", ncIdxType, ncdims[6])))
350 throw DudleyException(msgPrefix+"add_var(Points_Id)");
351 if (!(ids->put(mesh->Points->Id, num_Points)))
352 throw DudleyException(msgPrefix+"put(Points_Id)");
353
354 // Points_Tag
355 if (!(ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])))
356 throw DudleyException(msgPrefix+"add_var(Points_Tag)");
357 if (!(ids->put(mesh->Points->Tag, num_Points)))
358 throw DudleyException(msgPrefix+"put(Points_Tag)");
359
360 // Points_Owner
361 if (!(ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])))
362 throw DudleyException(msgPrefix+"add_var(Points_Owner)");
363 if (!(ids->put(mesh->Points->Owner, num_Points)))
364 throw DudleyException(msgPrefix+"put(Points_Owner)");
365
366 // Points_Color
367 if (!(ids = dataFile.add_var("Points_Color", ncIdxType, ncdims[6])))
368 throw DudleyException(msgPrefix+"add_var(Points_Color)");
369 if (!(ids->put(mesh->Points->Color, num_Points)))
370 throw DudleyException(msgPrefix+"put(Points_Color)");
371
372 // Points_Nodes
373 if (!(ids = dataFile.add_var("Points_Nodes", ncIdxType, ncdims[6])))
374 throw DudleyException(msgPrefix+"add_var(Points_Nodes)");
375 if (!(ids->put(mesh->Points->Nodes, num_Points)))
376 throw DudleyException(msgPrefix+"put(Points_Nodes)");
377 }
378
379 // // // // // TagMap // // // // //
380 if (num_Tags > 0) {
381 // Temp storage to gather node IDs
382 vector<int> Tags_keys;
383
384 // Copy tag data into temp arrays
385 TagMap::const_iterator it;
386 for (it = mesh->tagMap.begin(); it != mesh->tagMap.end(); it++) {
387 Tags_keys.push_back(it->second);
388 }
389
390 // Tags_keys
391 if (!(ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])))
392 throw DudleyException(msgPrefix+"add_var(Tags_keys)");
393 if (!(ids->put(&Tags_keys[0], num_Tags)))
394 throw DudleyException(msgPrefix+"put(Tags_keys)");
395
396 // Tags_names_*
397 // This is an array of strings, it should be stored as an array but
398 // instead I have hacked in one attribute per string because the NetCDF
399 // manual doesn't tell how to do an array of strings
400 int i = 0;
401 for (it = mesh->tagMap.begin(); it != mesh->tagMap.end(); it++, i++) {
402 stringstream ss;
403 ss << "Tags_name_" << i;
404 const string name(ss.str());
405 if (!dataFile.add_att(name.c_str(), it->first.c_str()))
406 throw DudleyException(msgPrefix+"add_att(Tags_names_XX)");
407 }
408 }
409
410 // Send token to next MPI process so he can take his turn
411 #ifdef ESYS_MPI
412 if (mpi_rank < mpi_size-1)
413 MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, getMPIComm());
414 #endif
415
416 // NetCDF file is closed by destructor of NcFile object
417
418 #else // USE_NETCDF
419 throw DudleyException("MeshAdapter::dump: not configured with netCDF. "
420 "Please contact your installation manager.");
421 #endif // USE_NETCDF
422 }
423
424 string MeshAdapter::getDescription() const
425 {
426 return "DudleyMesh";
427 }
428
429 string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
430 {
431 FunctionSpaceNamesMapType::iterator loc;
432 loc = m_functionSpaceTypeNames.find(functionSpaceType);
433 if (loc == m_functionSpaceTypeNames.end()) {
434 return "Invalid function space type code.";
435 } else {
436 return loc->second;
437 }
438 }
439
440 bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
441 {
442 FunctionSpaceNamesMapType::iterator loc;
443 loc = m_functionSpaceTypeNames.find(functionSpaceType);
444 return (loc != m_functionSpaceTypeNames.end());
445 }
446
447 void MeshAdapter::setFunctionSpaceTypeNames()
448 {
449 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
450 DegreesOfFreedom,"Dudley_DegreesOfFreedom [Solution(domain)]"));
451 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
452 Nodes,"Dudley_Nodes [ContinuousFunction(domain)]"));
453 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
454 Elements,"Dudley_Elements [Function(domain)]"));
455 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
456 ReducedElements,"Dudley_Reduced_Elements [ReducedFunction(domain)]"));
457 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
458 FaceElements,"Dudley_Face_Elements [FunctionOnBoundary(domain)]"));
459 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
460 ReducedFaceElements,"Dudley_Reduced_Face_Elements [ReducedFunctionOnBoundary(domain)]"));
461 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
462 Points,"Dudley_Points [DiracDeltaFunctions(domain)]"));
463 }
464
465 int MeshAdapter::getContinuousFunctionCode() const
466 {
467 return Nodes;
468 }
469
470 int MeshAdapter::getReducedContinuousFunctionCode() const
471 {
472 return Nodes;
473 }
474
475 int MeshAdapter::getFunctionCode() const
476 {
477 return Elements;
478 }
479
480 int MeshAdapter::getReducedFunctionCode() const
481 {
482 return ReducedElements;
483 }
484
485 int MeshAdapter::getFunctionOnBoundaryCode() const
486 {
487 return FaceElements;
488 }
489
490 int MeshAdapter::getReducedFunctionOnBoundaryCode() const
491 {
492 return ReducedFaceElements;
493 }
494
495 int MeshAdapter::getFunctionOnContactZeroCode() const
496 {
497 throw DudleyException("Dudley does not support contact elements.");
498 }
499
500 int MeshAdapter::getReducedFunctionOnContactZeroCode() const
501 {
502 throw DudleyException("Dudley does not support contact elements.");
503 }
504
505 int MeshAdapter::getFunctionOnContactOneCode() const
506 {
507 throw DudleyException("Dudley does not support contact elements.");
508 }
509
510 int MeshAdapter::getReducedFunctionOnContactOneCode() const
511 {
512 throw DudleyException("Dudley does not support contact elements.");
513 }
514
515 int MeshAdapter::getSolutionCode() const
516 {
517 return DegreesOfFreedom;
518 }
519
520 int MeshAdapter::getReducedSolutionCode() const
521 {
522 return DegreesOfFreedom;
523 }
524
525 int MeshAdapter::getDiracDeltaFunctionsCode() const
526 {
527 return Points;
528 }
529
530 int MeshAdapter::getDim() const
531 {
532 return m_dudleyMesh->getDim();
533 }
534
535 //
536 // Return the number of data points summed across all MPI processes
537 //
538 int MeshAdapter::getNumDataPointsGlobal() const
539 {
540 return m_dudleyMesh->Nodes->getGlobalNumNodes();
541 }
542
543 //
544 // return the number of data points per sample and the number of samples
545 // needed to represent data on a parts of the mesh.
546 //
547 pair<int,dim_t> MeshAdapter::getDataShape(int functionSpaceCode) const
548 {
549 int numDataPointsPerSample = 0;
550 dim_t numSamples = 0;
551 Mesh* mesh = getMesh();
552 switch (functionSpaceCode) {
553 case Nodes:
554 numDataPointsPerSample = 1;
555 numSamples = mesh->Nodes->getNumNodes();
556 break;
557 case Elements:
558 if (mesh->Elements) {
559 numSamples = mesh->Elements->numElements;
560 numDataPointsPerSample = mesh->Elements->numLocalDim+1;
561 }
562 break;
563 case ReducedElements:
564 if (mesh->Elements) {
565 numSamples = mesh->Elements->numElements;
566 numDataPointsPerSample =(mesh->Elements->numLocalDim==0)?0:1;
567 }
568 break;
569 case FaceElements:
570 if (mesh->FaceElements) {
571 numDataPointsPerSample = mesh->FaceElements->numLocalDim+1;
572 numSamples = mesh->FaceElements->numElements;
573 }
574 break;
575 case ReducedFaceElements:
576 if (mesh->FaceElements) {
577 numDataPointsPerSample = (mesh->FaceElements->numLocalDim==0)?0:1;
578 numSamples = mesh->FaceElements->numElements;
579 }
580 break;
581 case Points:
582 if (mesh->Points) {
583 numDataPointsPerSample = 1;
584 numSamples = mesh->Points->numElements;
585 }
586 break;
587 case DegreesOfFreedom:
588 if (mesh->Nodes) {
589 numDataPointsPerSample = 1;
590 numSamples = mesh->Nodes->getNumDegreesOfFreedom();
591 }
592 break;
593 default:
594 stringstream ss;
595 ss << "Invalid function space type: " << functionSpaceCode
596 << " for domain " << getDescription();
597 throw DudleyException(ss.str());
598 break;
599 }
600 return pair<int,dim_t>(numDataPointsPerSample,numSamples);
601 }
602
603 //
604 // adds linear PDE of second order into a given stiffness matrix and right hand side:
605 //
606 void MeshAdapter::addPDEToSystem(
607 escript::AbstractSystemMatrix& mat, escript::Data& rhs,
608 const escript::Data& A, const escript::Data& B, const escript::Data& C,
609 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
610 const escript::Data& d, const escript::Data& y,
611 const escript::Data& d_contact, const escript::Data& y_contact,
612 const escript::Data& d_dirac, const escript::Data& y_dirac) const
613 {
614 if (!d_contact.isEmpty() || !y_contact.isEmpty())
615 throw DudleyException("Dudley does not support contact elements");
616
617 #ifdef USE_TRILINOS
618 TrilinosMatrixAdapter* tm = dynamic_cast<TrilinosMatrixAdapter*>(&mat);
619 if (tm) {
620 tm->resumeFill();
621 }
622 #endif
623
624 Mesh* mesh = m_dudleyMesh.get();
625 Assemble_PDE(mesh->Nodes, mesh->Elements, mat.getPtr(), rhs,
626 A, B, C, D, X, Y);
627 Assemble_PDE(mesh->Nodes, mesh->FaceElements, mat.getPtr(), rhs,
628 escript::Data(), escript::Data(), escript::Data(), d,
629 escript::Data(), y);
630 Assemble_PDE(mesh->Nodes, mesh->Points, mat.getPtr(), rhs, escript::Data(),
631 escript::Data(), escript::Data(), d_dirac,
632 escript::Data(), y_dirac);
633
634 #ifdef USE_TRILINOS
635 if (tm) {
636 tm->fillComplete(true);
637 }
638 #endif
639 }
640
641 void MeshAdapter::addPDEToLumpedSystem(escript::Data& mat,
642 const escript::Data& D,
643 const escript::Data& d,
644 const escript::Data& d_dirac,
645 bool useHRZ) const
646 {
647 Mesh* mesh = m_dudleyMesh.get();
648 Assemble_LumpedSystem(mesh->Nodes, mesh->Elements, mat, D, useHRZ);
649 Assemble_LumpedSystem(mesh->Nodes, mesh->FaceElements, mat, d, useHRZ);
650 Assemble_LumpedSystem(mesh->Nodes, mesh->Points, mat, d_dirac, useHRZ);
651 }
652
653 //
654 // adds linear PDE of second order into the right hand side only
655 //
656 void MeshAdapter::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
657 const escript::Data& Y, const escript::Data& y,
658 const escript::Data& y_contact, const escript::Data& y_dirac) const
659 {
660 if (!y_contact.isEmpty())
661 throw DudleyException("Dudley does not support y_contact");
662
663 Mesh* mesh=m_dudleyMesh.get();
664
665 Assemble_PDE(mesh->Nodes, mesh->Elements, escript::ASM_ptr(), rhs,
666 escript::Data(), escript::Data(), escript::Data(),
667 escript::Data(), X, Y);
668
669 Assemble_PDE(mesh->Nodes, mesh->FaceElements, escript::ASM_ptr(), rhs,
670 escript::Data(), escript::Data(), escript::Data(),
671 escript::Data(), escript::Data(), y);
672
673 Assemble_PDE(mesh->Nodes, mesh->Points, escript::ASM_ptr(), rhs,
674 escript::Data(), escript::Data(), escript::Data(),
675 escript::Data(), escript::Data(), y_dirac);
676 }
677
678 //
679 // adds PDE of second order into a transport problem
680 //
681 void MeshAdapter::addPDEToTransportProblem(
682 escript::AbstractTransportProblem& tp, escript::Data& source,
683 const escript::Data& M, const escript::Data& A, const escript::Data& B,
684 const escript::Data& C, const escript::Data& D,const escript::Data& X,
685 const escript::Data& Y, const escript::Data& d, const escript::Data& y,
686 const escript::Data& d_contact,const escript::Data& y_contact,
687 const escript::Data& d_dirac, const escript::Data& y_dirac) const
688 {
689 if (!d_contact.isEmpty())
690 throw DudleyException("Dudley does not support d_contact");
691 if (!y_contact.isEmpty())
692 throw DudleyException("Dudley does not support y_contact");
693
694 paso::TransportProblem* ptp = dynamic_cast<paso::TransportProblem*>(&tp);
695 if (!ptp)
696 throw DudleyException("Dudley only accepts Paso transport problems");
697
698 source.expand();
699
700 Mesh* mesh = m_dudleyMesh.get();
701
702 Assemble_PDE(mesh->Nodes, mesh->Elements, ptp->borrowMassMatrix(), source,
703 escript::Data(), escript::Data(), escript::Data(), M,
704 escript::Data(), escript::Data());
705
706 Assemble_PDE(mesh->Nodes, mesh->Elements, ptp->borrowTransportMatrix(),
707 source, A, B, C, D, X, Y);
708
709 Assemble_PDE(mesh->Nodes, mesh->FaceElements, ptp->borrowTransportMatrix(),
710 source, escript::Data(), escript::Data(), escript::Data(), d,
711 escript::Data(), y);
712
713 Assemble_PDE(mesh->Nodes, mesh->Points, ptp->borrowTransportMatrix(),
714 source, escript::Data(), escript::Data(), escript::Data(),
715 d_dirac, escript::Data(), y_dirac);
716 }
717
718 //
719 // interpolates data between different function spaces:
720 //
721 void MeshAdapter::interpolateOnDomain(escript::Data& target,
722 const escript::Data& in) const
723 {
724 if (*in.getFunctionSpace().getDomain() != *this)
725 throw DudleyException("Illegal domain of interpolant.");
726 if (*target.getFunctionSpace().getDomain() != *this)
727 throw DudleyException("Illegal domain of interpolation target.");
728
729 Mesh* mesh = m_dudleyMesh.get();
730 switch (in.getFunctionSpace().getTypeCode()) {
731 case Nodes:
732 switch (target.getFunctionSpace().getTypeCode()) {
733 case Nodes:
734 case DegreesOfFreedom:
735 Assemble_CopyNodalData(mesh->Nodes, target, in);
736 break;
737 case Elements:
738 case ReducedElements:
739 Assemble_interpolate(mesh->Nodes, mesh->Elements, in, target);
740 break;
741 case FaceElements:
742 case ReducedFaceElements:
743 Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
744 break;
745 case Points:
746 Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
747 break;
748 default:
749 stringstream ss;
750 ss << "interpolateOnDomain: Dudley does not know anything "
751 "about function space type "
752 << target.getFunctionSpace().getTypeCode();
753 throw DudleyException(ss.str());
754 break;
755 }
756 break;
757 case Elements:
758 if (target.getFunctionSpace().getTypeCode() == Elements) {
759 Assemble_CopyElementData(mesh->Elements, target, in);
760 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
761 Assemble_AverageElementData(mesh->Elements, target, in);
762 } else {
763 throw DudleyException("No interpolation with data on elements possible.");
764 }
765 break;
766 case ReducedElements:
767 if (target.getFunctionSpace().getTypeCode() == ReducedElements) {
768 Assemble_CopyElementData(mesh->Elements, target, in);
769 } else {
770 throw DudleyException("No interpolation with data on elements "
771 "with reduced integration order possible.");
772 }
773 break;
774 case FaceElements:
775 if (target.getFunctionSpace().getTypeCode() == FaceElements) {
776 Assemble_CopyElementData(mesh->FaceElements, target, in);
777 } else if (target.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
778 Assemble_AverageElementData(mesh->FaceElements, target, in);
779 } else {
780 throw DudleyException("No interpolation with data on face elements possible.");
781 }
782 break;
783 case ReducedFaceElements:
784 if (target.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
785 Assemble_CopyElementData(mesh->FaceElements, target, in);
786 } else {
787 throw DudleyException("No interpolation with data on face "
788 "elements with reduced integration order possible.");
789 }
790 break;
791 case Points:
792 if (target.getFunctionSpace().getTypeCode() == Points) {
793 Assemble_CopyElementData(mesh->Points, target, in);
794 } else {
795 throw DudleyException("No interpolation with data on points possible.");
796 }
797 break;
798 case DegreesOfFreedom:
799 switch (target.getFunctionSpace().getTypeCode()) {
800 case DegreesOfFreedom:
801 Assemble_CopyNodalData(mesh->Nodes, target, in);
802 break;
803
804 case Nodes:
805 if (getMPISize() > 1) {
806 escript::Data temp = escript::Data(in);
807 temp.expand();
808 Assemble_CopyNodalData(mesh->Nodes, target, temp);
809 } else {
810 Assemble_CopyNodalData(mesh->Nodes, target, in);
811 }
812 break;
813 case Elements:
814 case ReducedElements:
815 if (getMPISize() > 1) {
816 escript::Data temp = escript::Data(in, continuousFunction(*this) );
817 Assemble_interpolate(mesh->Nodes, mesh->Elements, temp, target);
818 } else {
819 Assemble_interpolate(mesh->Nodes, mesh->Elements, in, target);
820 }
821 break;
822 case FaceElements:
823 case ReducedFaceElements:
824 if (getMPISize() > 1) {
825 escript::Data temp = escript::Data(in, continuousFunction(*this) );
826 Assemble_interpolate(mesh->Nodes, mesh->FaceElements, temp, target);
827 } else {
828 Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
829 }
830 break;
831 case Points:
832 if (getMPISize() > 1) {
833 //escript::Data temp=escript::Data(in, continuousFunction(*this) );
834 } else {
835 Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
836 }
837 break;
838 default:
839 stringstream ss;
840 ss << "interpolateOnDomain: Dudley does not know anything "
841 "about function space type "
842 << target.getFunctionSpace().getTypeCode();
843 throw DudleyException(ss.str());
844 break;
845 }
846 break;
847 default:
848 stringstream ss;
849 ss << "interpolateOnDomain: Dudley does not know anything about "
850 "function space type " << in.getFunctionSpace().getTypeCode();
851 throw DudleyException(ss.str());
852 break;
853 }
854 }
855
856 //
857 // copies the locations of sample points into x:
858 //
859 void MeshAdapter::setToX(escript::Data& arg) const
860 {
861 if (*arg.getFunctionSpace().getDomain() != *this)
862 throw DudleyException("setToX: Illegal domain of data point locations");
863
864 Mesh* mesh = m_dudleyMesh.get();
865 // in case of appropriate function space we can do the job directly:
866 if (arg.getFunctionSpace().getTypeCode() == Nodes) {
867 Assemble_NodeCoordinates(mesh->Nodes, arg);
868 } else {
869 escript::Data tmp_data = Vector(0., continuousFunction(*this), true);
870 Assemble_NodeCoordinates(mesh->Nodes, tmp_data);
871 // this is then interpolated onto arg:
872 interpolateOnDomain(arg, tmp_data);
873 }
874 }
875
876 //
877 // return the normal vectors at the location of data points as a Data object:
878 //
879 void MeshAdapter::setToNormal(escript::Data& normal) const
880 {
881 if (*normal.getFunctionSpace().getDomain() != *this)
882 throw ValueError("setToNormal: Illegal domain of normal locations");
883
884 Mesh* mesh=m_dudleyMesh.get();
885 if (normal.getFunctionSpace().getTypeCode() == FaceElements ||
886 normal.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
887 Assemble_getNormal(mesh->Nodes, mesh->FaceElements, normal);
888 } else {
889 stringstream ss;
890 ss << "setToNormal: Illegal function space type "
891 << normal.getFunctionSpace().getTypeCode();
892 throw ValueError(ss.str());
893 }
894 }
895
896 //
897 // interpolates data to other domain
898 //
899 void MeshAdapter::interpolateAcross(escript::Data& target,
900 const escript::Data& source) const
901 {
902 throw escript::NotImplementedError("Dudley does not allow interpolation "
903 "across domains.");
904 }
905
906 //
907 // calculates the integral of a function defined of arg:
908 //
909 void MeshAdapter::setToIntegrals(vector<double>& integrals,
910 const escript::Data& arg) const
911 {
912 if (*arg.getFunctionSpace().getDomain() != *this)
913 throw ValueError("setToIntegrals: Illegal domain of integration kernel");
914
915 Mesh* mesh = m_dudleyMesh.get();
916 switch (arg.getFunctionSpace().getTypeCode()) {
917 case Nodes: // fall through
918 case DegreesOfFreedom:
919 {
920 escript::Data temp(arg, escript::function(*this));
921 Assemble_integrate(mesh->Nodes, mesh->Elements, temp, integrals);
922 }
923 break;
924 case Elements: // fall through
925 case ReducedElements:
926 Assemble_integrate(mesh->Nodes,mesh->Elements, arg, integrals);
927 break;
928 case FaceElements: // fall through
929 case ReducedFaceElements:
930 Assemble_integrate(mesh->Nodes,mesh->FaceElements, arg, integrals);
931 break;
932 case Points:
933 throw ValueError("Integral of data on points is not supported.");
934 break;
935 default:
936 stringstream ss;
937 ss << "setToIntegrals: Dudley does not know anything about "
938 "function space type " << arg.getFunctionSpace().getTypeCode();
939 throw DudleyException(ss.str());
940 }
941 }
942
943 //
944 // calculates the gradient of arg:
945 //
946 void MeshAdapter::setToGradient(escript::Data& grad, const escript::Data& arg) const
947 {
948 if (*arg.getFunctionSpace().getDomain() != *this)
949 throw ValueError("setToGradient: Illegal domain of gradient argument");
950 if (*grad.getFunctionSpace().getDomain() != *this)
951 throw ValueError("setToGradient: Illegal domain of gradient");
952
953 Mesh* mesh = m_dudleyMesh.get();
954 escript::Data nodeData;
955 if (getMPISize() > 1) {
956 if (arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom) {
957 nodeData = escript::Data(arg, continuousFunction(*this));
958 } else {
959 nodeData = arg;
960 }
961 } else {
962 nodeData = arg;
963 }
964 switch (grad.getFunctionSpace().getTypeCode()) {
965 case Nodes:
966 throw DudleyException("Gradient at nodes is not supported.");
967 break;
968 case Elements:
969 Assemble_gradient(mesh->Nodes, mesh->Elements, grad, nodeData);
970 break;
971 case ReducedElements:
972 Assemble_gradient(mesh->Nodes, mesh->Elements, grad, nodeData);
973 break;
974 case FaceElements:
975 Assemble_gradient(mesh->Nodes,mesh->FaceElements, grad, nodeData);
976 break;
977 case ReducedFaceElements:
978 Assemble_gradient(mesh->Nodes, mesh->FaceElements, grad, nodeData);
979 break;
980 case Points:
981 throw DudleyException("Gradient at points is not supported.");
982 break;
983 case DegreesOfFreedom:
984 throw DudleyException("Gradient at degrees of freedom is not supported.");
985 break;
986 default:
987 stringstream ss;
988 ss << "Gradient: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
989 throw DudleyException(ss.str());
990 }
991 }
992
993 //
994 // returns the size of elements
995 //
996 void MeshAdapter::setToSize(escript::Data& size) const
997 {
998 Mesh* mesh=m_dudleyMesh.get();
999 switch (size.getFunctionSpace().getTypeCode()) {
1000 case Nodes:
1001 throw DudleyException("Size of nodes is not supported.");
1002 break;
1003 case Elements:
1004 Assemble_getSize(mesh->Nodes, mesh->Elements, size);
1005 break;
1006 case ReducedElements:
1007 Assemble_getSize(mesh->Nodes, mesh->Elements, size);
1008 break;
1009 case FaceElements:
1010 Assemble_getSize(mesh->Nodes, mesh->FaceElements, size);
1011 break;
1012 case ReducedFaceElements:
1013 Assemble_getSize(mesh->Nodes, mesh->FaceElements, size);
1014 break;
1015 case Points:
1016 throw DudleyException("Size of point elements is not supported.");
1017 break;
1018 case DegreesOfFreedom:
1019 throw DudleyException("Size of degrees of freedom is not supported.");
1020 break;
1021 default:
1022 stringstream ss;
1023 ss << "setToSize: Dudley does not know anything about function "
1024 "space type " << size.getFunctionSpace().getTypeCode();
1025 throw ValueError(ss.str());
1026 }
1027 }
1028
1029 //
1030 // sets the location of nodes
1031 //
1032 void MeshAdapter::setNewX(const escript::Data& newX)
1033 {
1034 if (*newX.getFunctionSpace().getDomain() != *this)
1035 throw DudleyException("Illegal domain of new point locations");
1036
1037 if (newX.getFunctionSpace() == continuousFunction(*this)) {
1038 m_dudleyMesh->setCoordinates(newX);
1039 } else {
1040 throw DudleyException("As of escript version 3.3 - setNewX only "
1041 "accepts ContinuousFunction arguments. Please interpolate.");
1042 }
1043 }
1044
1045 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1046 {
1047 #ifdef ESYS_MPI
1048 if (getMPISize() > 1) {
1049 Mesh* mesh = m_dudleyMesh.get();
1050 if (fs_code == DUDLEY_NODES) {
1051 const index_t myFirstNode = mesh->Nodes->getFirstNode();
1052 const index_t myLastNode = mesh->Nodes->getLastNode();
1053 const index_t k = mesh->Nodes->borrowGlobalNodesIndex()[id];
1054 return (myFirstNode <= k && k < myLastNode);
1055 } else {
1056 throw ValueError("ownSample: unsupported function space type");
1057 }
1058 }
1059 #endif
1060 return true;
1061 }
1062
1063 #ifdef USE_TRILINOS
1064 const_TrilinosGraph_ptr MeshAdapter::getTrilinosGraph() const
1065 {
1066 if (m_graph.is_null()) {
1067 m_graph = m_dudleyMesh->createTrilinosGraph();
1068 }
1069 return m_graph;
1070 }
1071 #endif
1072
1073 //
1074 // creates a stiffness matrix and initializes it with zeros
1075 //
1076 escript::ASM_ptr MeshAdapter::newSystemMatrix(int row_blocksize,
1077 const escript::FunctionSpace& row_functionspace,
1078 int column_blocksize,
1079 const escript::FunctionSpace& column_functionspace,
1080 int type) const
1081 {
1082 // is the domain right?
1083 if (*row_functionspace.getDomain() != *this)
1084 throw ValueError("domain of row function space does not match the domain of matrix generator.");
1085 if (*column_functionspace.getDomain() != *this)
1086 throw DudleyException("domain of column function space does not match the domain of matrix generator.");
1087
1088 // is the function space type right?
1089 if (row_functionspace.getTypeCode() != DegreesOfFreedom) {
1090 throw DudleyException("illegal function space type for system matrix rows.");
1091 }
1092 if (column_functionspace.getTypeCode() != DegreesOfFreedom) {
1093 throw DudleyException("illegal function space type for system matrix columns.");
1094 }
1095
1096 // generate matrix
1097 if (type & (int)SMT_TRILINOS) {
1098 #ifdef USE_TRILINOS
1099 const_TrilinosGraph_ptr graph(getTrilinosGraph());
1100 escript::ASM_ptr sm(new TrilinosMatrixAdapter(m_dudleyMesh->MPIInfo,
1101 row_blocksize, row_functionspace, graph));
1102 return sm;
1103 #else
1104 throw DudleyException("newSystemMatrix: dudley was not compiled "
1105 "with Trilinos support so the Trilinos solver stack cannot be "
1106 "used.");
1107 #endif
1108 } else if (type & (int)SMT_PASO) {
1109 paso::SystemMatrixPattern_ptr pattern(getMesh()->getPasoPattern());
1110 paso::SystemMatrix_ptr sm(new paso::SystemMatrix(type, pattern,
1111 row_blocksize, column_blocksize, false, row_functionspace,
1112 column_functionspace));
1113 return sm;
1114 } else {
1115 throw DudleyException("newSystemMatrix: unknown matrix type ID");
1116 }
1117 }
1118
1119 //
1120 // creates a TransportProblem
1121 //
1122 escript::ATP_ptr MeshAdapter::newTransportProblem(int blocksize,
1123 const escript::FunctionSpace& fs,
1124 int type) const
1125 {
1126 // is the domain right?
1127 if (*fs.getDomain() != *this)
1128 throw DudleyException("domain of function space does not match the domain of transport problem generator.");
1129 // is the function space type right
1130 if (fs.getTypeCode() != DegreesOfFreedom) {
1131 throw DudleyException("illegal function space type for system matrix rows.");
1132 }
1133
1134 // generate matrix
1135 paso::SystemMatrixPattern_ptr pattern(getMesh()->getPasoPattern());
1136 paso::TransportProblem_ptr transportProblem(new paso::TransportProblem(
1137 pattern, blocksize, fs));
1138 return transportProblem;
1139 }
1140
1141 //
1142 // returns true if data at the atom_type is considered as being cell centered:
1143 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1144 {
1145 switch (functionSpaceCode) {
1146 case Nodes:
1147 case DegreesOfFreedom:
1148 return false;
1149 break;
1150 case Elements:
1151 case FaceElements:
1152 case Points:
1153 case ReducedElements:
1154 case ReducedFaceElements:
1155 return true;
1156 break;
1157 default:
1158 stringstream ss;
1159 ss << "isCellOriented: Dudley does not know anything about "
1160 "function space type " << functionSpaceCode;
1161 throw ValueError(ss.str());
1162 }
1163 return false;
1164 }
1165
1166 bool
1167 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1168 {
1169 if (fs.empty())
1170 return false;
1171 // The idea is to use equivalence classes, i.e. types which can be
1172 // interpolated back and forth
1173 // class 1: DOF <-> Nodes
1174 // class 3: Points
1175 // class 4: Elements
1176 // class 5: ReducedElements
1177 // class 6: FaceElements
1178 // class 7: ReducedFaceElements
1179
1180 // There is also a set of lines. Interpolation is possible down a line but
1181 // not between lines.
1182 // class 1 and 2 belong to all lines so aren't considered.
1183 // line 0: class 3
1184 // line 1: class 4,5
1185 // line 2: class 6,7
1186
1187 // For classes with multiple members (class 1) we have vars to record
1188 // if there is at least one instance -> hasnodes is true if we have at
1189 // least one instance of Nodes.
1190 vector<int> hasclass(8);
1191 vector<int> hasline(3);
1192 bool hasnodes = false;
1193 for (int i = 0; i < fs.size(); ++i) {
1194 switch (fs[i]) {
1195 case Nodes:
1196 hasnodes = true; // fall through
1197 case DegreesOfFreedom:
1198 hasclass[1] = 1;
1199 break;
1200 case Points:
1201 hasline[0] = 1;
1202 hasclass[3] = 1;
1203 break;
1204 case Elements:
1205 hasclass[4] = 1;
1206 hasline[1] = 1;
1207 break;
1208 case ReducedElements:
1209 hasclass[5] = 1;
1210 hasline[1] = 1;
1211 break;
1212 case FaceElements:
1213 hasclass[6] = 1;
1214 hasline[2] = 1;
1215 break;
1216 case ReducedFaceElements:
1217 hasclass[7] = 1;
1218 hasline[2] = 1;
1219 break;
1220 default:
1221 return false;
1222 }
1223 }
1224 int totlines = hasline[0]+hasline[1]+hasline[2];
1225 // fail if we have more than one leaf group
1226 if (totlines > 1)
1227 // there are at least two branches we can't interpolate between
1228 return false;
1229
1230 if (totlines == 1) {
1231 if (hasline[0] == 1) // we have points
1232 resultcode = Points;
1233 else if (hasline[1] == 1) {
1234 if (hasclass[5]==1)
1235 resultcode=ReducedElements;
1236 else
1237 resultcode=Elements;
1238 } else if (hasline[2]==1) {
1239 if (hasclass[7]==1)
1240 resultcode=ReducedFaceElements;
1241 else
1242 resultcode=FaceElements;
1243 }
1244 } else { // totlines==0
1245 // something from class 1
1246 resultcode = (hasnodes ? Nodes : DegreesOfFreedom);
1247 }
1248 return true;
1249 }
1250
1251 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,
1252 int functionSpaceType_target) const
1253 {
1254 switch(functionSpaceType_source) {
1255 case Nodes:
1256 switch (functionSpaceType_target) {
1257 case Nodes:
1258 case DegreesOfFreedom:
1259 case Elements:
1260 case ReducedElements:
1261 case FaceElements:
1262 case ReducedFaceElements:
1263 case Points:
1264 return true;
1265 default:
1266 stringstream ss;
1267 ss << "Interpolation On Domain: Dudley does not know "
1268 "anything about function space type "
1269 << functionSpaceType_target;
1270 throw ValueError(ss.str());
1271 }
1272 break;
1273 case Elements:
1274 return (functionSpaceType_target == Elements ||
1275 functionSpaceType_target == ReducedElements);
1276 case ReducedElements:
1277 return (functionSpaceType_target == ReducedElements);
1278 case FaceElements:
1279 return (functionSpaceType_target == FaceElements ||
1280 functionSpaceType_target == ReducedFaceElements);
1281 case ReducedFaceElements:
1282 return (functionSpaceType_target == ReducedFaceElements);
1283 case Points:
1284 return (functionSpaceType_target == Points);
1285 case DegreesOfFreedom:
1286 switch (functionSpaceType_target) {
1287 case DegreesOfFreedom:
1288 case Nodes:
1289 case Elements:
1290 case ReducedElements:
1291 case Points:
1292 case FaceElements:
1293 case ReducedFaceElements:
1294 return true;
1295 default:
1296 stringstream ss;
1297 ss << "Interpolation On Domain: Dudley does not know "
1298 "anything about function space type "
1299 << functionSpaceType_target;
1300 throw DudleyException(ss.str());
1301 }
1302 break;
1303 default:
1304 stringstream ss;
1305 ss << "Interpolation On Domain: Dudley does not know anything "
1306 "about function space type " << functionSpaceType_source;
1307 throw DudleyException(ss.str());
1308 }
1309 return false;
1310 }
1311
1312 signed char MeshAdapter::preferredInterpolationOnDomain(
1313 int functionSpaceType_source,int functionSpaceType_target) const
1314 {
1315 if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1316 return 1;
1317 else if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1318 return -1;
1319 return 0;
1320 }
1321
1322 bool MeshAdapter::probeInterpolationAcross(int functionSpaceType_source,
1323 const AbstractDomain& targetDomain, int functionSpaceType_target) const
1324 {
1325 return false;
1326 }
1327
1328 bool MeshAdapter::operator==(const AbstractDomain& other) const
1329 {
1330 const MeshAdapter* temp = dynamic_cast<const MeshAdapter*>(&other);
1331 if (temp) {
1332 return (m_dudleyMesh == temp->m_dudleyMesh);
1333 }
1334 return false;
1335 }
1336
1337 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1338 {
1339 return !(operator==(other));
1340 }
1341
1342 int MeshAdapter::getSystemMatrixTypeId(const bp::object& options) const
1343 {
1344 const escript::SolverBuddy& sb = bp::extract<escript::SolverBuddy>(options);
1345
1346 int package = sb.getPackage();
1347 if (package == escript::SO_PACKAGE_TRILINOS) {
1348 #ifdef USE_TRILINOS
1349 return (int)SMT_TRILINOS;
1350 #else
1351 throw DudleyException("Trilinos requested but not built with Trilinos.");
1352 #endif
1353 }
1354 return (int)SMT_PASO | paso::SystemMatrix::getSystemMatrixTypeId(
1355 sb.getSolverMethod(), sb.getPreconditioner(), sb.getPackage(),
1356 sb.isSymmetric(), m_dudleyMesh->MPIInfo);
1357 }
1358
1359 int MeshAdapter::getTransportTypeId(int solver, int preconditioner,
1360 int package, bool symmetry) const
1361 {
1362 return paso::TransportProblem::getTypeId(solver, preconditioner, package,
1363 symmetry, getMPI());
1364 }
1365
1366 escript::Data MeshAdapter::getX() const
1367 {
1368 return continuousFunction(*this).getX();
1369 }
1370
1371 escript::Data MeshAdapter::getNormal() const
1372 {
1373 return functionOnBoundary(*this).getNormal();
1374 }
1375
1376 escript::Data MeshAdapter::getSize() const
1377 {
1378 return escript::function(*this).getSize();
1379 }
1380
1381 const index_t* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1382 {
1383 index_t* out = NULL;
1384 switch (functionSpaceType) {
1385 case Nodes:
1386 out = getMesh()->Nodes->Id;
1387 break;
1388 case Elements:
1389 out = getMesh()->Elements->Id;
1390 break;
1391 case ReducedElements:
1392 out = getMesh()->Elements->Id;
1393 break;
1394 case FaceElements:
1395 out = getMesh()->FaceElements->Id;
1396 break;
1397 case ReducedFaceElements:
1398 out = getMesh()->FaceElements->Id;
1399 break;
1400 case Points:
1401 out = getMesh()->Points->Id;
1402 break;
1403 case DegreesOfFreedom:
1404 out = getMesh()->Nodes->degreesOfFreedomId;
1405 break;
1406 default:
1407 stringstream ss;
1408 ss << "Invalid function space type: " << functionSpaceType
1409 << " for domain: " << getDescription();
1410 throw ValueError(ss.str());
1411 }
1412 return out;
1413 }
1414
1415 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, index_t sampleNo) const
1416 {
1417 int out = 0;
1418 switch (functionSpaceType) {
1419 case Nodes:
1420 out = getMesh()->Nodes->Tag[sampleNo];
1421 break;
1422 case Elements:
1423 out = getMesh()->Elements->Tag[sampleNo];
1424 break;
1425 case ReducedElements:
1426 out = getMesh()->Elements->Tag[sampleNo];
1427 break;
1428 case FaceElements:
1429 out = getMesh()->FaceElements->Tag[sampleNo];
1430 break;
1431 case ReducedFaceElements:
1432 out = getMesh()->FaceElements->Tag[sampleNo];
1433 break;
1434 case Points:
1435 out = getMesh()->Points->Tag[sampleNo];
1436 break;
1437 case DegreesOfFreedom:
1438 throw DudleyException("DegreesOfFreedom does not support tags.");
1439 break;
1440 default:
1441 stringstream ss;
1442 ss << "Invalid function space type: " << functionSpaceType
1443 << " for domain: " << getDescription();
1444 throw DudleyException(ss.str());
1445 }
1446 return out;
1447 }
1448
1449
1450 void MeshAdapter::setTags(int functionSpaceType, int newTag, const escript::Data& mask) const
1451 {
1452 switch (functionSpaceType) {
1453 case Nodes:
1454 getMesh()->Nodes->setTags(newTag, mask);
1455 break;
1456 case DegreesOfFreedom:
1457 throw DudleyException("DegreesOfFreedom does not support tags");
1458 break;
1459 case Elements: // fall through
1460 case ReducedElements:
1461 getMesh()->Elements->setTags(newTag, mask);
1462 break;
1463 case FaceElements:
1464 case ReducedFaceElements:
1465 getMesh()->FaceElements->setTags(newTag, mask);
1466 break;
1467 case Points:
1468 getMesh()->Points->setTags(newTag, mask);
1469 break;
1470 default:
1471 stringstream ss;
1472 ss << "Dudley does not know anything about function space type "
1473 << functionSpaceType;
1474 throw ValueError(ss.str());
1475 }
1476 }
1477
1478 void MeshAdapter::setTagMap(const string& name, int tag)
1479 {
1480 getMesh()->addTagMap(name, tag);
1481 }
1482
1483 int MeshAdapter::getTag(const string& name) const
1484 {
1485 return getMesh()->getTag(name);
1486 }
1487
1488 bool MeshAdapter::isValidTagName(const string& name) const
1489 {
1490 return getMesh()->isValidTagName(name);
1491 }
1492
1493 string MeshAdapter::showTagNames() const
1494 {
1495 stringstream ss;
1496 TagMap::const_iterator it = getMesh()->tagMap.begin();
1497 while (it != getMesh()->tagMap.end()) {
1498 ss << it->first;
1499 ++it;
1500 if (it != getMesh()->tagMap.end())
1501 ss << ", ";
1502 }
1503 return ss.str();
1504 }
1505
1506 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
1507 {
1508 switch (functionSpaceCode) {
1509 case Nodes:
1510 return getMesh()->Nodes->tagsInUse.size();
1511 case DegreesOfFreedom:
1512 throw ValueError("DegreesOfFreedom does not support tags");
1513 case Elements: // fall through
1514 case ReducedElements:
1515 return getMesh()->Elements->tagsInUse.size();
1516 case FaceElements: // fall through
1517 case ReducedFaceElements:
1518 return getMesh()->FaceElements->tagsInUse.size();
1519 case Points:
1520 return getMesh()->Points->tagsInUse.size();
1521 default:
1522 stringstream ss;
1523 ss << "Dudley does not know anything about function space type "
1524 << functionSpaceCode;
1525 throw ValueError(ss.str());
1526 }
1527 return 0;
1528 }
1529
1530 const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
1531 {
1532 switch (functionSpaceCode) {
1533 case Nodes:
1534 if (getMesh()->Nodes->tagsInUse.empty())
1535 return NULL;
1536 else
1537 return &getMesh()->Nodes->tagsInUse[0];
1538 case DegreesOfFreedom:
1539 throw DudleyException("DegreesOfFreedom does not support tags");
1540 case Elements: // fall through
1541 case ReducedElements:
1542 if (getMesh()->Elements->tagsInUse.empty())
1543 return NULL;
1544 else
1545 return &getMesh()->Elements->tagsInUse[0];
1546 case FaceElements: // fall through
1547 case ReducedFaceElements:
1548 if (getMesh()->FaceElements->tagsInUse.empty())
1549 return NULL;
1550 else
1551 return &getMesh()->FaceElements->tagsInUse[0];
1552 case Points:
1553 if (getMesh()->Points->tagsInUse.empty())
1554 return NULL;
1555 else
1556 return &getMesh()->Points->tagsInUse[0];
1557 default:
1558 stringstream ss;
1559 ss << "Dudley does not know anything about function space type "
1560 << functionSpaceCode;
1561 throw DudleyException(ss.str());
1562 }
1563 return NULL;
1564 }
1565
1566
1567 bool MeshAdapter::canTag(int functionSpaceCode) const
1568 {
1569 switch(functionSpaceCode) {
1570 case Nodes:
1571 case Elements:
1572 case ReducedElements:
1573 case FaceElements:
1574 case ReducedFaceElements:
1575 case Points:
1576 return true;
1577 default:
1578 return false;
1579 }
1580 }
1581
1582 MeshAdapter::StatusType MeshAdapter::getStatus() const
1583 {
1584 return getMesh()->getStatus();
1585 }
1586
1587 int MeshAdapter::getApproximationOrder(int functionSpaceCode) const
1588 {
1589 switch (functionSpaceCode) {
1590 case Nodes:
1591 case DegreesOfFreedom:
1592 return getMesh()->approximationOrder;
1593 case Elements:
1594 case FaceElements:
1595 case Points:
1596 return getMesh()->integrationOrder;
1597 case ReducedElements:
1598 case ReducedFaceElements:
1599 return getMesh()->reducedIntegrationOrder;
1600 default:
1601 stringstream ss;
1602 ss << "Dudley does not know anything about function space type "
1603 << functionSpaceCode;
1604 throw ValueError(ss.str());
1605 }
1606 return 0;
1607 }
1608
1609 bool MeshAdapter::supportsContactElements() const
1610 {
1611 return false;
1612 }
1613
1614 escript::Data MeshAdapter::randomFill(
1615 const escript::DataTypes::ShapeType& shape,
1616 const escript::FunctionSpace& what, long seed,
1617 const bp::tuple& filter) const
1618 {
1619 escript::Data towipe(0, shape, what, true);
1620 // since we just made this object, no sharing is possible and we don't
1621 // need to check for exclusive write
1622 escript::DataTypes::RealVectorType& dv(towipe.getExpandedVectorReference());
1623 escript::randomFillArray(seed, &dv[0], dv.size());
1624 return towipe;
1625 }
1626
1627 } // end of namespace
1628

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26