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

Contents of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4867 - (show annotations)
Fri Apr 11 12:29:49 2014 UTC (5 years, 4 months ago) by caltinay
File size: 91264 byte(s)
paso: starting to polish

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26