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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4498 - (show annotations)
Tue Jul 16 01:56:37 2013 UTC (5 years, 9 months ago) by caltinay
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 90622 byte(s)
Another special case fix.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26