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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6317 - (show annotations)
Thu Jun 23 04:48:29 2016 UTC (2 years, 3 months ago) by caltinay
File size: 97309 byte(s)
parts for feature #3 - automatically choose trilinos if:
1) package == default
2) method == direct
3) trilinos available
4) mpi size > 1

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26