/[escript]/trunk/finley/src/FinleyDomain.cpp
ViewVC logotype

Contents of /trunk/finley/src/FinleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 1 week ago) by uqaeller
File size: 126322 byte(s)
Updated the copyright header.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
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-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17
18 #include "FinleyDomain.h"
19 #include "Assemble.h"
20 #include "FinleyException.h"
21 #include "IndexList.h"
22 #include "ReferenceElements.h"
23
24 #include <escript/Data.h>
25 #include <escript/DataFactory.h>
26 #include <escript/Random.h>
27 #include <escript/SolverOptions.h>
28
29 #ifdef ESYS_HAVE_PASO
30 #include <paso/SystemMatrix.h>
31 #include <paso/Transport.h>
32 #endif
33
34 #ifdef ESYS_HAVE_TRILINOS
35 #include <trilinoswrap/TrilinosMatrixAdapter.h>
36
37 using esys_trilinos::TrilinosMatrixAdapter;
38 using esys_trilinos::const_TrilinosGraph_ptr;
39 #endif
40
41 #include <boost/scoped_array.hpp>
42
43 #ifdef ESYS_HAVE_NETCDF
44 #ifdef NETCDF4
45 #include <ncVar.h>
46 #include <ncDim.h>
47 #include <escript/NCHelper.h>
48 #else
49 #include <netcdfcpp.h>
50 #endif
51 #endif
52
53
54
55 using namespace std;
56 namespace bp = boost::python;
57 using escript::NotImplementedError;
58 using escript::ValueError;
59
60 #ifdef NETCDF4
61 using namespace netCDF;
62 #endif
63
64
65 namespace finley {
66
67 using escript::DataTypes::real_t;
68 using escript::DataTypes::cplx_t;
69
70 // define the static constants
71 FinleyDomain::FunctionSpaceNamesMapType FinleyDomain::m_functionSpaceTypeNames;
72
73 FinleyDomain::FinleyDomain(const string& name, int numDim, escript::JMPI jmpi) :
74 m_mpiInfo(jmpi),
75 m_name(name),
76 approximationOrder(-1),
77 reducedApproximationOrder(-1),
78 integrationOrder(-1),
79 reducedIntegrationOrder(-1),
80 m_elements(NULL),
81 m_faceElements(NULL),
82 m_contactElements(NULL),
83 m_points(NULL)
84 {
85 // allocate node table
86 m_nodes = new NodeFile(numDim, m_mpiInfo);
87 setFunctionSpaceTypeNames();
88 }
89
90 FinleyDomain::FinleyDomain(const FinleyDomain& in) :
91 m_mpiInfo(in.m_mpiInfo),
92 m_name(in.m_name),
93 approximationOrder(in.approximationOrder),
94 reducedApproximationOrder(in.reducedApproximationOrder),
95 integrationOrder(in.integrationOrder),
96 reducedIntegrationOrder(in.reducedIntegrationOrder),
97 m_nodes(in.m_nodes),
98 m_elements(in.m_elements),
99 m_faceElements(in.m_faceElements),
100 m_contactElements(in.m_contactElements),
101 m_points(in.m_points)
102 {
103 setFunctionSpaceTypeNames();
104 }
105
106 FinleyDomain::~FinleyDomain()
107 {
108 delete m_nodes;
109 delete m_elements;
110 delete m_faceElements;
111 delete m_contactElements;
112 delete m_points;
113 }
114
115 void FinleyDomain::MPIBarrier() const
116 {
117 #ifdef ESYS_MPI
118 MPI_Barrier(getMPIComm());
119 #endif
120 }
121
122 void FinleyDomain::setElements(ElementFile* elements)
123 {
124 delete m_elements;
125 m_elements = elements;
126 }
127
128 void FinleyDomain::setFaceElements(ElementFile* elements)
129 {
130 delete m_faceElements;
131 m_faceElements = elements;
132 }
133
134 void FinleyDomain::setContactElements(ElementFile* elements)
135 {
136 delete m_contactElements;
137 m_contactElements = elements;
138 }
139
140 void FinleyDomain::setPoints(ElementFile* elements)
141 {
142 delete m_points;
143 m_points = elements;
144 }
145
146 void FinleyDomain::setOrders()
147 {
148 const int ORDER_MAX = 9999999;
149 int locals[4] = { ORDER_MAX, ORDER_MAX, ORDER_MAX, ORDER_MAX };
150
151 if (m_elements != NULL && m_elements->numElements > 0) {
152 locals[0] = std::min(locals[0], m_elements->referenceElementSet->referenceElement->BasisFunctions->Type->numOrder);
153 locals[1] = std::min(locals[1], m_elements->referenceElementSet->referenceElement->LinearBasisFunctions->Type->numOrder);
154 locals[2] = std::min(locals[2], m_elements->referenceElementSet->referenceElement->integrationOrder);
155 locals[3] = std::min(locals[3], m_elements->referenceElementSet->referenceElementReducedQuadrature->integrationOrder);
156 }
157 if (m_faceElements != NULL && m_faceElements->numElements > 0) {
158 locals[0] = std::min(locals[0], m_faceElements->referenceElementSet->referenceElement->BasisFunctions->Type->numOrder);
159 locals[1] = std::min(locals[1], m_faceElements->referenceElementSet->referenceElement->LinearBasisFunctions->Type->numOrder);
160 locals[2] = std::min(locals[2], m_faceElements->referenceElementSet->referenceElement->integrationOrder);
161 locals[3] = std::min(locals[3], m_faceElements->referenceElementSet->referenceElementReducedQuadrature->integrationOrder);
162 }
163 if (m_contactElements != NULL && m_contactElements->numElements > 0) {
164 locals[0] = std::min(locals[0], m_contactElements->referenceElementSet->referenceElement->BasisFunctions->Type->numOrder);
165 locals[1] = std::min(locals[1], m_contactElements->referenceElementSet->referenceElement->LinearBasisFunctions->Type->numOrder);
166 locals[2] = std::min(locals[2], m_contactElements->referenceElementSet->referenceElement->integrationOrder);
167 locals[3] = std::min(locals[3], m_contactElements->referenceElementSet->referenceElementReducedQuadrature->integrationOrder);
168 }
169
170 #ifdef ESYS_MPI
171 int globals[4];
172 MPI_Allreduce(locals, globals, 4, MPI_INT, MPI_MIN, m_mpiInfo->comm);
173 approximationOrder = (globals[0] < ORDER_MAX ? globals[0] : -1);
174 reducedApproximationOrder = (globals[1] < ORDER_MAX ? globals[1] : -1);
175 integrationOrder = (globals[2] < ORDER_MAX ? globals[2] : -1);
176 reducedIntegrationOrder = (globals[3] < ORDER_MAX ? globals[3] : -1);
177 #else
178 approximationOrder = (locals[0] < ORDER_MAX ? locals[0] : -1);
179 reducedApproximationOrder = (locals[1] < ORDER_MAX ? locals[1] : -1);
180 integrationOrder = (locals[2] < ORDER_MAX ? locals[2] : -1);
181 reducedIntegrationOrder = (locals[3] < ORDER_MAX ? locals[3] : -1);
182 #endif
183 }
184
185 void FinleyDomain::createMappings(const IndexVector& dofDist,
186 const IndexVector& nodeDist)
187 {
188 std::vector<short> maskReducedNodes(m_nodes->getNumNodes(), -1);
189 markNodes(maskReducedNodes, 0, true);
190 IndexVector indexReducedNodes = util::packMask(maskReducedNodes);
191 m_nodes->createNodeMappings(indexReducedNodes, dofDist, nodeDist);
192 }
193
194 void FinleyDomain::markNodes(vector<short>& mask, index_t offset,
195 bool useLinear) const
196 {
197 m_elements->markNodes(mask, offset, useLinear);
198 m_faceElements->markNodes(mask, offset, useLinear);
199 m_contactElements->markNodes(mask, offset, useLinear);
200 m_points->markNodes(mask, offset, useLinear);
201 }
202
203 void FinleyDomain::relabelElementNodes(const IndexVector& newNode, index_t offset)
204 {
205 m_elements->relabelNodes(newNode, offset);
206 m_faceElements->relabelNodes(newNode, offset);
207 m_contactElements->relabelNodes(newNode, offset);
208 m_points->relabelNodes(newNode, offset);
209 }
210
211 #ifdef NETCDF4
212
213 void FinleyDomain::dump(const string& fileName) const
214 {
215 #ifdef ESYS_HAVE_NETCDF
216 NcDim ncdims[12];
217 NcVar ids;
218 index_t* index_ptr;
219 #ifdef ESYS_INDEXTYPE_LONG
220 NcType ncIdxType = ncInt64; // "ncLong" is deprecated
221 #else
222 NcType ncIdxType = ncInt;
223 #endif
224 int num_Tags = 0;
225 int mpi_size = getMPISize();
226 int mpi_rank = getMPIRank();
227 int numDim = m_nodes->numDim;
228 dim_t numNodes = m_nodes->getNumNodes();
229 dim_t num_Elements = m_elements->numElements;
230 dim_t num_FaceElements = m_faceElements->numElements;
231 dim_t num_ContactElements = m_contactElements->numElements;
232 dim_t num_Points = m_points->numElements;
233 int num_Elements_numNodes = m_elements->numNodes;
234 int num_FaceElements_numNodes = m_faceElements->numNodes;
235 int num_ContactElements_numNodes = m_contactElements->numNodes;
236 #ifdef ESYS_MPI
237 MPI_Status status;
238 #endif
239
240 // Incoming token indicates it's my turn to write
241 #ifdef ESYS_MPI
242 if (mpi_rank > 0)
243 MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, getMPIComm(), &status);
244 #endif
245
246 const string newFileName(m_mpiInfo->appendRankToFileName(fileName));
247
248 // Figure out how much storage is required for tags
249 num_Tags = m_tagMap.size();
250
251 NcFile dataFile;
252 try
253 {
254 dataFile.open(newFileName.c_str(), NcFile::FileMode::replace, NcFile::FileFormat::nc4);
255 }
256 catch (exceptions::NcException* e)
257 {
258 throw FinleyException("Error - FinleyDomain:: opening of netCDF file for output failed.");
259 }
260
261
262 string msgPrefix("Error in FinleyDomain::dump: NetCDF operation failed - ");
263 // Define dimensions (num_Elements and dim_Elements are identical,
264 // dim_Elements only appears if > 0)
265 if ((ncdims[0] = dataFile.addDim("numNodes", numNodes)).isNull() )
266 throw FinleyException(msgPrefix+"add_dim(numNodes)");
267 if ((ncdims[1] = dataFile.addDim("numDim", numDim)).isNull() )
268 throw FinleyException(msgPrefix+"add_dim(numDim)");
269 if ((ncdims[2] = dataFile.addDim("mpi_size_plus_1", mpi_size+1)).isNull() )
270 throw FinleyException(msgPrefix+"add_dim(mpi_size)");
271 if (num_Elements > 0)
272 if ((ncdims[3] = dataFile.addDim("dim_Elements", num_Elements)).isNull() )
273 throw FinleyException(msgPrefix+"add_dim(dim_Elements)");
274 if (num_FaceElements > 0)
275 if ((ncdims[4] = dataFile.addDim("dim_FaceElements", num_FaceElements)).isNull() )
276 throw FinleyException(msgPrefix+"add_dim(dim_FaceElements)");
277 if (num_ContactElements > 0)
278 if ((ncdims[5] = dataFile.addDim("dim_ContactElements", num_ContactElements)).isNull() )
279 throw FinleyException(msgPrefix+"add_dim(dim_ContactElements)");
280 if (num_Points > 0)
281 if ((ncdims[6] = dataFile.addDim("dim_Points", num_Points)).isNull() )
282 throw FinleyException(msgPrefix+"add_dim(dim_Points)");
283 if (num_Elements > 0)
284 if ((ncdims[7] = dataFile.addDim("dim_Elements_Nodes", num_Elements_numNodes)).isNull() )
285 throw FinleyException(msgPrefix+"add_dim(dim_Elements_Nodes)");
286 if (num_FaceElements > 0)
287 if ((ncdims[8] = dataFile.addDim("dim_FaceElements_numNodes", num_FaceElements_numNodes)).isNull() )
288 throw FinleyException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
289 if (num_ContactElements > 0)
290 if ((ncdims[9] = dataFile.addDim("dim_ContactElements_numNodes", num_ContactElements_numNodes)).isNull() )
291 throw FinleyException(msgPrefix+"add_dim(dim_ContactElements_numNodes)");
292 if (num_Tags > 0)
293 if ((ncdims[10] = dataFile.addDim("dim_Tags", num_Tags)).isNull() )
294 throw FinleyException(msgPrefix+"add_dim(dim_Tags)");
295
296 // Attributes: MPI size, MPI rank, Name, order, reduced_order
297 NcInt ni;
298 if (dataFile.putAtt("index_size", ni, (int)sizeof(index_t)).isNull())
299 throw FinleyException(msgPrefix+"putAtt(index_size)");
300 if (dataFile.putAtt("mpi_size", ni, mpi_size).isNull())
301 throw FinleyException(msgPrefix+"putAtt(mpi_size)");
302 if (dataFile.putAtt("mpi_rank", ni, mpi_rank).isNull())
303 throw FinleyException(msgPrefix+"putAtt(mpi_rank)");
304 if (dataFile.putAtt("Name", m_name).isNull())
305 throw FinleyException(msgPrefix+"putAtt(Name)");
306 if (dataFile.putAtt("numDim", ni, numDim).isNull())
307 throw FinleyException(msgPrefix+"putAtt(order)");
308 if (dataFile.putAtt("order", ni, integrationOrder).isNull())
309 throw FinleyException(msgPrefix+"putAtt(order)");
310 if (dataFile.putAtt("reduced_order", ni, reducedIntegrationOrder).isNull())
311 throw FinleyException(msgPrefix+"putAtt(reduced_order)");
312 if (dataFile.putAtt("numNodes", ni, numNodes).isNull())
313 throw FinleyException(msgPrefix+"putAtt(numNodes)");
314 if (dataFile.putAtt("num_Elements", ni, num_Elements).isNull())
315 throw FinleyException(msgPrefix+"putAtt(num_Elements)");
316 if (dataFile.putAtt("num_FaceElements", ni, num_FaceElements).isNull())
317 throw FinleyException(msgPrefix+"putAtt(num_FaceElements)");
318 if (dataFile.putAtt("num_ContactElements", ni, num_ContactElements).isNull())
319 throw FinleyException(msgPrefix+"putAtt(num_ContactElements)");
320 if (dataFile.putAtt("num_Points", ni, num_Points).isNull())
321 throw FinleyException(msgPrefix+"putAtt(num_Points)");
322 if (dataFile.putAtt("num_Elements_numNodes", ni, num_Elements_numNodes).isNull())
323 throw FinleyException(msgPrefix+"putAtt(num_Elements_numNodes)");
324 if (dataFile.putAtt("num_FaceElements_numNodes", ni, num_FaceElements_numNodes).isNull())
325 throw FinleyException(msgPrefix+"putAtt(num_FaceElements_numNodes)");
326 if (dataFile.putAtt("num_ContactElements_numNodes", ni, num_ContactElements_numNodes).isNull())
327 throw FinleyException(msgPrefix+"putAtt(num_ContactElements_numNodes)");
328 if (dataFile.putAtt("Elements_TypeId", ni, m_elements->referenceElementSet->referenceElement->Type->TypeId).isNull() )
329 throw FinleyException(msgPrefix+"putAtt(Elements_TypeId)");
330 if (dataFile.putAtt("FaceElements_TypeId", ni, m_faceElements->referenceElementSet->referenceElement->Type->TypeId).isNull() )
331 throw FinleyException(msgPrefix+"putAtt(FaceElements_TypeId)");
332 if (dataFile.putAtt("ContactElements_TypeId", ni, m_contactElements->referenceElementSet->referenceElement->Type->TypeId).isNull() )
333 throw FinleyException(msgPrefix+"putAtt(ContactElements_TypeId)");
334 if (dataFile.putAtt("Points_TypeId", ni, m_points->referenceElementSet->referenceElement->Type->TypeId).isNull() )
335 throw FinleyException(msgPrefix+"putAtt(Points_TypeId)");
336 if (dataFile.putAtt("num_Tags", ni, num_Tags).isNull())
337 throw FinleyException(msgPrefix+"putAtt(num_Tags)");
338
339 // // // // // Nodes // // // // //
340
341 // Nodes nodeDistribution
342 if ((ids = dataFile.addVar("Nodes_NodeDistribution", ncIdxType, ncdims[2])).isNull() )
343 throw FinleyException(msgPrefix+"add_var(Nodes_NodeDistribution)");
344 index_ptr = &m_nodes->nodesDistribution->first_component[0];
345 ids.putVar(index_ptr); // don't think I need to specify bounds here but it should be mpi_size+1
346
347 // Nodes degreesOfFreedomDistribution
348 if (( ids = dataFile.addVar("Nodes_DofDistribution", ncIdxType, ncdims[2])).isNull() )
349 throw FinleyException(msgPrefix+"add_var(Nodes_DofDistribution)");
350 index_ptr = &m_nodes->degreesOfFreedomDistribution->first_component[0];
351 ids.putVar(index_ptr); // don't think I need to specify bounds here but it should be mpi_size+1
352
353 // Only write nodes if non-empty because NetCDF doesn't like empty arrays
354 // (it treats them as NC_UNLIMITED)
355 if (numNodes > 0) {
356 // Nodes Id
357 if (( ids = dataFile.addVar("Nodes_Id", ncIdxType, ncdims[0])).isNull() )
358 throw FinleyException(msgPrefix+"add_var(Nodes_Id)");
359 ids.putVar(&m_nodes->Id[0]); // should be numNodes values written
360
361 // Nodes Tag
362 if (( ids = dataFile.addVar("Nodes_Tag", ncInt, ncdims[0])).isNull() )
363 throw FinleyException(msgPrefix+"add_var(Nodes_Tag)");
364 ids.putVar(&m_nodes->Tag[0]); // should be numNodes values written
365
366 // Nodes gDOF
367 if (( ids = dataFile.addVar("Nodes_gDOF", ncIdxType, ncdims[0])).isNull() )
368 throw FinleyException(msgPrefix+"add_var(Nodes_gDOF)");
369 ids.putVar(&m_nodes->globalDegreesOfFreedom[0]); // should be numNodes values written
370
371 // Nodes global node index
372 if (( ids = dataFile.addVar("Nodes_gNI", ncIdxType, ncdims[0])).isNull() )
373 throw FinleyException(msgPrefix+"add_var(Nodes_gNI)");
374 ids.putVar(&m_nodes->globalNodesIndex[0]); // should be numNodes values written
375
376 // Nodes grDof
377 if (( ids = dataFile.addVar("Nodes_grDfI", ncIdxType, ncdims[0])).isNull() )
378 throw FinleyException(msgPrefix+"add_var(Nodes_grDfI)");
379 ids.putVar(&m_nodes->globalReducedDOFIndex[0]); // should be numNodes values written
380
381 // Nodes grNI
382 if (( ids = dataFile.addVar("Nodes_grNI", ncIdxType, ncdims[0])).isNull() )
383 throw FinleyException(msgPrefix+"add_var(Nodes_grNI)");
384 ids.putVar(&m_nodes->globalReducedNodesIndex[0]); // should be numNodes values written
385
386 // Nodes Coordinates
387 vector<NcDim> ncds;
388 ncds.push_back(ncdims[0]);
389 ncds.push_back(ncdims[1]);
390 if (( ids = dataFile.addVar("Nodes_Coordinates", ncDouble, ncds)).isNull() )
391 throw FinleyException(msgPrefix+"add_var(Nodes_Coordinates)");
392 ids.putVar(m_nodes->Coordinates); // should be (numNodes, numDim) values written
393 }
394
395 // // // // // Elements // // // // //
396 if (num_Elements > 0) {
397 // Elements_Id
398 if (( ids = dataFile.addVar("Elements_Id", ncIdxType, ncdims[3])).isNull() )
399 throw FinleyException(msgPrefix+"add_var(Elements_Id)");
400 ids.putVar(m_elements->Id); // write num_Elements values
401
402 // Elements_Tag
403 if (( ids = dataFile.addVar("Elements_Tag", ncInt, ncdims[3])).isNull() )
404 throw FinleyException(msgPrefix+"add_var(Elements_Tag)");
405 ids.putVar(m_elements->Tag); // write num_Elements values
406
407 // Elements_Owner
408 if (( ids = dataFile.addVar("Elements_Owner", ncInt, ncdims[3])).isNull() )
409 throw FinleyException(msgPrefix+"add_var(Elements_Owner)");
410 ids.putVar(m_elements->Owner); // write num_Elements values
411
412 // Elements_Color
413 if (( ids = dataFile.addVar("Elements_Color", ncInt, ncdims[3])).isNull() )
414 throw FinleyException(msgPrefix+"add_var(Elements_Color)");
415 ids.putVar(m_elements->Color); // write num_Elements values
416
417 // Elements_Nodes
418 vector<NcDim> dv;
419 dv.push_back(ncdims[3]);
420 dv.push_back(ncdims[7]);
421 if (( ids = dataFile.addVar("Elements_Nodes", ncIdxType, dv) ).isNull() )
422 throw FinleyException(msgPrefix+"add_var(Elements_Nodes)");
423 ids.putVar(&m_elements->Nodes[0]); //(, num_Elements, num_Elements_numNodes) values written
424 }
425
426 // // // // // Face_Elements // // // // //
427 if (num_FaceElements > 0) {
428 // FaceElements_Id
429 if ((ids = dataFile.addVar("FaceElements_Id", ncIdxType, ncdims[4])).isNull())
430 throw FinleyException(msgPrefix+"add_var(FaceElements_Id)");
431 ids.putVar(m_faceElements->Id); // num_FaceElements values written
432
433 // FaceElements_Tag
434 if ((ids = dataFile.addVar("FaceElements_Tag", ncInt, ncdims[4])).isNull())
435 throw FinleyException(msgPrefix+"add_var(FaceElements_Tag)");
436 ids.putVar(m_faceElements->Tag); // num_FaceElements
437
438 // FaceElements_Owner
439 if ((ids = dataFile.addVar("FaceElements_Owner", ncInt, ncdims[4])).isNull())
440 throw FinleyException(msgPrefix+"add_var(FaceElements_Owner)");
441 ids.putVar(m_faceElements->Owner); // num_FaceElements
442
443 // FaceElements_Color
444 if ((ids = dataFile.addVar("FaceElements_Color", ncIdxType, ncdims[4])).isNull())
445 throw FinleyException(msgPrefix+"add_var(FaceElements_Color)");
446 ids.putVar(m_faceElements->Color); // num_FaceElements
447
448 // FaceElements_Nodes
449 vector<NcDim> dv;
450 dv.push_back(ncdims[4]);
451 dv.push_back(ncdims[8]);
452 if ((ids = dataFile.addVar("FaceElements_Nodes", ncIdxType, dv)).isNull())
453 throw FinleyException(msgPrefix+"add_var(FaceElements_Nodes)");
454 ids.putVar(m_faceElements->Nodes); // num_FaceElements_numNodes
455
456 }
457
458 // // // // // Contact_Elements // // // // //
459 if (num_ContactElements > 0) {
460 // ContactElements_Id
461 if ((ids = dataFile.addVar("ContactElements_Id", ncIdxType, ncdims[5])).isNull())
462 throw FinleyException(msgPrefix+"add_var(ContactElements_Id)");
463 ids.putVar(m_contactElements->Id); // num_ContactElements values written
464
465 // ContactElements_Tag
466 if ((ids = dataFile.addVar("ContactElements_Tag", ncInt, ncdims[5])).isNull())
467 throw FinleyException(msgPrefix+"add_var(ContactElements_Tag)");
468 ids.putVar(m_contactElements->Tag); // num_ContactElements values written
469
470 // ContactElements_Owner
471 if ((ids = dataFile.addVar("ContactElements_Owner", ncInt, ncdims[5])).isNull())
472 throw FinleyException(msgPrefix+"add_var(ContactElements_Owner)");
473 ids.putVar(m_contactElements->Owner); // num_ContactElements values written
474
475 // ContactElements_Color
476 if ((ids = dataFile.addVar("ContactElements_Color", ncInt, ncdims[5])).isNull())
477 throw FinleyException(msgPrefix+"add_var(ContactElements_Color)");
478 ids.putVar(m_contactElements->Color); // num_ContactElements values written
479
480 // ContactElements_Nodes
481 vector<NcDim> dv;
482 dv.push_back(ncdims[5]);
483 dv.push_back(ncdims[9]);
484 if ((ids = dataFile.addVar("ContactElements_Nodes", ncIdxType, dv)).isNull())
485 throw FinleyException(msgPrefix+"add_var(ContactElements_Nodes)");
486 ids.putVar(m_contactElements->Nodes); // (num_ContactElements, num_ContactElements_numNodes) values written
487 }
488
489 // // // // // Points // // // // //
490 if (num_Points > 0) {
491 // Points_Id
492 if ((ids = dataFile.addVar("Points_Id", ncIdxType, ncdims[6])).isNull())
493 throw FinleyException(msgPrefix+"add_var(Points_Id)");
494 ids.putVar(m_points->Id); // num_Points
495
496 // Points_Tag
497 if ((ids = dataFile.addVar("Points_Tag", ncInt, ncdims[6])).isNull())
498 throw FinleyException(msgPrefix+"add_var(Points_Tag)");
499 ids.putVar(m_points->Tag); // num_Points
500
501 // Points_Owner
502 if ((ids = dataFile.addVar("Points_Owner", ncInt, ncdims[6])).isNull())
503 throw FinleyException(msgPrefix+"add_var(Points_Owner)");
504 ids.putVar(m_points->Owner); // num_Points
505
506 // Points_Color
507 if ((ids = dataFile.addVar("Points_Color", ncIdxType, ncdims[6])).isNull())
508 throw FinleyException(msgPrefix+"add_var(Points_Color)");
509 ids.putVar(m_points->Color); // num_Points
510
511 // Points_Nodes
512 if ((ids = dataFile.addVar("Points_Nodes", ncIdxType, ncdims[6])).isNull())
513 throw FinleyException(msgPrefix+"add_var(Points_Nodes)");
514 ids.putVar(m_points->Nodes); // num_Points)))
515
516 }
517
518 // // // // // TagMap // // // // //
519 if (num_Tags > 0) {
520 // Temp storage to gather node IDs
521 vector<int> Tags_keys;
522
523 // Copy tag data into temp arrays
524 TagMap::const_iterator it;
525 for (it = m_tagMap.begin(); it != m_tagMap.end(); it++) {
526 Tags_keys.push_back(it->second);
527 }
528
529 // Tags_keys
530 if ((ids = dataFile.addVar("Tags_keys", ncInt, ncdims[10])).isNull())
531 throw FinleyException(msgPrefix+"add_var(Tags_keys)");
532 ids.putVar(&Tags_keys[0]); // num_Tags
533
534 // Tags_names_*
535 // This is an array of strings, it should be stored as an array but
536 // instead I have hacked in one attribute per string because the NetCDF
537 // manual doesn't tell how to do an array of strings
538 int i = 0;
539 for (it = m_tagMap.begin(); it != m_tagMap.end(); it++, i++) {
540 stringstream ss;
541 ss << "Tags_name_" << i;
542 const string name(ss.str());
543 if (dataFile.putAtt(name.c_str(), it->first).isNull())
544 throw FinleyException(msgPrefix+"add_att(Tags_names_X)");
545 }
546 }
547
548 // Send token to next MPI process so he can take his turn
549 #ifdef ESYS_MPI
550 if (mpi_rank < mpi_size-1)
551 MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, getMPIComm());
552 #endif
553
554 // NetCDF file is closed by destructor of NcFile object
555
556 #else
557 throw FinleyException("FinleyDomain::dump: not configured with netCDF. "
558 "Please contact your installation manager.");
559 #endif // ESYS_HAVE_NETCDF
560 }
561
562 #else
563
564 void FinleyDomain::dump(const string& fileName) const
565 {
566 #ifdef ESYS_HAVE_NETCDF
567 const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
568 NcVar* ids;
569 index_t* index_ptr;
570 #ifdef ESYS_INDEXTYPE_LONG
571 NcType ncIdxType = ncLong;
572 #else
573 NcType ncIdxType = ncInt;
574 #endif
575 int num_Tags = 0;
576 int mpi_size = getMPISize();
577 int mpi_rank = getMPIRank();
578 int numDim = m_nodes->numDim;
579 dim_t numNodes = m_nodes->getNumNodes();
580 dim_t num_Elements = m_elements->numElements;
581 dim_t num_FaceElements = m_faceElements->numElements;
582 dim_t num_ContactElements = m_contactElements->numElements;
583 dim_t num_Points = m_points->numElements;
584 int num_Elements_numNodes = m_elements->numNodes;
585 int num_FaceElements_numNodes = m_faceElements->numNodes;
586 int num_ContactElements_numNodes = m_contactElements->numNodes;
587 #ifdef ESYS_MPI
588 MPI_Status status;
589 #endif
590
591 // Incoming token indicates it's my turn to write
592 #ifdef ESYS_MPI
593 if (mpi_rank > 0)
594 MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, getMPIComm(), &status);
595 #endif
596
597 const string newFileName(m_mpiInfo->appendRankToFileName(fileName));
598
599 // Figure out how much storage is required for tags
600 num_Tags = m_tagMap.size();
601
602 // NetCDF error handler
603 NcError err(NcError::verbose_nonfatal);
604 // Create the file
605 NcFile dataFile(newFileName.c_str(), NcFile::Replace);
606 string msgPrefix("Error in FinleyDomain::dump: NetCDF operation failed - ");
607 // check if writing was successful
608 if (!dataFile.is_valid())
609 throw FinleyException(msgPrefix + "Open file for output");
610
611 // Define dimensions (num_Elements and dim_Elements are identical,
612 // dim_Elements only appears if > 0)
613 if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
614 throw FinleyException(msgPrefix+"add_dim(numNodes)");
615 if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
616 throw FinleyException(msgPrefix+"add_dim(numDim)");
617 if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )
618 throw FinleyException(msgPrefix+"add_dim(mpi_size)");
619 if (num_Elements > 0)
620 if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
621 throw FinleyException(msgPrefix+"add_dim(dim_Elements)");
622 if (num_FaceElements > 0)
623 if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
624 throw FinleyException(msgPrefix+"add_dim(dim_FaceElements)");
625 if (num_ContactElements > 0)
626 if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
627 throw FinleyException(msgPrefix+"add_dim(dim_ContactElements)");
628 if (num_Points > 0)
629 if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
630 throw FinleyException(msgPrefix+"add_dim(dim_Points)");
631 if (num_Elements > 0)
632 if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
633 throw FinleyException(msgPrefix+"add_dim(dim_Elements_Nodes)");
634 if (num_FaceElements > 0)
635 if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
636 throw FinleyException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
637 if (num_ContactElements > 0)
638 if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
639 throw FinleyException(msgPrefix+"add_dim(dim_ContactElements_numNodes)");
640 if (num_Tags > 0)
641 if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
642 throw FinleyException(msgPrefix+"add_dim(dim_Tags)");
643
644 // Attributes: MPI size, MPI rank, Name, order, reduced_order
645 if (!dataFile.add_att("index_size", (int)sizeof(index_t)))
646 throw FinleyException(msgPrefix+"add_att(index_size)");
647 if (!dataFile.add_att("mpi_size", mpi_size))
648 throw FinleyException(msgPrefix+"add_att(mpi_size)");
649 if (!dataFile.add_att("mpi_rank", mpi_rank))
650 throw FinleyException(msgPrefix+"add_att(mpi_rank)");
651 if (!dataFile.add_att("Name", m_name.c_str()))
652 throw FinleyException(msgPrefix+"add_att(Name)");
653 if (!dataFile.add_att("numDim", numDim))
654 throw FinleyException(msgPrefix+"add_att(order)");
655 if (!dataFile.add_att("order", integrationOrder))
656 throw FinleyException(msgPrefix+"add_att(order)");
657 if (!dataFile.add_att("reduced_order", reducedIntegrationOrder))
658 throw FinleyException(msgPrefix+"add_att(reduced_order)");
659 if (!dataFile.add_att("numNodes", numNodes))
660 throw FinleyException(msgPrefix+"add_att(numNodes)");
661 if (!dataFile.add_att("num_Elements", num_Elements))
662 throw FinleyException(msgPrefix+"add_att(num_Elements)");
663 if (!dataFile.add_att("num_FaceElements", num_FaceElements))
664 throw FinleyException(msgPrefix+"add_att(num_FaceElements)");
665 if (!dataFile.add_att("num_ContactElements", num_ContactElements))
666 throw FinleyException(msgPrefix+"add_att(num_ContactElements)");
667 if (!dataFile.add_att("num_Points", num_Points))
668 throw FinleyException(msgPrefix+"add_att(num_Points)");
669 if (!dataFile.add_att("num_Elements_numNodes", num_Elements_numNodes))
670 throw FinleyException(msgPrefix+"add_att(num_Elements_numNodes)");
671 if (!dataFile.add_att("num_FaceElements_numNodes", num_FaceElements_numNodes))
672 throw FinleyException(msgPrefix+"add_att(num_FaceElements_numNodes)");
673 if (!dataFile.add_att("num_ContactElements_numNodes", num_ContactElements_numNodes))
674 throw FinleyException(msgPrefix+"add_att(num_ContactElements_numNodes)");
675 if (!dataFile.add_att("Elements_TypeId", m_elements->referenceElementSet->referenceElement->Type->TypeId) )
676 throw FinleyException(msgPrefix+"add_att(Elements_TypeId)");
677 if (!dataFile.add_att("FaceElements_TypeId", m_faceElements->referenceElementSet->referenceElement->Type->TypeId) )
678 throw FinleyException(msgPrefix+"add_att(FaceElements_TypeId)");
679 if (!dataFile.add_att("ContactElements_TypeId", m_contactElements->referenceElementSet->referenceElement->Type->TypeId) )
680 throw FinleyException(msgPrefix+"add_att(ContactElements_TypeId)");
681 if (!dataFile.add_att("Points_TypeId", m_points->referenceElementSet->referenceElement->Type->TypeId) )
682 throw FinleyException(msgPrefix+"add_att(Points_TypeId)");
683 if (!dataFile.add_att("num_Tags", num_Tags))
684 throw FinleyException(msgPrefix+"add_att(num_Tags)");
685
686 // // // // // Nodes // // // // //
687
688 // Nodes nodeDistribution
689 if (! (ids = dataFile.add_var("Nodes_NodeDistribution", ncIdxType, ncdims[2])) )
690 throw FinleyException(msgPrefix+"add_var(Nodes_NodeDistribution)");
691 index_ptr = &m_nodes->nodesDistribution->first_component[0];
692 if (! (ids->put(index_ptr, mpi_size+1)) )
693 throw FinleyException(msgPrefix+"put(Nodes_NodeDistribution)");
694
695 // Nodes degreesOfFreedomDistribution
696 if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncIdxType, ncdims[2])) )
697 throw FinleyException(msgPrefix+"add_var(Nodes_DofDistribution)");
698 index_ptr = &m_nodes->degreesOfFreedomDistribution->first_component[0];
699 if (! (ids->put(index_ptr, mpi_size+1)) )
700 throw FinleyException(msgPrefix+"put(Nodes_DofDistribution)");
701
702 // Only write nodes if non-empty because NetCDF doesn't like empty arrays
703 // (it treats them as NC_UNLIMITED)
704 if (numNodes > 0) {
705 // Nodes Id
706 if (! ( ids = dataFile.add_var("Nodes_Id", ncIdxType, ncdims[0])) )
707 throw FinleyException(msgPrefix+"add_var(Nodes_Id)");
708 if (! (ids->put(&m_nodes->Id[0], numNodes)) )
709 throw FinleyException(msgPrefix+"put(Nodes_Id)");
710
711 // Nodes Tag
712 if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
713 throw FinleyException(msgPrefix+"add_var(Nodes_Tag)");
714 if (! (ids->put(&m_nodes->Tag[0], numNodes)) )
715 throw FinleyException(msgPrefix+"put(Nodes_Tag)");
716
717 // Nodes gDOF
718 if (! ( ids = dataFile.add_var("Nodes_gDOF", ncIdxType, ncdims[0])) )
719 throw FinleyException(msgPrefix+"add_var(Nodes_gDOF)");
720 if (! (ids->put(&m_nodes->globalDegreesOfFreedom[0], numNodes)) )
721 throw FinleyException(msgPrefix+"put(Nodes_gDOF)");
722
723 // Nodes global node index
724 if (! ( ids = dataFile.add_var("Nodes_gNI", ncIdxType, ncdims[0])) )
725 throw FinleyException(msgPrefix+"add_var(Nodes_gNI)");
726 if (! (ids->put(&m_nodes->globalNodesIndex[0], numNodes)) )
727 throw FinleyException(msgPrefix+"put(Nodes_gNI)");
728
729 // Nodes grDof
730 if (! ( ids = dataFile.add_var("Nodes_grDfI", ncIdxType, ncdims[0])) )
731 throw FinleyException(msgPrefix+"add_var(Nodes_grDfI)");
732 if (! (ids->put(&m_nodes->globalReducedDOFIndex[0], numNodes)) )
733 throw FinleyException(msgPrefix+"put(Nodes_grDfI)");
734
735 // Nodes grNI
736 if (! ( ids = dataFile.add_var("Nodes_grNI", ncIdxType, ncdims[0])) )
737 throw FinleyException(msgPrefix+"add_var(Nodes_grNI)");
738 if (! (ids->put(&m_nodes->globalReducedNodesIndex[0], numNodes)) )
739 throw FinleyException(msgPrefix+"put(Nodes_grNI)");
740
741 // Nodes Coordinates
742 if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
743 throw FinleyException(msgPrefix+"add_var(Nodes_Coordinates)");
744 if (! (ids->put(m_nodes->Coordinates, numNodes, numDim)) )
745 throw FinleyException(msgPrefix+"put(Nodes_Coordinates)");
746 }
747
748 // // // // // Elements // // // // //
749 if (num_Elements > 0) {
750 // Elements_Id
751 if (! ( ids = dataFile.add_var("Elements_Id", ncIdxType, ncdims[3])) )
752 throw FinleyException(msgPrefix+"add_var(Elements_Id)");
753 if (! (ids->put(m_elements->Id, num_Elements)) )
754 throw FinleyException(msgPrefix+"put(Elements_Id)");
755
756 // Elements_Tag
757 if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
758 throw FinleyException(msgPrefix+"add_var(Elements_Tag)");
759 if (! (ids->put(m_elements->Tag, num_Elements)) )
760 throw FinleyException(msgPrefix+"put(Elements_Tag)");
761
762 // Elements_Owner
763 if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
764 throw FinleyException(msgPrefix+"add_var(Elements_Owner)");
765 if (! (ids->put(m_elements->Owner, num_Elements)) )
766 throw FinleyException(msgPrefix+"put(Elements_Owner)");
767
768 // Elements_Color
769 if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
770 throw FinleyException(msgPrefix+"add_var(Elements_Color)");
771 if (! (ids->put(m_elements->Color, num_Elements)) )
772 throw FinleyException(msgPrefix+"put(Elements_Color)");
773
774 // Elements_Nodes
775 if (! ( ids = dataFile.add_var("Elements_Nodes", ncIdxType, ncdims[3], ncdims[7]) ) )
776 throw FinleyException(msgPrefix+"add_var(Elements_Nodes)");
777 if (! (ids->put(&m_elements->Nodes[0], num_Elements, num_Elements_numNodes)) )
778 throw FinleyException(msgPrefix+"put(Elements_Nodes)");
779 }
780
781 // // // // // Face_Elements // // // // //
782 if (num_FaceElements > 0) {
783 // FaceElements_Id
784 if (!(ids = dataFile.add_var("FaceElements_Id", ncIdxType, ncdims[4])))
785 throw FinleyException(msgPrefix+"add_var(FaceElements_Id)");
786 if (!(ids->put(m_faceElements->Id, num_FaceElements)))
787 throw FinleyException(msgPrefix+"put(FaceElements_Id)");
788
789 // FaceElements_Tag
790 if (!(ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])))
791 throw FinleyException(msgPrefix+"add_var(FaceElements_Tag)");
792 if (!(ids->put(m_faceElements->Tag, num_FaceElements)))
793 throw FinleyException(msgPrefix+"put(FaceElements_Tag)");
794
795 // FaceElements_Owner
796 if (!(ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])))
797 throw FinleyException(msgPrefix+"add_var(FaceElements_Owner)");
798 if (!(ids->put(m_faceElements->Owner, num_FaceElements)))
799 throw FinleyException(msgPrefix+"put(FaceElements_Owner)");
800
801 // FaceElements_Color
802 if (!(ids = dataFile.add_var("FaceElements_Color", ncIdxType, ncdims[4])))
803 throw FinleyException(msgPrefix+"add_var(FaceElements_Color)");
804 if (!(ids->put(m_faceElements->Color, num_FaceElements)))
805 throw FinleyException(msgPrefix+"put(FaceElements_Color)");
806
807 // FaceElements_Nodes
808 if (!(ids = dataFile.add_var("FaceElements_Nodes", ncIdxType, ncdims[4], ncdims[8])))
809 throw FinleyException(msgPrefix+"add_var(FaceElements_Nodes)");
810 if (!(ids->put(m_faceElements->Nodes, num_FaceElements, num_FaceElements_numNodes)))
811 throw FinleyException(msgPrefix+"put(FaceElements_Nodes)");
812 }
813
814 // // // // // Contact_Elements // // // // //
815 if (num_ContactElements > 0) {
816 // ContactElements_Id
817 if (!(ids = dataFile.add_var("ContactElements_Id", ncIdxType, ncdims[5])))
818 throw FinleyException(msgPrefix+"add_var(ContactElements_Id)");
819 if (!(ids->put(m_contactElements->Id, num_ContactElements)))
820 throw FinleyException(msgPrefix+"put(ContactElements_Id)");
821
822 // ContactElements_Tag
823 if (!(ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])))
824 throw FinleyException(msgPrefix+"add_var(ContactElements_Tag)");
825 if (!(ids->put(m_contactElements->Tag, num_ContactElements)))
826 throw FinleyException(msgPrefix+"put(ContactElements_Tag)");
827
828 // ContactElements_Owner
829 if (!(ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])))
830 throw FinleyException(msgPrefix+"add_var(ContactElements_Owner)");
831 if (!(ids->put(m_contactElements->Owner, num_ContactElements)))
832 throw FinleyException(msgPrefix+"put(ContactElements_Owner)");
833
834 // ContactElements_Color
835 if (!(ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])))
836 throw FinleyException(msgPrefix+"add_var(ContactElements_Color)");
837 if (!(ids->put(m_contactElements->Color, num_ContactElements)))
838 throw FinleyException(msgPrefix+"put(ContactElements_Color)");
839
840 // ContactElements_Nodes
841 if (!(ids = dataFile.add_var("ContactElements_Nodes", ncIdxType, ncdims[5], ncdims[9])))
842 throw FinleyException(msgPrefix+"add_var(ContactElements_Nodes)");
843 if (!(ids->put(m_contactElements->Nodes, num_ContactElements, num_ContactElements_numNodes)))
844 throw FinleyException(msgPrefix+"put(ContactElements_Nodes)");
845 }
846
847 // // // // // Points // // // // //
848 if (num_Points > 0) {
849 // Points_Id
850 if (!(ids = dataFile.add_var("Points_Id", ncIdxType, ncdims[6])))
851 throw FinleyException(msgPrefix+"add_var(Points_Id)");
852 if (!(ids->put(m_points->Id, num_Points)))
853 throw FinleyException(msgPrefix+"put(Points_Id)");
854
855 // Points_Tag
856 if (!(ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])))
857 throw FinleyException(msgPrefix+"add_var(Points_Tag)");
858 if (!(ids->put(m_points->Tag, num_Points)))
859 throw FinleyException(msgPrefix+"put(Points_Tag)");
860
861 // Points_Owner
862 if (!(ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])))
863 throw FinleyException(msgPrefix+"add_var(Points_Owner)");
864 if (!(ids->put(m_points->Owner, num_Points)))
865 throw FinleyException(msgPrefix+"put(Points_Owner)");
866
867 // Points_Color
868 if (!(ids = dataFile.add_var("Points_Color", ncIdxType, ncdims[6])))
869 throw FinleyException(msgPrefix+"add_var(Points_Color)");
870 if (!(ids->put(m_points->Color, num_Points)))
871 throw FinleyException(msgPrefix+"put(Points_Color)");
872
873 // Points_Nodes
874 if (!(ids = dataFile.add_var("Points_Nodes", ncIdxType, ncdims[6])))
875 throw FinleyException(msgPrefix+"add_var(Points_Nodes)");
876 if (!(ids->put(m_points->Nodes, num_Points)))
877 throw FinleyException(msgPrefix+"put(Points_Nodes)");
878 }
879
880 // // // // // TagMap // // // // //
881 if (num_Tags > 0) {
882 // Temp storage to gather node IDs
883 vector<int> Tags_keys;
884
885 // Copy tag data into temp arrays
886 TagMap::const_iterator it;
887 for (it = m_tagMap.begin(); it != m_tagMap.end(); it++) {
888 Tags_keys.push_back(it->second);
889 }
890
891 // Tags_keys
892 if (!(ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])))
893 throw FinleyException(msgPrefix+"add_var(Tags_keys)");
894 if (!(ids->put(&Tags_keys[0], num_Tags)))
895 throw FinleyException(msgPrefix+"put(Tags_keys)");
896
897 // Tags_names_*
898 // This is an array of strings, it should be stored as an array but
899 // instead I have hacked in one attribute per string because the NetCDF
900 // manual doesn't tell how to do an array of strings
901 int i = 0;
902 for (it = m_tagMap.begin(); it != m_tagMap.end(); it++, i++) {
903 stringstream ss;
904 ss << "Tags_name_" << i;
905 const string name(ss.str());
906 if (!dataFile.add_att(name.c_str(), it->first.c_str()))
907 throw FinleyException(msgPrefix+"add_att(Tags_names_X)");
908 }
909 }
910
911 // Send token to next MPI process so he can take his turn
912 #ifdef ESYS_MPI
913 if (mpi_rank < mpi_size-1)
914 MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, getMPIComm());
915 #endif
916
917 // NetCDF file is closed by destructor of NcFile object
918
919 #else
920 throw FinleyException("FinleyDomain::dump: not configured with netCDF. "
921 "Please contact your installation manager.");
922 #endif // ESYS_HAVE_NETCDF
923 }
924
925 #endif
926
927 string FinleyDomain::getDescription() const
928 {
929 return "FinleyMesh";
930 }
931
932 string FinleyDomain::functionSpaceTypeAsString(int functionSpaceType) const
933 {
934 FunctionSpaceNamesMapType::iterator loc;
935 loc = m_functionSpaceTypeNames.find(functionSpaceType);
936 if (loc == m_functionSpaceTypeNames.end()) {
937 return "Invalid function space type code.";
938 } else {
939 return loc->second;
940 }
941 }
942
943 bool FinleyDomain::isValidFunctionSpaceType(int functionSpaceType) const
944 {
945 FunctionSpaceNamesMapType::iterator loc;
946 loc = m_functionSpaceTypeNames.find(functionSpaceType);
947 return (loc != m_functionSpaceTypeNames.end());
948 }
949
950 void FinleyDomain::setFunctionSpaceTypeNames()
951 {
952 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
953 DegreesOfFreedom,"Finley_DegreesOfFreedom [Solution(domain)]"));
954 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
955 ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom [ReducedSolution(domain)]"));
956 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
957 Nodes,"Finley_Nodes [ContinuousFunction(domain)]"));
958 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
959 ReducedNodes,"Finley_Reduced_Nodes [ReducedContinuousFunction(domain)]"));
960 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
961 Elements,"Finley_Elements [Function(domain)]"));
962 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
963 ReducedElements,"Finley_Reduced_Elements [ReducedFunction(domain)]"));
964 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
965 FaceElements,"Finley_Face_Elements [FunctionOnBoundary(domain)]"));
966 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
967 ReducedFaceElements,"Finley_Reduced_Face_Elements [ReducedFunctionOnBoundary(domain)]"));
968 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
969 Points,"Finley_Points [DiracDeltaFunctions(domain)]"));
970 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
971 ContactElementsZero,"Finley_Contact_Elements_0 [FunctionOnContactZero(domain)]"));
972 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
973 ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0 [ReducedFunctionOnContactZero(domain)]"));
974 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
975 ContactElementsOne,"Finley_Contact_Elements_1 [FunctionOnContactOne(domain)]"));
976 m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
977 ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1 [ReducedFunctionOnContactOne(domain)]"));
978 }
979
980 int FinleyDomain::getContinuousFunctionCode() const
981 {
982 return Nodes;
983 }
984
985 int FinleyDomain::getReducedContinuousFunctionCode() const
986 {
987 return ReducedNodes;
988 }
989
990 int FinleyDomain::getFunctionCode() const
991 {
992 return Elements;
993 }
994
995 int FinleyDomain::getReducedFunctionCode() const
996 {
997 return ReducedElements;
998 }
999
1000 int FinleyDomain::getFunctionOnBoundaryCode() const
1001 {
1002 return FaceElements;
1003 }
1004
1005 int FinleyDomain::getReducedFunctionOnBoundaryCode() const
1006 {
1007 return ReducedFaceElements;
1008 }
1009
1010 int FinleyDomain::getFunctionOnContactZeroCode() const
1011 {
1012 return ContactElementsZero;
1013 }
1014
1015 int FinleyDomain::getReducedFunctionOnContactZeroCode() const
1016 {
1017 return ReducedContactElementsZero;
1018 }
1019
1020 int FinleyDomain::getFunctionOnContactOneCode() const
1021 {
1022 return ContactElementsOne;
1023 }
1024
1025 int FinleyDomain::getReducedFunctionOnContactOneCode() const
1026 {
1027 return ReducedContactElementsOne;
1028 }
1029
1030 int FinleyDomain::getSolutionCode() const
1031 {
1032 return DegreesOfFreedom;
1033 }
1034
1035 int FinleyDomain::getReducedSolutionCode() const
1036 {
1037 return ReducedDegreesOfFreedom;
1038 }
1039
1040 int FinleyDomain::getDiracDeltaFunctionsCode() const
1041 {
1042 return Points;
1043 }
1044
1045 //
1046 // Return the number of data points summed across all MPI processes
1047 //
1048 dim_t FinleyDomain::getNumDataPointsGlobal() const
1049 {
1050 return m_nodes->getGlobalNumNodes();
1051 }
1052
1053 //
1054 // return the number of data points per sample and the number of samples
1055 // needed to represent data on a parts of the mesh.
1056 //
1057 pair<int,dim_t> FinleyDomain::getDataShape(int functionSpaceCode) const
1058 {
1059 int numDataPointsPerSample = 0;
1060 dim_t numSamples = 0;
1061 switch (functionSpaceCode) {
1062 case Nodes:
1063 numDataPointsPerSample = 1;
1064 numSamples = m_nodes->getNumNodes();
1065 break;
1066 case ReducedNodes:
1067 numDataPointsPerSample = 1;
1068 numSamples = m_nodes->getNumReducedNodes();
1069 break;
1070 case Elements:
1071 if (m_elements) {
1072 numSamples = m_elements->numElements;
1073 numDataPointsPerSample = m_elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
1074 }
1075 break;
1076 case ReducedElements:
1077 if (m_elements) {
1078 numSamples = m_elements->numElements;
1079 numDataPointsPerSample = m_elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
1080 }
1081 break;
1082 case FaceElements:
1083 if (m_faceElements) {
1084 numSamples = m_faceElements->numElements;
1085 numDataPointsPerSample = m_faceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
1086 }
1087 break;
1088 case ReducedFaceElements:
1089 if (m_faceElements) {
1090 numSamples = m_faceElements->numElements;
1091 numDataPointsPerSample = m_faceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
1092 }
1093 break;
1094 case Points:
1095 if (m_points) {
1096 numSamples = m_points->numElements;
1097 numDataPointsPerSample = 1;
1098 }
1099 break;
1100 case ContactElementsZero:
1101 case ContactElementsOne:
1102 if (m_contactElements) {
1103 numSamples = m_contactElements->numElements;
1104 numDataPointsPerSample = m_contactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
1105 }
1106 break;
1107 case ReducedContactElementsZero:
1108 case ReducedContactElementsOne:
1109 if (m_contactElements) {
1110 numSamples = m_contactElements->numElements;
1111 numDataPointsPerSample = m_contactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
1112 }
1113 break;
1114 case DegreesOfFreedom:
1115 if (m_nodes) {
1116 numSamples = m_nodes->getNumDegreesOfFreedom();
1117 numDataPointsPerSample = 1;
1118 }
1119 break;
1120 case ReducedDegreesOfFreedom:
1121 if (m_nodes) {
1122 numSamples = m_nodes->getNumReducedDegreesOfFreedom();
1123 numDataPointsPerSample = 1;
1124 }
1125 break;
1126 default:
1127 stringstream ss;
1128 ss << "Invalid function space type: " << functionSpaceCode
1129 << " for domain " << getDescription();
1130 throw ValueError(ss.str());
1131 }
1132 return pair<int,dim_t>(numDataPointsPerSample, numSamples);
1133 }
1134
1135 //
1136 // adds linear PDE of second order into a given stiffness matrix and
1137 // right hand side
1138 //
1139 void FinleyDomain::addPDEToSystem(
1140 escript::AbstractSystemMatrix& mat, escript::Data& rhs,
1141 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1142 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1143 const escript::Data& d, const escript::Data& y,
1144 const escript::Data& d_contact, const escript::Data& y_contact,
1145 const escript::Data& d_dirac, const escript::Data& y_dirac) const
1146 {
1147 #ifdef ESYS_HAVE_TRILINOS
1148 TrilinosMatrixAdapter* tm = dynamic_cast<TrilinosMatrixAdapter*>(&mat);
1149 if (tm) {
1150 tm->resumeFill();
1151 }
1152 #endif
1153
1154 Assemble_PDE(m_nodes, m_elements, mat.getPtr(), rhs, A, B, C, D, X, Y);
1155 Assemble_PDE(m_nodes, m_faceElements, mat.getPtr(), rhs,
1156 escript::Data(), escript::Data(), escript::Data(), d,
1157 escript::Data(), y);
1158 Assemble_PDE(m_nodes, m_contactElements, mat.getPtr(), rhs,
1159 escript::Data(), escript::Data(), escript::Data(), d_contact,
1160 escript::Data(), y_contact);
1161 Assemble_PDE(m_nodes, m_points, mat.getPtr(), rhs, escript::Data(),
1162 escript::Data(), escript::Data(), d_dirac,
1163 escript::Data(), y_dirac);
1164
1165 #ifdef ESYS_HAVE_TRILINOS
1166 if (tm) {
1167 tm->fillComplete(true);
1168 }
1169 #endif
1170 }
1171
1172 void FinleyDomain::addPDEToLumpedSystem(escript::Data& mat,
1173 const escript::Data& D,
1174 const escript::Data& d,
1175 const escript::Data& d_dirac,
1176 bool useHRZ) const
1177 {
1178 Assemble_LumpedSystem(m_nodes, m_elements, mat, D, useHRZ);
1179 Assemble_LumpedSystem(m_nodes, m_faceElements, mat, d, useHRZ);
1180 Assemble_LumpedSystem(m_nodes, m_points, mat, d_dirac, useHRZ);
1181 }
1182
1183 //
1184 // adds linear PDE of second order into the right hand side only
1185 //
1186 void FinleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
1187 const escript::Data& Y, const escript::Data& y,
1188 const escript::Data& y_contact, const escript::Data& y_dirac) const
1189 {
1190 Assemble_PDE(m_nodes, m_elements, escript::ASM_ptr(), rhs,
1191 escript::Data(), escript::Data(), escript::Data(),
1192 escript::Data(), X, Y);
1193
1194 Assemble_PDE(m_nodes, m_faceElements, escript::ASM_ptr(), rhs,
1195 escript::Data(), escript::Data(), escript::Data(),
1196 escript::Data(), escript::Data(), y);
1197
1198 Assemble_PDE(m_nodes, m_contactElements, escript::ASM_ptr(),
1199 rhs, escript::Data(), escript::Data(), escript::Data(),
1200 escript::Data(), escript::Data(), y_contact);
1201
1202 Assemble_PDE(m_nodes, m_points, escript::ASM_ptr(), rhs,
1203 escript::Data(), escript::Data(), escript::Data(),
1204 escript::Data(), escript::Data(), y_dirac);
1205 }
1206
1207 //
1208 // adds PDE of second order into a transport problem
1209 //
1210 void FinleyDomain::addPDEToTransportProblem(
1211 escript::AbstractTransportProblem& tp, escript::Data& source,
1212 const escript::Data& M, const escript::Data& A, const escript::Data& B,
1213 const escript::Data& C, const escript::Data& D, const escript::Data& X,
1214 const escript::Data& Y, const escript::Data& d, const escript::Data& y,
1215 const escript::Data& d_contact, const escript::Data& y_contact,
1216 const escript::Data& d_dirac, const escript::Data& y_dirac) const
1217 {
1218 #ifdef ESYS_HAVE_PASO
1219 paso::TransportProblem* ptp = dynamic_cast<paso::TransportProblem*>(&tp);
1220 if (!ptp)
1221 throw ValueError("Finley only supports Paso transport problems.");
1222
1223 source.expand();
1224
1225 escript::ASM_ptr mm(boost::static_pointer_cast<escript::AbstractSystemMatrix>(
1226 ptp->borrowMassMatrix()));
1227 escript::ASM_ptr tm(boost::static_pointer_cast<escript::AbstractSystemMatrix>(
1228 ptp->borrowTransportMatrix()));
1229
1230 Assemble_PDE(m_nodes, m_elements, mm, source, escript::Data(),
1231 escript::Data(), escript::Data(), M, escript::Data(),
1232 escript::Data());
1233 Assemble_PDE(m_nodes, m_elements, tm, source, A, B, C, D, X, Y);
1234 Assemble_PDE(m_nodes, m_faceElements, tm, source, escript::Data(),
1235 escript::Data(), escript::Data(), d, escript::Data(), y);
1236 Assemble_PDE(m_nodes, m_contactElements, tm, source,
1237 escript::Data(), escript::Data(), escript::Data(), d_contact,
1238 escript::Data(), y_contact);
1239 Assemble_PDE(m_nodes, m_points, tm, source, escript::Data(),
1240 escript::Data(), escript::Data(), d_dirac, escript::Data(),
1241 y_dirac);
1242 #else
1243 throw FinleyException("Transport problems require the Paso library which "
1244 "is not available.");
1245 #endif
1246 }
1247
1248 //
1249 // interpolates data between different function spaces
1250 //
1251 void FinleyDomain::interpolateOnDomain(escript::Data& target,
1252 const escript::Data& in) const
1253 {
1254 if (*in.getFunctionSpace().getDomain() != *this)
1255 throw ValueError("Illegal domain of interpolant.");
1256 if (*target.getFunctionSpace().getDomain() != *this)
1257 throw ValueError("Illegal domain of interpolation target.");
1258
1259 switch (in.getFunctionSpace().getTypeCode()) {
1260 case Nodes:
1261 switch (target.getFunctionSpace().getTypeCode()) {
1262 case Nodes:
1263 case ReducedNodes:
1264 case DegreesOfFreedom:
1265 case ReducedDegreesOfFreedom:
1266 if (in.isComplex())
1267 Assemble_CopyNodalData<cplx_t>(m_nodes, target, in);
1268 else
1269 Assemble_CopyNodalData<real_t>(m_nodes, target, in);
1270 break;
1271 case Elements:
1272 case ReducedElements:
1273 if (in.isComplex())
1274 Assemble_interpolate<cplx_t>(m_nodes, m_elements, in, target);
1275 else
1276 Assemble_interpolate<real_t>(m_nodes, m_elements, in, target);
1277 break;
1278 case FaceElements:
1279 case ReducedFaceElements:
1280 if (in.isComplex())
1281 Assemble_interpolate<cplx_t>(m_nodes, m_faceElements, in, target);
1282 else
1283 Assemble_interpolate<real_t>(m_nodes, m_faceElements, in, target);
1284 break;
1285 case Points:
1286 if (in.isComplex())
1287 Assemble_interpolate<cplx_t>(m_nodes, m_points, in, target);
1288 else
1289 Assemble_interpolate<real_t>(m_nodes, m_points, in, target);
1290 break;
1291 case ContactElementsZero:
1292 case ReducedContactElementsZero:
1293 case ContactElementsOne:
1294 case ReducedContactElementsOne:
1295 if (in.isComplex())
1296 Assemble_interpolate<cplx_t>(m_nodes, m_contactElements, in, target);
1297 else
1298 Assemble_interpolate<real_t>(m_nodes, m_contactElements, in, target);
1299 break;
1300 default:
1301 stringstream ss;
1302 ss << "interpolateOnDomain: Finley does not know anything "
1303 "about function space type "
1304 << target.getFunctionSpace().getTypeCode();
1305 throw ValueError(ss.str());
1306 }
1307 break;
1308 case ReducedNodes:
1309 switch(target.getFunctionSpace().getTypeCode()) {
1310 case Nodes:
1311 case ReducedNodes:
1312 case DegreesOfFreedom:
1313 case ReducedDegreesOfFreedom:
1314 if (in.isComplex())
1315 Assemble_CopyNodalData<cplx_t>(m_nodes, target, in);
1316 else
1317 Assemble_CopyNodalData<real_t>(m_nodes, target, in);
1318 break;
1319 case Elements:
1320 case ReducedElements:
1321 if (in.isComplex())
1322 Assemble_interpolate<cplx_t>(m_nodes, m_elements, in, target);
1323 else
1324 Assemble_interpolate<real_t>(m_nodes, m_elements, in, target);
1325 break;
1326 case FaceElements:
1327 case ReducedFaceElements:
1328 if (in.isComplex())
1329 Assemble_interpolate<cplx_t>(m_nodes, m_faceElements, in, target);
1330 else
1331 Assemble_interpolate<real_t>(m_nodes, m_faceElements, in, target);
1332 break;
1333 case Points:
1334 if (in.isComplex())
1335 Assemble_interpolate<cplx_t>(m_nodes, m_points, in, target);
1336 else
1337 Assemble_interpolate<real_t>(m_nodes, m_points, in, target);
1338 break;
1339 case ContactElementsZero:
1340 case ReducedContactElementsZero:
1341 case ContactElementsOne:
1342 case ReducedContactElementsOne:
1343 if (in.isComplex())
1344 Assemble_interpolate<cplx_t>(m_nodes, m_contactElements, in, target);
1345 else
1346 Assemble_interpolate<real_t>(m_nodes, m_contactElements, in, target);
1347 break;
1348 default:
1349 stringstream ss;
1350 ss << "interpolateOnDomain: Finley does not know anything "
1351 "about function space type "
1352 << target.getFunctionSpace().getTypeCode();
1353 throw ValueError(ss.str());
1354 }
1355 break;
1356 case Elements:
1357 if (target.getFunctionSpace().getTypeCode() == Elements) {
1358 if (in.isComplex())
1359 Assemble_CopyElementData<cplx_t>(m_elements, target, in);
1360 else
1361 Assemble_CopyElementData<real_t>(m_elements, target, in);
1362 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1363 if (in.isComplex())
1364 Assemble_AverageElementData<cplx_t>(m_elements, target, in);
1365 else
1366 Assemble_AverageElementData<real_t>(m_elements, target, in);
1367 } else {
1368 throw ValueError("No interpolation with data on elements possible.");
1369 }
1370 break;
1371 case ReducedElements:
1372 if (target.getFunctionSpace().getTypeCode() == ReducedElements) {
1373 if (in.isComplex())
1374 Assemble_CopyElementData<cplx_t>(m_elements, target, in);
1375 else
1376 Assemble_CopyElementData<real_t>(m_elements, target, in);
1377 } else if (target.getFunctionSpace().getTypeCode() == Elements) {
1378 if (in.isComplex())
1379 Assemble_CopyElementData<cplx_t>(m_elements, target, in);
1380 else
1381 Assemble_CopyElementData<real_t>(m_elements, target, in);
1382 } else {
1383 throw ValueError("No interpolation with data on elements "
1384 "with reduced integration order possible.");
1385 }
1386 break;
1387 case FaceElements:
1388 if (target.getFunctionSpace().getTypeCode() == FaceElements) {
1389 if (in.isComplex())
1390 Assemble_CopyElementData<cplx_t>(m_faceElements, target, in);
1391 else
1392 Assemble_CopyElementData<real_t>(m_faceElements, target, in);
1393 } else if (target.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1394 if (in.isComplex())
1395 Assemble_AverageElementData<cplx_t>(m_faceElements, target, in);
1396 else
1397 Assemble_AverageElementData<real_t>(m_faceElements, target, in);
1398 } else {
1399 throw ValueError("No interpolation with data on face elements possible.");
1400 }
1401 break;
1402 case ReducedFaceElements:
1403 if (target.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1404 if (in.isComplex())
1405 Assemble_CopyElementData<cplx_t>(m_faceElements, target, in);
1406 else
1407 Assemble_CopyElementData<real_t>(m_faceElements, target, in);
1408 } else {
1409 throw ValueError("No interpolation with data on face "
1410 "elements with reduced integration order possible.");
1411 }
1412 break;
1413 case Points:
1414 if (target.getFunctionSpace().getTypeCode() == Points) {
1415 if (in.isComplex())
1416 Assemble_CopyElementData<cplx_t>(m_points, target, in);
1417 else
1418 Assemble_CopyElementData<real_t>(m_points, target, in);
1419 } else {
1420 throw ValueError("No interpolation with data on points possible.");
1421 }
1422 break;
1423 case ContactElementsZero:
1424 case ContactElementsOne:
1425 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1426 if (in.isComplex())
1427 Assemble_CopyElementData<cplx_t>(m_contactElements, target, in);
1428 else
1429 Assemble_CopyElementData<real_t>(m_contactElements, target, in);
1430 } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1431 if (in.isComplex())
1432 Assemble_AverageElementData<cplx_t>(m_contactElements, target, in);
1433 else
1434 Assemble_AverageElementData<real_t>(m_contactElements, target, in);
1435 } else {
1436 throw ValueError("No interpolation with data on contact elements possible.");
1437 }
1438 break;
1439 case ReducedContactElementsZero:
1440 case ReducedContactElementsOne:
1441 if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1442 if (in.isComplex())
1443 Assemble_CopyElementData<cplx_t>(m_contactElements, target, in);
1444 else
1445 Assemble_CopyElementData<real_t>(m_contactElements, target, in);
1446 } else {
1447 throw ValueError("No interpolation with data on contact elements with reduced integration order possible.");
1448 }
1449 break;
1450 case DegreesOfFreedom:
1451 switch (target.getFunctionSpace().getTypeCode()) {
1452 case ReducedDegreesOfFreedom:
1453 case DegreesOfFreedom:
1454 if (in.isComplex())
1455 Assemble_CopyNodalData<cplx_t>(m_nodes, target, in);
1456 else
1457 Assemble_CopyNodalData<real_t>(m_nodes, target, in);
1458 break;
1459
1460 case Nodes:
1461 case ReducedNodes:
1462 if (getMPISize() > 1) {
1463 escript::Data temp(in);
1464 temp.expand();
1465 if (in.isComplex())
1466 Assemble_CopyNodalData<cplx_t>(m_nodes, target, temp);
1467 else
1468 Assemble_CopyNodalData<real_t>(m_nodes, target, temp);
1469 } else {
1470 if (in.isComplex())
1471 Assemble_CopyNodalData<cplx_t>(m_nodes, target, in);
1472 else
1473 Assemble_CopyNodalData<real_t>(m_nodes, target, in);
1474 }
1475 break;
1476 case Elements:
1477 case ReducedElements:
1478 if (getMPISize() > 1) {
1479 escript::Data temp(in, continuousFunction(*this));
1480 if (in.isComplex())
1481 Assemble_interpolate<cplx_t>(m_nodes, m_elements, temp, target);
1482 else
1483 Assemble_interpolate<real_t>(m_nodes, m_elements, temp, target);
1484 } else {
1485 if (in.isComplex())
1486 Assemble_interpolate<cplx_t>(m_nodes, m_elements, in, target);
1487 else
1488 Assemble_interpolate<real_t>(m_nodes, m_elements, in, target);
1489 }
1490 break;
1491 case FaceElements:
1492 case ReducedFaceElements:
1493 if (getMPISize() > 1) {
1494 escript::Data temp(in, continuousFunction(*this));
1495 if (in.isComplex())
1496 Assemble_interpolate<cplx_t>(m_nodes, m_faceElements, temp, target);
1497 else
1498 Assemble_interpolate<real_t>(m_nodes, m_faceElements, temp, target);
1499 } else {
1500 if (in.isComplex())
1501 Assemble_interpolate<cplx_t>(m_nodes, m_faceElements, in, target);
1502 else
1503 Assemble_interpolate<real_t>(m_nodes, m_faceElements, in, target);
1504 }
1505 break;
1506 case Points:
1507 if (getMPISize() > 1) {
1508 //escript::Data temp(in, continuousFunction(*this));
1509 } else {
1510 if (in.isComplex())
1511 Assemble_interpolate<cplx_t>(m_nodes, m_points, in, target);
1512 else
1513 Assemble_interpolate<real_t>(m_nodes, m_points, in, target);
1514 }
1515 break;
1516 case ContactElementsZero:
1517 case ContactElementsOne:
1518 case ReducedContactElementsZero:
1519 case ReducedContactElementsOne:
1520 if (getMPISize() > 1) {
1521 escript::Data temp(in, continuousFunction(*this));
1522 if (in.isComplex())
1523 Assemble_interpolate<cplx_t>(m_nodes, m_contactElements, temp, target);
1524 else
1525 Assemble_interpolate<real_t>(m_nodes, m_contactElements, temp, target);
1526 } else {
1527 if (in.isComplex())
1528 Assemble_interpolate<cplx_t>(m_nodes, m_contactElements, in, target);
1529 else
1530 Assemble_interpolate<real_t>(m_nodes, m_contactElements, in, target);
1531 }
1532 break;
1533 default:
1534 stringstream ss;
1535 ss << "interpolateOnDomain: Finley does not know anything "
1536 "about function space type "
1537 << target.getFunctionSpace().getTypeCode();
1538 throw ValueError(ss.str());
1539 }
1540 break;
1541 case ReducedDegreesOfFreedom:
1542 switch (target.getFunctionSpace().getTypeCode()) {
1543 case Nodes:
1544 throw ValueError("Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1545 case ReducedNodes:
1546 if (getMPISize() > 1) {
1547 escript::Data in2(in);
1548 in2.expand();
1549 if (in.isComplex())
1550 Assemble_CopyNodalData<cplx_t>(m_nodes, target, in2);
1551 else
1552 Assemble_CopyNodalData<real_t>(m_nodes, target, in2);
1553 } else {
1554 if (in.isComplex())
1555 Assemble_CopyNodalData<cplx_t>(m_nodes, target, in);
1556 else
1557 Assemble_CopyNodalData<real_t>(m_nodes, target, in);
1558 }
1559 break;
1560 case DegreesOfFreedom:
1561 throw ValueError("Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1562 case ReducedDegreesOfFreedom:
1563 if (in.isComplex())
1564 Assemble_CopyNodalData<cplx_t>(m_nodes, target, in);
1565 else
1566 Assemble_CopyNodalData<real_t>(m_nodes, target, in);
1567 break;
1568 case Elements:
1569 case ReducedElements:
1570 if (getMPISize() > 1) {
1571 escript::Data in2(in, reducedContinuousFunction(*this));
1572 if (in.isComplex())
1573 Assemble_interpolate<cplx_t>(m_nodes, m_elements, in2, target);
1574 else
1575 Assemble_interpolate<real_t>(m_nodes, m_elements, in2, target);
1576 } else {
1577 if (in.isComplex())
1578 Assemble_interpolate<cplx_t>(m_nodes, m_elements, in, target);
1579 else
1580 Assemble_interpolate<real_t>(m_nodes, m_elements, in, target);
1581 }
1582 break;
1583 case FaceElements:
1584 case ReducedFaceElements:
1585 if (getMPISize() > 1) {
1586 escript::Data in2(in, reducedContinuousFunction(*this));
1587 if (in.isComplex())
1588 Assemble_interpolate<cplx_t>(m_nodes, m_faceElements, in2, target);
1589 else
1590 Assemble_interpolate<real_t>(m_nodes, m_faceElements, in2, target);
1591 } else {
1592 if (in.isComplex())
1593 Assemble_interpolate<cplx_t>(m_nodes, m_faceElements, in, target);
1594 else
1595 Assemble_interpolate<real_t>(m_nodes, m_faceElements, in, target);
1596 }
1597 break;
1598 case Points:
1599 if (getMPISize() > 1) {
1600 escript::Data in2(in, reducedContinuousFunction(*this));
1601 if (in.isComplex())
1602 Assemble_interpolate<cplx_t>(m_nodes, m_points, in2, target);
1603 else
1604 Assemble_interpolate<real_t>(m_nodes, m_points, in2, target);
1605 } else {
1606 if (in.isComplex())
1607 Assemble_interpolate<cplx_t>(m_nodes, m_points, in, target);
1608 else
1609 Assemble_interpolate<real_t>(m_nodes, m_points, in, target);
1610 }
1611 break;
1612 case ContactElementsZero:
1613 case ContactElementsOne:
1614 case ReducedContactElementsZero:
1615 case ReducedContactElementsOne:
1616 if (getMPISize()>1) {
1617 escript::Data in2(in, reducedContinuousFunction(*this));
1618 if (in.isComplex())
1619 Assemble_interpolate<cplx_t>(m_nodes, m_contactElements, in2, target);
1620 else
1621 Assemble_interpolate<real_t>(m_nodes, m_contactElements, in2, target);
1622 } else {
1623 if (in.isComplex())
1624 Assemble_interpolate<cplx_t>(m_nodes, m_contactElements, in, target);
1625 else
1626 Assemble_interpolate<real_t>(m_nodes, m_contactElements, in, target);
1627 }
1628 break;
1629 default:
1630 stringstream ss;
1631 ss << "interpolateOnDomain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1632 throw ValueError(ss.str());
1633 }
1634 break;
1635 default:
1636 stringstream ss;
1637 ss << "interpolateOnDomain: Finley does not know anything about "
1638 "function space type " << in.getFunctionSpace().getTypeCode();
1639 throw ValueError(ss.str());
1640 }
1641 }
1642
1643 //
1644 // copies the locations of sample points into x
1645 //
1646 void FinleyDomain::setToX(escript::Data& arg) const
1647 {
1648 if (*arg.getFunctionSpace().getDomain() != *this)
1649 throw ValueError("setToX: Illegal domain of data point locations");
1650
1651 // in case of appropriate function space we can do the job directly:
1652 if (arg.getFunctionSpace().getTypeCode() == Nodes) {
1653 Assemble_NodeCoordinates(m_nodes, arg);
1654 } else {
1655 escript::Data tmp_data = Vector(0., continuousFunction(*this), true);
1656 Assemble_NodeCoordinates(m_nodes, tmp_data);
1657 // this is then interpolated onto arg:
1658 interpolateOnDomain(arg, tmp_data);
1659 }
1660 }
1661
1662 //
1663 // return the normal vectors at the location of data points as a Data object
1664 //
1665 void FinleyDomain::setToNormal(escript::Data& normal) const
1666 {
1667 if (*normal.getFunctionSpace().getDomain() != *this)
1668 throw ValueError("setToNormal: Illegal domain of normal locations");
1669
1670 if (normal.getFunctionSpace().getTypeCode() == FaceElements ||
1671 normal.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1672 Assemble_getNormal(m_nodes, m_faceElements, normal);
1673 } else if (normal.getFunctionSpace().getTypeCode() == ContactElementsOne ||
1674 normal.getFunctionSpace().getTypeCode() == ContactElementsZero ||
1675 normal.getFunctionSpace().getTypeCode() == ReducedContactElementsOne ||
1676 normal.getFunctionSpace().getTypeCode() == ReducedContactElementsZero) {
1677 Assemble_getNormal(m_nodes, m_contactElements, normal);
1678 } else {
1679 stringstream ss;
1680 ss << "setToNormal: Illegal function space type "
1681 << normal.getFunctionSpace().getTypeCode();
1682 throw ValueError(ss.str());
1683 }
1684 }
1685
1686 //
1687 // interpolates data to other domain
1688 //
1689 void FinleyDomain::interpolateAcross(escript::Data& /*target*/,
1690 const escript::Data& /*source*/) const
1691 {
1692 throw NotImplementedError("Finley does not allow interpolation across "
1693 "domains.");
1694 }
1695
1696 //
1697 // calculates the integral of a function defined on arg
1698 //
1699 void FinleyDomain::setToIntegrals(vector<real_t>& integrals,
1700 const escript::Data& arg) const
1701 {
1702 setToIntegralsWorker<real_t>(integrals, arg);
1703 }
1704
1705 void FinleyDomain::setToIntegrals(vector<cplx_t>& integrals,
1706 const escript::Data& arg) const
1707 {
1708 setToIntegralsWorker<cplx_t>(integrals, arg);
1709 }
1710
1711 template<typename Scalar>
1712 void FinleyDomain::setToIntegralsWorker(vector<Scalar>& integrals,
1713 const escript::Data& arg) const
1714 {
1715 if (*arg.getFunctionSpace().getDomain() != *this)
1716 throw ValueError("setToIntegrals: Illegal domain of integration kernel");
1717
1718 switch (arg.getFunctionSpace().getTypeCode()) {
1719 case Nodes:
1720 case ReducedNodes:
1721 case DegreesOfFreedom:
1722 case ReducedDegreesOfFreedom:
1723 {
1724 escript::Data temp(arg, escript::function(*this));
1725 Assemble_integrate(m_nodes, m_elements, temp, &integrals[0]);
1726 }
1727 break;
1728 case Points:
1729 case Elements:
1730 case ReducedElements:
1731 Assemble_integrate(m_nodes, m_elements, arg, &integrals[0]);
1732 break;
1733 case FaceElements:
1734 case ReducedFaceElements:
1735 Assemble_integrate(m_nodes, m_faceElements, arg, &integrals[0]);
1736 break;
1737 case ContactElementsZero:
1738 case ReducedContactElementsZero:
1739 case ContactElementsOne:
1740 case ReducedContactElementsOne:
1741 Assemble_integrate(m_nodes, m_contactElements, arg, &integrals[0]);
1742 break;
1743 default:
1744 stringstream ss;
1745 ss << "setToIntegrals: Finley does not know anything about "
1746 "function space type " << arg.getFunctionSpace().getTypeCode();
1747 throw ValueError(ss.str());
1748 }
1749 }
1750
1751 //
1752 // calculates the gradient of arg
1753 //
1754 void FinleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
1755 {
1756 if (*arg.getFunctionSpace().getDomain() != *this)
1757 throw ValueError("setToGradient: Illegal domain of gradient argument");
1758 if (*grad.getFunctionSpace().getDomain() != *this)
1759 throw ValueError("setToGradient: Illegal domain of gradient");
1760 if (grad.isComplex() != arg.isComplex())
1761 throw ValueError("setToGradient: Complexity of input and output must match");
1762
1763 escript::Data nodeData;
1764 if (getMPISize() > 1) {
1765 if (arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom) {
1766 nodeData = escript::Data(arg, continuousFunction(*this));
1767 } else if (arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom) {
1768 nodeData = escript::Data(arg, reducedContinuousFunction(*this));
1769 } else {
1770 nodeData = arg;
1771 }
1772 } else {
1773 nodeData = arg;
1774 }
1775 switch (grad.getFunctionSpace().getTypeCode()) {
1776 case Nodes:
1777 throw ValueError("Gradient at nodes is not supported.");
1778 case ReducedNodes:
1779 throw ValueError("Gradient at reduced nodes is not supported.");
1780 case Elements:
1781 case ReducedElements:
1782 if (arg.isComplex())
1783 Assemble_gradient<cplx_t>(m_nodes, m_elements, grad, nodeData);
1784 else
1785 Assemble_gradient<real_t>(m_nodes, m_elements, grad, nodeData);
1786 break;
1787 case FaceElements:
1788 case ReducedFaceElements:
1789 if (arg.isComplex())
1790 Assemble_gradient<cplx_t>(m_nodes, m_faceElements, grad, nodeData);
1791 else
1792 Assemble_gradient<real_t>(m_nodes, m_faceElements, grad, nodeData);
1793 break;
1794 case ContactElementsZero:
1795 case ReducedContactElementsZero:
1796 case ContactElementsOne:
1797 case ReducedContactElementsOne:
1798 if (arg.isComplex())
1799 Assemble_gradient<cplx_t>(m_nodes, m_contactElements, grad, nodeData);
1800 else
1801 Assemble_gradient<real_t>(m_nodes, m_contactElements, grad, nodeData);
1802 break;
1803 case Points:
1804 throw ValueError("Gradient at points is not supported.");
1805 case DegreesOfFreedom:
1806 throw ValueError("Gradient at degrees of freedom is not supported.");
1807 case ReducedDegreesOfFreedom:
1808 throw ValueError("Gradient at reduced degrees of freedom is not supported.");
1809 default:
1810 stringstream ss;
1811 ss << "Gradient: Finley does not know anything about function "
1812 "space type " << arg.getFunctionSpace().getTypeCode();
1813 throw ValueError(ss.str());
1814 }
1815 }
1816
1817 //
1818 // returns the size of elements
1819 //
1820 void FinleyDomain::setToSize(escript::Data& size) const
1821 {
1822 switch (size.getFunctionSpace().getTypeCode()) {
1823 case Nodes:
1824 throw ValueError("Size of nodes is not supported.");
1825 case ReducedNodes:
1826 throw ValueError("Size of reduced nodes is not supported.");
1827 case Elements:
1828 case ReducedElements:
1829 Assemble_getSize(m_nodes, m_elements, size);
1830 break;
1831 case FaceElements:
1832 case ReducedFaceElements:
1833 Assemble_getSize(m_nodes, m_faceElements, size);
1834 break;
1835 case ContactElementsZero:
1836 case ContactElementsOne:
1837 case ReducedContactElementsZero:
1838 case ReducedContactElementsOne:
1839 Assemble_getSize(m_nodes, m_contactElements, size);
1840 break;
1841 case Points:
1842 throw ValueError("Size of point elements is not supported.");
1843 case DegreesOfFreedom:
1844 throw ValueError("Size of degrees of freedom is not supported.");
1845 case ReducedDegreesOfFreedom:
1846 throw ValueError("Size of reduced degrees of freedom is not supported.");
1847 default:
1848 stringstream ss;
1849 ss << "setToSize: Finley does not know anything about function "
1850 "space type " << size.getFunctionSpace().getTypeCode();
1851 throw ValueError(ss.str());
1852 }
1853 }
1854
1855 //
1856 // sets the location of nodes
1857 //
1858 void FinleyDomain::setNewX(const escript::Data& newX)
1859 {
1860 if (*newX.getFunctionSpace().getDomain() != *this)
1861 throw ValueError("Illegal domain of new point locations");
1862
1863 if (newX.getFunctionSpace() == continuousFunction(*this)) {
1864 m_nodes->setCoordinates(newX);
1865 } else {
1866 throw ValueError("As of escript version 3.3 setNewX only accepts "
1867 "ContinuousFunction arguments. Please interpolate.");
1868 }
1869 }
1870
1871 bool FinleyDomain::ownSample(int fs_code, index_t id) const
1872 {
1873
1874 #ifdef ESYS_MPI
1875 if (getMPISize() > 1 && fs_code != FINLEY_DEGREES_OF_FREEDOM &&
1876 fs_code != FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
1877 if (fs_code == Nodes || fs_code == Elements ||
1878 fs_code == ReducedNodes || fs_code == ReducedElements ||
1879 fs_code == FaceElements || fs_code == ReducedFaceElements) {
1880 const index_t myFirstNode = m_nodes->getFirstNode();
1881 const index_t myLastNode = m_nodes->getLastNode();
1882 const index_t k = m_nodes->borrowGlobalNodesIndex()[id];
1883 return (myFirstNode <= k && k < myLastNode);
1884 } else {
1885 std::ostringstream oss;
1886 oss << "ownSample: unsupported function space type (" << fs_code <<")";
1887 throw ValueError(oss.str());
1888 }
1889 }
1890 #endif
1891 return true;
1892 }
1893
1894 //
1895 // creates a stiffness matrix and initializes it with zeros
1896 //
1897 escript::ASM_ptr FinleyDomain::newSystemMatrix(int row_blocksize,
1898 const escript::FunctionSpace& row_functionspace,
1899 int column_blocksize,
1900 const escript::FunctionSpace& column_functionspace,
1901 int type) const
1902 {
1903 // is the domain right?
1904 if (*row_functionspace.getDomain() != *this)
1905 throw ValueError("domain of row function space does not match the domain of matrix generator.");
1906 if (*column_functionspace.getDomain() != *this)
1907 throw ValueError("domain of column function space does not match the domain of matrix generator.");
1908
1909 bool reduceRowOrder = false;
1910 bool reduceColOrder = false;
1911 // is the function space type right?
1912 if (row_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1913 reduceRowOrder = true;
1914 } else if (row_functionspace.getTypeCode() != DegreesOfFreedom) {
1915 throw ValueError("illegal function space type for system matrix rows.");
1916 }
1917 if (column_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1918 reduceColOrder = true;
1919 } else if (column_functionspace.getTypeCode() != DegreesOfFreedom) {
1920 throw ValueError("illegal function space type for system matrix columns.");
1921 }
1922
1923 // generate matrix
1924 if (type & (int)SMT_TRILINOS) {
1925 #ifdef ESYS_HAVE_TRILINOS
1926 if (reduceRowOrder != reduceColOrder)
1927 throw ValueError("element order of matrix rows and columns must "
1928 "match when using Trilinos");
1929 const_TrilinosGraph_ptr graph(getTrilinosGraph(reduceRowOrder));
1930 bool isComplex = (type & (int)SMT_COMPLEX);
1931 bool unroll = (type & (int)SMT_UNROLL);
1932 escript::ASM_ptr sm(new TrilinosMatrixAdapter(m_mpiInfo, row_blocksize,
1933 row_functionspace, graph, isComplex, unroll));
1934 return sm;
1935 #else
1936 throw FinleyException("newSystemMatrix: finley was not compiled "
1937 "with Trilinos support so the Trilinos solver stack cannot be "
1938 "used.");
1939 #endif
1940 } else if (type & (int)SMT_PASO) {
1941 #ifdef ESYS_HAVE_PASO
1942 paso::SystemMatrixPattern_ptr pattern(getPasoPattern(
1943 reduceRowOrder, reduceColOrder));
1944 paso::SystemMatrix_ptr sm(new paso::SystemMatrix(type, pattern,
1945 row_blocksize, column_blocksize, false, row_functionspace,
1946 column_functionspace));
1947 return sm;
1948 #else
1949 throw FinleyException("newSystemMatrix: finley was not compiled "
1950 "with Paso support so the Paso solver stack cannot be used.");
1951 #endif
1952 } else {
1953 throw FinleyException("newSystemMatrix: unknown matrix type ID");
1954 }
1955 }
1956
1957 //
1958 // creates a TransportProblem
1959 //
1960 escript::ATP_ptr FinleyDomain::newTransportProblem(int blocksize,
1961 const escript::FunctionSpace& fs,
1962 int type) const
1963 {
1964 // is the domain right?
1965 if (*fs.getDomain() != *this)
1966 throw ValueError("domain of function space does not match the domain of transport problem generator.");
1967
1968 #ifdef ESYS_HAVE_PASO
1969 // is the function space type right
1970 bool reduceOrder = false;
1971 if (fs.getTypeCode() == ReducedDegreesOfFreedom) {
1972 reduceOrder = true;
1973 } else if (fs.getTypeCode() != DegreesOfFreedom) {
1974 throw ValueError("illegal function space type for transport problem.");
1975 }
1976
1977 // generate transport problem
1978 paso::SystemMatrixPattern_ptr pattern(getPasoPattern(
1979 reduceOrder, reduceOrder));
1980 paso::TransportProblem_ptr transportProblem(new paso::TransportProblem(
1981 pattern, blocksize, fs));
1982 return transportProblem;
1983 #else
1984 throw FinleyException("Transport problems require the Paso library which "
1985 "is not available.");
1986 #endif
1987 }
1988
1989 //
1990 // returns true if data on functionSpaceCode is considered as being cell centered
1991 //
1992 bool FinleyDomain::isCellOriented(int functionSpaceCode) const
1993 {
1994 switch (functionSpaceCode) {
1995 case Nodes:
1996 case DegreesOfFreedom:
1997 case ReducedDegreesOfFreedom:
1998 return false;
1999 case Elements:
2000 case FaceElements:
2001 case Points:
2002 case ContactElementsZero:
2003 case ContactElementsOne:
2004 case ReducedElements:
2005 case ReducedFaceElements:
2006 case ReducedContactElementsZero:
2007 case ReducedContactElementsOne:
2008 return true;
2009 }
2010 stringstream ss;
2011 ss << "isCellOriented: Finley does not know anything about "
2012 "function space type " << functionSpaceCode;
2013 throw ValueError(ss.str());
2014 }
2015
2016 bool
2017 FinleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
2018 {
2019 if (fs.empty())
2020 return false;
2021 // The idea is to use equivalence classes, i.e. types which can be
2022 // interpolated back and forth
2023 // class 1: DOF <-> Nodes
2024 // class 2: ReducedDOF <-> ReducedNodes
2025 // class 3: Points
2026 // class 4: Elements
2027 // class 5: ReducedElements
2028 // class 6: FaceElements
2029 // class 7: ReducedFaceElements
2030 // class 8: ContactElementZero <-> ContactElementOne
2031 // class 9: ReducedContactElementZero <-> ReducedContactElementOne
2032
2033 // There is also a set of lines. Interpolation is possible down a line but
2034 // not between lines.
2035 // class 1 and 2 belong to all lines so aren't considered.
2036 // line 0: class 3
2037 // line 1: class 4,5
2038 // line 2: class 6,7
2039 // line 3: class 8,9
2040
2041 // For classes with multiple members (e.g. class 2) we have vars to record
2042 // if there is at least one instance.
2043 // e.g. hasnodes is true if we have at least one instance of Nodes.
2044 vector<int> hasclass(10);
2045 vector<int> hasline(4);
2046 bool hasnodes = false;
2047 bool hasrednodes = false;
2048 bool hascez = false;
2049 bool hasrcez = false;
2050 for (int i = 0; i < fs.size(); ++i) {
2051 switch (fs[i]) {
2052 case Nodes:
2053 hasnodes = true; // fall through
2054 case DegreesOfFreedom:
2055 hasclass[1] = 1;
2056 break;
2057 case ReducedNodes:
2058 hasrednodes = true; // fall through
2059 case ReducedDegreesOfFreedom:
2060 hasclass[2] = 1;
2061 break;
2062 case Points:
2063 hasline[0] = 1;
2064 hasclass[3] = 1;
2065 break;
2066 case Elements:
2067 hasclass[4] = 1;
2068 hasline[1] = 1;
2069 break;
2070 case ReducedElements:
2071 hasclass[5] = 1;
2072 hasline[1] = 1;
2073 break;
2074 case FaceElements:
2075 hasclass[6] = 1;
2076 hasline[2] = 1;
2077 break;
2078 case ReducedFaceElements:
2079 hasclass[7] = 1;
2080 hasline[2] = 1;
2081 break;
2082 case ContactElementsZero:
2083 hascez = true; // fall through
2084 case ContactElementsOne:
2085 hasclass[8] = 1;
2086 hasline[3] = 1;
2087 break;
2088 case ReducedContactElementsZero:
2089 hasrcez = true; // fall through
2090 case ReducedContactElementsOne:
2091 hasclass[9] = 1;
2092 hasline[3] = 1;
2093 break;
2094 default:
2095 return false;
2096 }
2097 }
2098 int totlines = hasline[0]+hasline[1]+hasline[2]+hasline[3];
2099
2100 // fail if we have more than one leaf group
2101 if (totlines > 1)
2102 // there are at least two branches we can't interpolate between
2103 return false;
2104 else if (totlines == 1) {
2105 if (hasline[0] == 1) // we have points
2106 resultcode = Points;
2107 else if (hasline[1] == 1) {
2108 if (hasclass[5] == 1)
2109 resultcode=ReducedElements;
2110 else
2111 resultcode=Elements;
2112 } else if (hasline[2] == 1) {
2113 if (hasclass[7] == 1)
2114 resultcode=ReducedFaceElements;
2115 else
2116 resultcode=FaceElements;
2117 } else { // so we must be in line3
2118 if (hasclass[9] == 1) {
2119 // need something from class 9
2120 resultcode = (hasrcez ? ReducedContactElementsZero : ReducedContactElementsOne);
2121 } else {
2122 // something from class 8
2123 resultcode = (hascez ? ContactElementsZero : ContactElementsOne);
2124 }
2125 }
2126 } else { // totlines==0
2127 if (hasclass[2] == 1) {
2128 // something from class 2
2129 resultcode = (hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
2130 } else {
2131 // something from class 1
2132 resultcode = (hasnodes ? Nodes : DegreesOfFreedom);
2133 }
2134 }
2135 return true;
2136 }
2137
2138 bool FinleyDomain::probeInterpolationOnDomain(int functionSpaceType_source,
2139 int functionSpaceType_target) const
2140 {
2141 switch(functionSpaceType_source) {
2142 case Nodes:
2143 switch (functionSpaceType_target) {
2144 case Nodes:
2145 case ReducedNodes:
2146 case ReducedDegreesOfFreedom:
2147 case DegreesOfFreedom:
2148 case Elements:
2149 case ReducedElements:
2150 case FaceElements:
2151 case ReducedFaceElements:
2152 case Points:
2153 case ContactElementsZero:
2154 case ReducedContactElementsZero:
2155 case ContactElementsOne:
2156 case ReducedContactElementsOne:
2157 return true;
2158 default:
2159 stringstream ss;
2160 ss << "Interpolation On Domain: Finley does not know "
2161 "anything about function space type "
2162 << functionSpaceType_target;
2163 throw ValueError(ss.str());
2164 }
2165 case ReducedNodes:
2166 switch(functionSpaceType_target) {
2167 case ReducedNodes:
2168 case ReducedDegreesOfFreedom:
2169 case Elements:
2170 case ReducedElements:
2171 case FaceElements:
2172 case ReducedFaceElements:
2173 case Points:
2174 case ContactElementsZero:
2175 case ReducedContactElementsZero:
2176 case ContactElementsOne:
2177 case ReducedContactElementsOne:
2178 return true;
2179 case Nodes:
2180 case DegreesOfFreedom:
2181 return false;
2182 default:
2183 stringstream ss;
2184 ss << "Interpolation On Domain: Finley does not know "
2185 "anything about function space type "
2186 << functionSpaceType_target;
2187 throw ValueError(ss.str());
2188 }
2189 case Elements:
2190 return (functionSpaceType_target == Elements ||
2191 functionSpaceType_target == ReducedElements);
2192 case ReducedElements:
2193 return (functionSpaceType_target == ReducedElements);
2194 case FaceElements:
2195 return (functionSpaceType_target == FaceElements ||
2196 functionSpaceType_target == ReducedFaceElements);
2197 case ReducedFaceElements:
2198 return (functionSpaceType_target == ReducedFaceElements);
2199 case Points:
2200 return (functionSpaceType_target == Points);
2201 case ContactElementsZero:
2202 case ContactElementsOne:
2203 return (functionSpaceType_target == ContactElementsZero ||
2204 functionSpaceType_target == ContactElementsOne ||
2205 functionSpaceType_target == ReducedContactElementsZero ||
2206 functionSpaceType_target == ReducedContactElementsOne);
2207 case ReducedContactElementsZero:
2208 case ReducedContactElementsOne:
2209 return (functionSpaceType_target == ReducedContactElementsZero ||
2210 functionSpaceType_target == ReducedContactElementsOne);
2211 case DegreesOfFreedom:
2212 switch (functionSpaceType_target) {
2213 case ReducedDegreesOfFreedom:
2214 case DegreesOfFreedom:
2215 case Nodes:
2216 case ReducedNodes:
2217 case Elements:
2218 case ReducedElements:
2219 case Points:
2220 case FaceElements:
2221 case ReducedFaceElements:
2222 case ContactElementsZero:
2223 case ReducedContactElementsZero:
2224 case ContactElementsOne:
2225 case ReducedContactElementsOne:
2226 return true;
2227 default:
2228 stringstream ss;
2229 ss << "Interpolation On Domain: Finley does not know "
2230 "anything about function space type "
2231 << functionSpaceType_target;
2232 throw ValueError(ss.str());
2233 }
2234 case ReducedDegreesOfFreedom:
2235 switch(functionSpaceType_target) {
2236 case ReducedDegreesOfFreedom:
2237 case ReducedNodes:
2238 case Elements:
2239 case ReducedElements:
2240 case FaceElements:
2241 case ReducedFaceElements:
2242 case Points:
2243 case ContactElementsZero:
2244 case ReducedContactElementsZero:
2245 case ContactElementsOne:
2246 case ReducedContactElementsOne:
2247 return true;
2248 case Nodes:
2249 case DegreesOfFreedom:
2250 return false;
2251 default:
2252 stringstream ss;
2253 ss << "Interpolation On Domain: Finley does not know "
2254 "anything about function space type "
2255 << functionSpaceType_target;
2256 throw ValueError(ss.str());
2257 }
2258 }
2259 stringstream ss;
2260 ss << "Interpolation On Domain: Finley does not know anything "
2261 "about function space type " << functionSpaceType_source;
2262 throw ValueError(ss.str());
2263 }
2264
2265 signed char FinleyDomain::preferredInterpolationOnDomain(
2266 int functionSpaceType_source, int functionSpaceType_target) const
2267 {
2268 if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
2269 return 1;
2270 if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
2271 return -1;
2272
2273 return 0;
2274 }
2275
2276 bool FinleyDomain::probeInterpolationAcross(int /*source*/,
2277 const AbstractDomain& /*targetDomain*/, int /*target*/) const
2278 {
2279 return false;
2280 }
2281
2282 bool FinleyDomain::operator==(const AbstractDomain& other) const
2283 {
2284 const FinleyDomain* temp = dynamic_cast<const FinleyDomain*>(&other);
2285 if (temp) {
2286 return (m_nodes == temp->m_nodes &&
2287 m_elements == temp->m_elements &&
2288 m_faceElements == temp->m_faceElements &&
2289 m_contactElements == temp->m_contactElements &&
2290 m_points == temp->m_points);
2291 }
2292 return false;
2293 }
2294
2295 bool FinleyDomain::operator!=(const AbstractDomain& other) const
2296 {
2297 return !(operator==(other));
2298 }
2299
2300 int FinleyDomain::getSystemMatrixTypeId(const bp::object& options) const
2301 {
2302 const escript::SolverBuddy& sb = bp::extract<escript::SolverBuddy>(options);
2303
2304 int package = sb.getPackage();
2305 escript::SolverOptions method = sb.getSolverMethod();
2306 #ifdef ESYS_HAVE_TRILINOS
2307 bool isDirect = escript::isDirectSolver(method);
2308 #endif
2309
2310 // the configuration of finley should have taken care that we have either
2311 // paso or trilinos so here's how we prioritize
2312 #if defined(ESYS_HAVE_PASO) && defined(ESYS_HAVE_TRILINOS)
2313 // we have Paso & Trilinos so use Trilinos for parallel direct solvers and
2314 // for complex problems
2315 if (package == escript::SO_DEFAULT) {
2316 if ((method == escript::SO_METHOD_DIRECT && getMPISize() > 1)
2317 || isDirect
2318 || sb.isComplex()) {
2319 package = escript::SO_PACKAGE_TRILINOS;
2320 }
2321 }
2322 #endif
2323 #ifdef ESYS_HAVE_PASO
2324 if (package == escript::SO_DEFAULT)
2325 package = escript::SO_PACKAGE_PASO;
2326 #endif
2327 #ifdef ESYS_HAVE_TRILINOS
2328 if (package == escript::SO_DEFAULT)
2329 package = escript::SO_PACKAGE_TRILINOS;
2330 #endif
2331 if (package == escript::SO_PACKAGE_TRILINOS) {
2332 #ifdef ESYS_HAVE_TRILINOS
2333 int type = (int)SMT_TRILINOS;
2334 if (sb.isComplex())
2335 type |= (int)SMT_COMPLEX;
2336 // This is required because MueLu (AMG) and Amesos2 (direct) do not
2337 // support block matrices at this point. Remove if they ever do...
2338 if (sb.getPreconditioner() == escript::SO_PRECONDITIONER_AMG ||
2339 sb.getPreconditioner() == escript::SO_PRECONDITIONER_ILUT ||
2340 isDirect) {
2341 type |= (int)SMT_UNROLL;
2342 }
2343 return type;
2344 #else
2345 throw FinleyException("Trilinos requested but not built with Trilinos.");
2346 #endif
2347 }
2348 #ifdef ESYS_HAVE_PASO
2349 if (sb.isComplex()) {
2350 throw NotImplementedError("Paso does not support complex-valued matrices");
2351 }
2352 return (int)SMT_PASO | paso::SystemMatrix::getSystemMatrixTypeId(
2353 method, sb.getPreconditioner(), sb.getPackage(),
2354 sb.isSymmetric(), m_mpiInfo);
2355 #else
2356 throw FinleyException("Unable to find a working solver library!");
2357 #endif
2358 }
2359
2360 int FinleyDomain::getTransportTypeId(int solver, int preconditioner,
2361 int package, bool symmetry) const
2362 {
2363 #ifdef ESYS_HAVE_PASO
2364 return paso::TransportProblem::getTypeId(solver, preconditioner, package,
2365 symmetry, getMPI());
2366 #else
2367 throw FinleyException("Transport solvers require Paso but finley was not "
2368 "compiled with Paso!");
2369 #endif
2370 }
2371
2372 escript::Data FinleyDomain::getX() const
2373 {
2374 return continuousFunction(*this).getX();
2375 }
2376
2377 #ifdef ESYS_HAVE_BOOST_NUMPY
2378 boost::python::numpy::ndarray FinleyDomain::getNumpyX() const
2379 {
2380 return continuousFunction(*this).getNumpyX();
2381 }
2382
2383 boost::python::numpy::ndarray FinleyDomain::getConnectivityInfo() const
2384 {
2385 // Initialise boost numpy
2386 boost::python::numpy::initialize();
2387
2388 // Get the node information
2389 escript::DataTypes::index_t * nodedata = m_elements->Nodes;
2390
2391 // Work out how many elements there are
2392 int numberOfElements = m_elements->numElements;
2393
2394 // Work out how many data points there are per element
2395 int numNodesPerElement = m_elements->numNodes;
2396
2397 // Initialise the ndarray
2398 boost::python::tuple arrayshape = boost::python::make_tuple(numberOfElements, numNodesPerElement);
2399 boost::python::numpy::dtype datatype = boost::python::numpy::dtype::get_builtin<double>();
2400 boost::python::numpy::ndarray dataArray = boost::python::numpy::zeros(arrayshape, datatype);
2401
2402 // Initialise variables
2403 std::string localmsg;
2404 std::vector<const escript::DataTypes::real_t*> samplesR(1);
2405
2406 // Copy the information over
2407 // #pragma omp parallel for
2408 for (int k = 0; k < numberOfElements; k++) {
2409 for (int j = 0; j < numNodesPerElement; j++) {
2410 dataArray[k][j] = nodedata[j+k*numNodesPerElement];
2411 }
2412 }
2413
2414 // Print out the ndarray to the console - used during debugging
2415 // std::cout << "Finished array:\n" << bp::extract<char const *>(bp::str(dataArray)) << std::endl;
2416
2417 return dataArray;
2418 }
2419 #endif
2420
2421 int FinleyDomain::getVTKElementType() const
2422 {
2423 const_ReferenceElementSet_ptr refElement = m_elements->referenceElementSet;
2424 const_ReferenceElement_ptr borrowedRefElement = refElement->borrowReferenceElement(false);
2425 const ReferenceElementInfo* type = borrowedRefElement->Type;
2426
2427 // From vtkCellType.h
2428 // #define VTK_EMPTY_CELL 0
2429 // #define VTK_VERTEX 1
2430 // #define VTK_POLY_VERTEX 2
2431 // #define VTK_LINE 3
2432 // #define VTK_POLY_LINE 4
2433 // #define VTK_TRIANGLE 5
2434 // #define VTK_TRIANGLE_STRIP 6
2435 // #define VTK_POLYGON 7
2436 // #define VTK_PIXEL 8
2437 // #define VTK_QUAD 9
2438 // #define VTK_TETRA 10
2439 // #define VTK_VOXEL 11
2440 // #define VTK_HEXAHEDRON 12
2441 // #define VTK_WEDGE 13
2442 // #define VTK_PYRAMID 14
2443
2444 if(std::strcmp(type->Name, "Tri3")
2445 || std::strcmp(type->Name, "Tri6")
2446 || std::strcmp(type->Name, "Tri9")
2447 || std::strcmp(type->Name, "Tri10"))
2448 {
2449 return 5; //VTK_TRIANGLE
2450 }
2451 else if (std::strcmp(type->Name, "Rec4")
2452 || std::strcmp(type->Name, "Rec8")
2453 || std::strcmp(type->Name, "Rec9")
2454 || std::strcmp(type->Name, "Rec12")
2455 || std::strcmp(type->Name, "Rec16"))
2456 {
2457 return 8; // VTK_PIXEL
2458 }
2459 else if (std::strcmp(type->Name, "Tet4")
2460 || std::strcmp(type->Name, "Tet10")
2461 || std::strcmp(type->Name, "Tet16"))
2462 {
2463 return 10; // VTK_TETRA
2464 }
2465 else if (std::strcmp(type->Name, "Hex8")
2466 || std::strcmp(type->Name, "Hex20")
2467 || std::strcmp(type->Name, "Hex27")
2468 || std::strcmp(type->Name, "Hex32"))
2469 {
2470 return 11; // VTK_VOXEL
2471 }
2472 else
2473 {
2474 throw FinleyException("getVTKElementType: unknown element type");
2475 }
2476 }
2477
2478 escript::Data FinleyDomain::getNormal() const
2479 {
2480 return functionOnBoundary(*this).getNormal();
2481 }
2482
2483 escript::Data FinleyDomain::getSize() const
2484 {
2485 return escript::function(*this).getSize();
2486 }
2487
2488 const index_t* FinleyDomain::borrowSampleReferenceIDs(int functionSpaceType) const
2489 {
2490 index_t* out = NULL;
2491 switch (functionSpaceType) {
2492 case Nodes:
2493 out = m_nodes->Id;
2494 break;
2495 case ReducedNodes:
2496 out = m_nodes->reducedNodesId;
2497 break;
2498 case Elements:
2499 case ReducedElements:
2500 out = m_elements->Id;
2501 break;
2502 case FaceElements:
2503 case ReducedFaceElements:
2504 out = m_faceElements->Id;
2505 break;
2506 case Points:
2507 out = m_points->Id;
2508 break;
2509 case ContactElementsZero:
2510 case ReducedContactElementsZero:
2511 case ContactElementsOne:
2512 case ReducedContactElementsOne:
2513 out = m_contactElements->Id;
2514 break;
2515 case DegreesOfFreedom:
2516 out = m_nodes->degreesOfFreedomId;
2517 break;
2518 case ReducedDegreesOfFreedom:
2519 out = m_nodes->reducedDegreesOfFreedomId;
2520 break;
2521 default:
2522 stringstream ss;
2523 ss << "Invalid function space type: " << functionSpaceType
2524 << " for domain: " << getDescription();
2525 throw ValueError(ss.str());
2526 }
2527 return out;
2528 }
2529 int FinleyDomain::getTagFromSampleNo(int functionSpaceType, index_t sampleNo) const
2530 {
2531 int out = 0;
2532 switch (functionSpaceType) {
2533 case Nodes:
2534 out = m_nodes->Tag[sampleNo];
2535 break;
2536 case ReducedNodes:
2537 throw ValueError("ReducedNodes does not support tags.");
2538 case Elements:
2539 case ReducedElements:
2540 out = m_elements->Tag[sampleNo];
2541 break;
2542 case FaceElements:
2543 case ReducedFaceElements:
2544 out = m_faceElements->Tag[sampleNo];
2545 break;
2546 case Points:
2547 out = m_points->Tag[sampleNo];
2548 break;
2549 case ContactElementsZero:
2550 case ReducedContactElementsZero:
2551 case ContactElementsOne:
2552 case ReducedContactElementsOne:
2553 out = m_contactElements->Tag[sampleNo];
2554 break;
2555 case DegreesOfFreedom:
2556 throw ValueError("DegreesOfFreedom does not support tags.");
2557 case ReducedDegreesOfFreedom:
2558 throw ValueError("ReducedDegreesOfFreedom does not support tags.");
2559 default:
2560 stringstream ss;
2561 ss << "Invalid function space type: " << functionSpaceType
2562 << " for domain: " << getDescription();
2563 throw ValueError(ss.str());
2564 }
2565 return out;
2566 }
2567
2568
2569 void FinleyDomain::setTags(int functionSpaceType, int newTag, const escript::Data& mask) const
2570 {
2571 switch (functionSpaceType) {
2572 case Nodes:
2573 m_nodes->setTags(newTag, mask);
2574 break;
2575 case ReducedNodes:
2576 throw ValueError("ReducedNodes does not support tags");
2577 case DegreesOfFreedom:
2578 throw ValueError("DegreesOfFreedom does not support tags");
2579 case ReducedDegreesOfFreedom:
2580 throw ValueError("ReducedDegreesOfFreedom does not support tags");
2581 case Elements:
2582 case ReducedElements:
2583 m_elements->setTags(newTag, mask);
2584 break;
2585 case FaceElements:
2586 case ReducedFaceElements:
2587 m_faceElements->setTags(newTag, mask);
2588 break;
2589 case Points:
2590 m_points->setTags(newTag, mask);
2591 break;
2592 case ContactElementsZero:
2593 case ReducedContactElementsZero:
2594 case ContactElementsOne:
2595 case ReducedContactElementsOne:
2596 m_contactElements->setTags(newTag, mask);
2597 break;
2598 default:
2599 stringstream ss;
2600 ss << "Finley does not know anything about function space type "
2601 << functionSpaceType;
2602 throw ValueError(ss.str());
2603 }
2604 }
2605
2606 void FinleyDomain::setTagMap(const string& name, int tag)
2607 {
2608 m_tagMap[name] = tag;
2609 }
2610
2611 int FinleyDomain::getTag(const string& name) const
2612 {
2613 TagMap::const_iterator it = m_tagMap.find(name);
2614 if (it == m_tagMap.end()) {
2615 stringstream ss;
2616 ss << "getTag: unknown tag name " << name << ".";
2617 throw escript::ValueError(ss.str());
2618 }
2619 return it->second;
2620 }
2621
2622 bool FinleyDomain::isValidTagName(const string& name) const
2623 {
2624 return (m_tagMap.count(name) > 0);
2625 }
2626
2627 string FinleyDomain::showTagNames() const
2628 {
2629 stringstream ss;
2630 TagMap::const_iterator it = m_tagMap.begin();
2631 while (it != m_tagMap.end()) {
2632 ss << it->first;
2633 ++it;
2634 if (it != m_tagMap.end())
2635 ss << ", ";
2636 }
2637 return ss.str();
2638 }
2639
2640 int FinleyDomain::getNumberOfTagsInUse(int functionSpaceCode) const
2641 {
2642 switch (functionSpaceCode) {
2643 case Nodes:
2644 return m_nodes->tagsInUse.size();
2645 case ReducedNodes:
2646 throw ValueError("ReducedNodes does not support tags");
2647 case DegreesOfFreedom:
2648 throw ValueError("DegreesOfFreedom does not support tags");
2649 case ReducedDegreesOfFreedom:
2650 throw ValueError("ReducedDegreesOfFreedom does not support tags");
2651 case Elements:
2652 case ReducedElements:
2653 return m_elements->tagsInUse.size();
2654 case FaceElements:
2655 case ReducedFaceElements:
2656 return m_faceElements->tagsInUse.size();
2657 case Points:
2658 return m_points->tagsInUse.size();
2659 case ContactElementsZero:
2660 case ReducedContactElementsZero:
2661 case ContactElementsOne:
2662 case ReducedContactElementsOne:
2663 return m_contactElements->tagsInUse.size();
2664 }
2665 stringstream ss;
2666 ss << "Finley does not know anything about function space type "
2667 << functionSpaceCode;
2668 throw ValueError(ss.str());
2669 }
2670
2671 const int* FinleyDomain::borrowListOfTagsInUse(int functionSpaceCode) const
2672 {
2673 switch (functionSpaceCode) {
2674 case Nodes:
2675 if (m_nodes->tagsInUse.empty())
2676 return NULL;
2677 else
2678 return &m_nodes->tagsInUse[0];
2679 case ReducedNodes:
2680 throw ValueError("ReducedNodes does not support tags");
2681 case DegreesOfFreedom:
2682 throw ValueError("DegreesOfFreedom does not support tags");
2683 case ReducedDegreesOfFreedom:
2684 throw ValueError("ReducedDegreesOfFreedom does not support tags");
2685 case Elements:
2686 case ReducedElements:
2687 if (m_elements->tagsInUse.empty())
2688 return NULL;
2689 else
2690 return &m_elements->tagsInUse[0];
2691 case FaceElements:
2692 case ReducedFaceElements:
2693 if (m_faceElements->tagsInUse.empty())
2694 return NULL;
2695 else
2696 return &m_faceElements->tagsInUse[0];
2697 case Points:
2698 if (m_points->tagsInUse.empty())
2699 return NULL;
2700 else
2701 return &m_points->tagsInUse[0];
2702 case ContactElementsZero:
2703 case ReducedContactElementsZero:
2704 case ContactElementsOne:
2705 case ReducedContactElementsOne:
2706 if (m_contactElements->tagsInUse.empty())
2707 return NULL;
2708 else
2709 return &m_contactElements->tagsInUse[0];
2710 }
2711 stringstream ss;
2712 ss << "Finley does not know anything about function space type "
2713 << functionSpaceCode;
2714 throw ValueError(ss.str());
2715 }
2716
2717 bool FinleyDomain::canTag(int functionSpaceCode) const
2718 {
2719 switch(functionSpaceCode) {
2720 case Nodes:
2721 case Elements:
2722 case ReducedElements:
2723 case FaceElements:
2724 case ReducedFaceElements:
2725 case Points:
2726 case ContactElementsZero:
2727 case ReducedContactElementsZero:
2728 case ContactElementsOne:
2729 case ReducedContactElementsOne:
2730 return true;
2731 default:
2732 return false;
2733 }
2734 }
2735
2736 FinleyDomain::StatusType FinleyDomain::getStatus() const
2737 {
2738 return m_nodes->status;
2739 }
2740
2741 int FinleyDomain::getApproximationOrder(int functionSpaceCode) const
2742 {
2743 switch(functionSpaceCode) {
2744 case Nodes:
2745 case DegreesOfFreedom:
2746 return approximationOrder;
2747 case ReducedNodes:
2748 case ReducedDegreesOfFreedom:
2749 return reducedApproximationOrder;
2750 case Elements:
2751 case FaceElements:
2752 case Points:
2753 case ContactElementsZero:
2754 case ContactElementsOne:
2755 return integrationOrder;
2756 case ReducedElements:
2757 case ReducedFaceElements:
2758 case ReducedContactElementsZero:
2759 case ReducedContactElementsOne:
2760 return reducedIntegrationOrder;
2761 }
2762 stringstream ss;
2763 ss << "Finley does not know anything about function space type "
2764 << functionSpaceCode;
2765 throw ValueError(ss.str());
2766 }
2767
2768 escript::Data FinleyDomain::randomFill(
2769 const escript::DataTypes::ShapeType& shape,
2770 const escript::FunctionSpace& what, long seed,
2771 const bp::tuple& filter) const
2772 {
2773 escript::Data towipe(0, shape, what, true);
2774 // since we just made this object, no sharing is possible and we don't
2775 // need to check for exclusive write
2776 escript::DataTypes::RealVectorType& dv(towipe.getExpandedVectorReference());
2777 escript::randomFillArray(seed, &dv[0], dv.size());
2778 return towipe;
2779 }
2780
2781 /// prepares the mesh for further use
2782 void FinleyDomain::prepare(bool optimize)
2783 {
2784 setOrders();
2785
2786 // first step is to distribute the elements according to a global
2787 // distribution of DOF
2788 IndexVector distribution(m_mpiInfo->size + 1);
2789
2790 // first we create dense labeling for the DOFs
2791 dim_t newGlobalNumDOFs = m_nodes->createDenseDOFLabeling();
2792
2793 // create a distribution of the global DOFs and determine the MPI rank
2794 // controlling the DOFs on this processor
2795 m_mpiInfo->setDistribution(0, newGlobalNumDOFs - 1, &distribution[0]);
2796
2797 // now the mesh is re-distributed according to the distribution vector
2798 // this will redistribute the Nodes and Elements including overlap and
2799 // will create an element colouring but will not create any mappings
2800 // (see later in this function)
2801 distributeByRankOfDOF(distribution);
2802
2803 // at this stage we are able to start an optimization of the DOF
2804 // distribution using ParaMetis. On return distribution is altered and
2805 // new DOF IDs have been assigned
2806 if (optimize && m_mpiInfo->size > 1) {
2807 optimizeDOFDistribution(distribution);
2808 distributeByRankOfDOF(distribution);
2809 }
2810 // the local labelling of the degrees of freedom is optimized
2811 if (optimize) {
2812 optimizeDOFLabeling(distribution);
2813 }
2814
2815 // rearrange elements with the aim of bringing elements closer to memory
2816 // locations of the nodes (distributed shared memory!):
2817 optimizeElementOrdering();
2818
2819 // create the global indices
2820 std::vector<short> maskReducedNodes(m_nodes->getNumNodes(), -1);
2821 IndexVector nodeDistribution(m_mpiInfo->size + 1);
2822 markNodes(maskReducedNodes, 0, true);
2823 IndexVector indexReducedNodes = util::packMask(maskReducedNodes);
2824
2825 m_nodes->createDenseNodeLabeling(nodeDistribution, distribution);
2826 // created reduced DOF labeling
2827 m_nodes->createDenseReducedLabeling(maskReducedNodes, false);
2828 // created reduced node labeling
2829 m_nodes->createDenseReducedLabeling(maskReducedNodes, true);
2830 // create the missing mappings
2831 m_nodes->createNodeMappings(indexReducedNodes, distribution, nodeDistribution);
2832
2833 updateTagList();
2834 }
2835
2836 /// redistributes the Nodes and Elements including overlap
2837 /// according to the DOF distribution. It will create an element colouring
2838 /// but will not create any mappings.
2839 void FinleyDomain::distributeByRankOfDOF(const std::vector<index_t>& dof_distribution)
2840 {
2841 std::vector<int> mpiRankOfDOF(m_nodes->getNumNodes());
2842 m_nodes->assignMPIRankToDOFs(mpiRankOfDOF, dof_distribution);
2843
2844 // first, the elements are redistributed according to mpiRankOfDOF
2845 // at the input the Node tables refer to the local labeling of the nodes
2846 // while at the output they refer to the global labeling which is rectified
2847 // in the next step
2848 m_elements->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->Id);
2849 m_faceElements->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->Id);
2850 m_contactElements->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->Id);
2851 m_points->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->Id);
2852
2853 resolveNodeIds();
2854
2855 // create a local labeling of the DOFs
2856 const std::pair<index_t,index_t> dof_range(m_nodes->getDOFRange());
2857 const index_t len = dof_range.second-dof_range.first+1;
2858 // local mask for used nodes
2859 std::vector<index_t> localDOF_mask(len, -1);
2860 std::vector<index_t> localDOF_map(m_nodes->getNumNodes(), -1);
2861
2862 #pragma omp parallel for
2863 for (index_t n = 0; n < m_nodes->getNumNodes(); n++) {
2864 #ifdef BOUNDS_CHECK
2865 ESYS_ASSERT(m_nodes->globalDegreesOfFreedom[n]-dof_range.first < len, "BOUNDS_CHECK");
2866 ESYS_ASSERT(m_nodes->globalDegreesOfFreedom[n]-dof_range.first >= 0, "BOUNDS_CHECK");
2867 #endif
2868 localDOF_mask[m_nodes->globalDegreesOfFreedom[n]-dof_range.first] = n;
2869 }
2870
2871 index_t numDOFs = 0;
2872 for (index_t n = 0; n < len; n++) {
2873 const index_t k = localDOF_mask[n];
2874 if (k >= 0) {
2875 localDOF_mask[n] = numDOFs;
2876 numDOFs++;
2877 }
2878 }
2879 #pragma omp parallel for
2880 for (index_t n = 0; n < m_nodes->getNumNodes(); n++) {
2881 const index_t k = localDOF_mask[m_nodes->globalDegreesOfFreedom[n]-dof_range.first];
2882 localDOF_map[n] = k;
2883 }
2884 // create element coloring
2885 createColoring(localDOF_map);
2886 }
2887
2888 /// optimizes the labeling of the DOFs on each processor
2889 void FinleyDomain::optimizeDOFLabeling(const IndexVector& distribution)
2890 {
2891 // this method relies on Pattern::reduceBandwidth so requires PASO
2892 // at the moment
2893 #ifdef ESYS_HAVE_PASO
2894 const int myRank = getMPIRank();
2895 const int mpiSize = getMPISize();
2896 const index_t myFirstVertex = distribution[myRank];
2897 const index_t myLastVertex = distribution[myRank+1];
2898 const dim_t myNumVertices = myLastVertex-myFirstVertex;
2899 dim_t len = 0;
2900 for (int p = 0; p < mpiSize; ++p)
2901 len=std::max(len, distribution[p+1]-distribution[p]);
2902
2903 boost::scoped_array<IndexList> index_list(new IndexList[myNumVertices]);
2904 boost::scoped_array<index_t> newGlobalDOFID(new index_t[len]);
2905
2906 // create the adjacency structure xadj and adjncy
2907 #pragma omp parallel
2908 {
2909 // insert contributions from element matrices into columns index
2910 IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list.get(),
2911 myFirstVertex, myLastVertex, m_elements,
2912 m_nodes->globalDegreesOfFreedom,
2913 m_nodes->globalDegreesOfFreedom);
2914 IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list.get(),
2915 myFirstVertex, myLastVertex, m_faceElements,
2916 m_nodes->globalDegreesOfFreedom,
2917 m_nodes->globalDegreesOfFreedom);
2918 IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list.get(),
2919 myFirstVertex, myLastVertex, m_contactElements,
2920 m_nodes->globalDegreesOfFreedom,
2921 m_nodes->globalDegreesOfFreedom);
2922 IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list.get(),
2923 myFirstVertex, myLastVertex, m_points,
2924 m_nodes->globalDegreesOfFreedom,
2925 m_nodes->globalDegreesOfFreedom);
2926 }
2927 // create the local matrix pattern
2928 paso::Pattern_ptr pattern = paso::Pattern::fromIndexListArray(0,
2929 myNumVertices, index_list.get(), myFirstVertex, myLastVertex,
2930 -myFirstVertex);
2931
2932 pattern->reduceBandwidth(&newGlobalDOFID[0]);
2933
2934 // shift new labeling to create a global id