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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6058 - (show annotations)
Thu Mar 10 06:51:55 2016 UTC (3 years, 1 month ago) by caltinay
File size: 90211 byte(s)
added getPtr() to AbstractSystemMatrix so we can now use shared systemmatrix
pointers rather than circumventing them.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26