/[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 4612 - (show annotations)
Thu Jan 9 05:30:27 2014 UTC (5 years, 3 months ago) by gross
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 90623 byte(s)
modifications to sonic seismic
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
758 Assemble_PDE(mesh->Nodes, mesh->FaceElements, S, rhs,
759 escript::Data(), escript::Data(), escript::Data(), d,
760 escript::Data(), y);
761 checkFinleyError();
762
763 Assemble_PDE(mesh->Nodes, mesh->ContactElements, S, rhs,
764 escript::Data(), escript::Data(), escript::Data(), d_contact,
765 escript::Data(), y_contact);
766 checkFinleyError();
767
768 Assemble_PDE(mesh->Nodes, mesh->Points, S, rhs, escript::Data(),
769 escript::Data(), escript::Data(), d_dirac, escript::Data(), y_dirac);
770 checkFinleyError();
771 }
772
773 void MeshAdapter::addPDEToLumpedSystem(escript::Data& mat,
774 const escript::Data& D, const escript::Data& d,
775 const escript::Data& d_dirac, const bool useHRZ) const
776 {
777 Mesh* mesh=m_finleyMesh.get();
778 Assemble_LumpedSystem(mesh->Nodes, mesh->Elements, mat, D, useHRZ);
779 checkFinleyError();
780
781 Assemble_LumpedSystem(mesh->Nodes, mesh->FaceElements, mat, d, useHRZ);
782 checkFinleyError();
783
784 Assemble_LumpedSystem(mesh->Nodes, mesh->Points, mat, d_dirac, useHRZ);
785 checkFinleyError();
786 }
787
788 //
789 // adds linear PDE of second order into the right hand side only
790 //
791 void MeshAdapter::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
792 const escript::Data& Y, const escript::Data& y,
793 const escript::Data& y_contact, const escript::Data& y_dirac) const
794 {
795 Mesh* mesh=m_finleyMesh.get();
796 Assemble_PDE(mesh->Nodes, mesh->Elements, 0, rhs,
797 escript::Data(), escript::Data(), escript::Data(), escript::Data(),
798 X, Y);
799 checkFinleyError();
800
801 Assemble_PDE(mesh->Nodes, mesh->FaceElements, 0, rhs,
802 escript::Data(), escript::Data(), escript::Data(), escript::Data(),
803 escript::Data(), y);
804 checkFinleyError();
805
806 Assemble_PDE(mesh->Nodes, mesh->ContactElements, 0, rhs,
807 escript::Data(), escript::Data(), escript::Data(),
808 escript::Data(), escript::Data(), y_contact);
809 checkFinleyError();
810
811 Assemble_PDE(mesh->Nodes, mesh->Points, 0, rhs,
812 escript::Data(), escript::Data(), escript::Data(), escript::Data(),
813 escript::Data(), y_dirac);
814 checkFinleyError();
815 }
816
817 //
818 // adds PDE of second order into a transport problem
819 //
820 void MeshAdapter::addPDEToTransportProblem(
821 escript::AbstractTransportProblem& tp, escript::Data& source,
822 const escript::Data& M, const escript::Data& A, const escript::Data& B,
823 const escript::Data& C, const escript::Data& D, const escript::Data& X,
824 const escript::Data& Y, const escript::Data& d, const escript::Data& y,
825 const escript::Data& d_contact, const escript::Data& y_contact,
826 const escript::Data& d_dirac, const escript::Data& y_dirac) const
827 {
828 TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
829 if (!tpa)
830 throw FinleyAdapterException("finley only supports Paso transport problems.");
831
832 source.expand();
833
834 Mesh* mesh=m_finleyMesh.get();
835 Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
836
837 Assemble_PDE(mesh->Nodes, mesh->Elements, _tp->mass_matrix, source,
838 escript::Data(), escript::Data(), escript::Data(),
839 M, escript::Data(), escript::Data());
840 checkFinleyError();
841
842 Assemble_PDE(mesh->Nodes, mesh->Elements, _tp->transport_matrix,
843 source, A, B, C, D, X, Y);
844 checkFinleyError();
845
846 Assemble_PDE(mesh->Nodes, mesh->FaceElements, _tp->transport_matrix,
847 source, escript::Data(), escript::Data(),
848 escript::Data(), d, escript::Data(), y);
849 checkFinleyError();
850
851 Assemble_PDE(mesh->Nodes, mesh->ContactElements,
852 _tp->transport_matrix, source, escript::Data(),
853 escript::Data(), escript::Data(), d_contact,
854 escript::Data(), y_contact);
855 checkFinleyError();
856
857 Assemble_PDE(mesh->Nodes, mesh->Points, _tp->transport_matrix,
858 source, escript::Data(), escript::Data(),
859 escript::Data(), d_dirac, escript::Data(), y_dirac);
860 checkFinleyError();
861 }
862
863 //
864 // interpolates data between different function spaces
865 //
866 void MeshAdapter::interpolateOnDomain(escript::Data& target, const escript::Data& in) const
867 {
868 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
869 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
870 if (inDomain!=*this)
871 throw FinleyAdapterException("Error - Illegal domain of interpolant.");
872 if (targetDomain!=*this)
873 throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
874
875 Mesh* mesh=m_finleyMesh.get();
876 switch(in.getFunctionSpace().getTypeCode()) {
877 case Nodes:
878 switch(target.getFunctionSpace().getTypeCode()) {
879 case Nodes:
880 case ReducedNodes:
881 case DegreesOfFreedom:
882 case ReducedDegreesOfFreedom:
883 Assemble_CopyNodalData(mesh->Nodes, target, in);
884 break;
885 case Elements:
886 case ReducedElements:
887 Assemble_interpolate(mesh->Nodes, mesh->Elements, in,target);
888 break;
889 case FaceElements:
890 case ReducedFaceElements:
891 Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
892 break;
893 case Points:
894 Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
895 break;
896 case ContactElementsZero:
897 case ReducedContactElementsZero:
898 case ContactElementsOne:
899 case ReducedContactElementsOne:
900 Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in, target);
901 break;
902 default:
903 stringstream temp;
904 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
905 throw FinleyAdapterException(temp.str());
906 break;
907 }
908 break;
909 case ReducedNodes:
910 switch(target.getFunctionSpace().getTypeCode()) {
911 case Nodes:
912 case ReducedNodes:
913 case DegreesOfFreedom:
914 case ReducedDegreesOfFreedom:
915 Assemble_CopyNodalData(mesh->Nodes, target, in);
916 break;
917 case Elements:
918 case ReducedElements:
919 Assemble_interpolate(mesh->Nodes, mesh->Elements, in, target);
920 break;
921 case FaceElements:
922 case ReducedFaceElements:
923 Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
924 break;
925 case Points:
926 Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
927 break;
928 case ContactElementsZero:
929 case ReducedContactElementsZero:
930 case ContactElementsOne:
931 case ReducedContactElementsOne:
932 Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in, target);
933 break;
934 default:
935 stringstream temp;
936 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
937 throw FinleyAdapterException(temp.str());
938 break;
939 }
940 break;
941 case Elements:
942 if (target.getFunctionSpace().getTypeCode()==Elements) {
943 Assemble_CopyElementData(mesh->Elements, target, in);
944 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
945 Assemble_AverageElementData(mesh->Elements, target, in);
946 } else {
947 throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
948 }
949 break;
950 case ReducedElements:
951 if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
952 Assemble_CopyElementData(mesh->Elements, target, in);
953 } else {
954 throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
955 }
956 break;
957 case FaceElements:
958 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
959 Assemble_CopyElementData(mesh->FaceElements, target, in);
960 } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
961 Assemble_AverageElementData(mesh->FaceElements, target, in);
962 } else {
963 throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
964 }
965 break;
966 case ReducedFaceElements:
967 if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
968 Assemble_CopyElementData(mesh->FaceElements, target, in);
969 } else {
970 throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
971 }
972 break;
973 case Points:
974 if (target.getFunctionSpace().getTypeCode()==Points) {
975 Assemble_CopyElementData(mesh->Points, target, in);
976 } else {
977 throw FinleyAdapterException("Error - No interpolation with data on points possible.");
978 }
979 break;
980 case ContactElementsZero:
981 case ContactElementsOne:
982 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
983 Assemble_CopyElementData(mesh->ContactElements, target, in);
984 } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
985 Assemble_AverageElementData(mesh->ContactElements, target, in);
986 } else {
987 throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
988 }
989 break;
990 case ReducedContactElementsZero:
991 case ReducedContactElementsOne:
992 if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
993 Assemble_CopyElementData(mesh->ContactElements, target, in);
994 } else {
995 throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
996 }
997 break;
998 case DegreesOfFreedom:
999 switch(target.getFunctionSpace().getTypeCode()) {
1000 case ReducedDegreesOfFreedom:
1001 case DegreesOfFreedom:
1002 Assemble_CopyNodalData(mesh->Nodes, target, in);
1003 break;
1004
1005 case Nodes:
1006 case ReducedNodes:
1007 if (getMPISize()>1) {
1008 escript::Data in2=escript::Data(in);
1009 in2.expand();
1010 Assemble_CopyNodalData(mesh->Nodes, target, in2);
1011 } else {
1012 Assemble_CopyNodalData(mesh->Nodes, target, in);
1013 }
1014 break;
1015 case Elements:
1016 case ReducedElements:
1017 if (getMPISize()>1) {
1018 escript::Data in2=escript::Data(in, continuousFunction(*this));
1019 Assemble_interpolate(mesh->Nodes, mesh->Elements, in2, target);
1020 } else {
1021 Assemble_interpolate(mesh->Nodes, mesh->Elements, in, target);
1022 }
1023 break;
1024 case FaceElements:
1025 case ReducedFaceElements:
1026 if (getMPISize()>1) {
1027 escript::Data in2=escript::Data(in, continuousFunction(*this));
1028 Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in2, target);
1029 } else {
1030 Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
1031 }
1032 break;
1033 case Points:
1034 if (getMPISize()>1) {
1035 //escript::Data in2=escript::Data(in, continuousFunction(*this) );
1036 } else {
1037 Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
1038 }
1039 break;
1040 case ContactElementsZero:
1041 case ContactElementsOne:
1042 case ReducedContactElementsZero:
1043 case ReducedContactElementsOne:
1044 if (getMPISize()>1) {
1045 escript::Data in2=escript::Data(in, continuousFunction(*this));
1046 Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in2, target);
1047 } else {
1048 Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in, target);
1049 }
1050 break;
1051 default:
1052 stringstream temp;
1053 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1054 throw FinleyAdapterException(temp.str());
1055 }
1056 break;
1057 case ReducedDegreesOfFreedom:
1058 switch(target.getFunctionSpace().getTypeCode()) {
1059 case Nodes:
1060 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1061 case ReducedNodes:
1062 if (getMPISize()>1) {
1063 escript::Data in2=escript::Data(in);
1064 in2.expand();
1065 Assemble_CopyNodalData(mesh->Nodes, target, in2);
1066 } else {
1067 Assemble_CopyNodalData(mesh->Nodes, target, in);
1068 }
1069 break;
1070 case DegreesOfFreedom:
1071 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1072 break;
1073 case ReducedDegreesOfFreedom:
1074 Assemble_CopyNodalData(mesh->Nodes, target, in);
1075 break;
1076 case Elements:
1077 case ReducedElements:
1078 if (getMPISize()>1) {
1079 escript::Data in2=escript::Data(in, reducedContinuousFunction(*this) );
1080 Assemble_interpolate(mesh->Nodes, mesh->Elements, in2, target);
1081 } else {
1082 Assemble_interpolate(mesh->Nodes, mesh->Elements, in, target);
1083 }
1084 break;
1085 case FaceElements:
1086 case ReducedFaceElements:
1087 if (getMPISize()>1) {
1088 escript::Data in2=escript::Data(in, reducedContinuousFunction(*this) );
1089 Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in2, target);
1090 } else {
1091 Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
1092 }
1093 break;
1094 case Points:
1095 if (getMPISize()>1) {
1096 escript::Data in2=escript::Data(in, reducedContinuousFunction(*this));
1097 Assemble_interpolate(mesh->Nodes, mesh->Points, in2, target);
1098 } else {
1099 Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
1100 }
1101 break;
1102 case ContactElementsZero:
1103 case ContactElementsOne:
1104 case ReducedContactElementsZero:
1105 case ReducedContactElementsOne:
1106 if (getMPISize()>1) {
1107 escript::Data in2=escript::Data(in, reducedContinuousFunction(*this));
1108 Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in2, target);
1109 } else {
1110 Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in, target);
1111 }
1112 break;
1113 default:
1114 stringstream temp;
1115 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1116 throw FinleyAdapterException(temp.str());
1117 break;
1118 }
1119 break;
1120 default:
1121 stringstream temp;
1122 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1123 throw FinleyAdapterException(temp.str());
1124 break;
1125 }
1126 checkFinleyError();
1127 }
1128
1129 //
1130 // copies the locations of sample points into x
1131 //
1132 void MeshAdapter::setToX(escript::Data& arg) const
1133 {
1134 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1135 if (argDomain!=*this)
1136 throw FinleyAdapterException("Error - Illegal domain of data point locations");
1137 Mesh* mesh=m_finleyMesh.get();
1138 // in case of appropriate function space we can do the job directly:
1139 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1140 Assemble_NodeCoordinates(mesh->Nodes, arg);
1141 } else {
1142 escript::Data tmp_data=Vector(0., continuousFunction(*this), true);
1143 Assemble_NodeCoordinates(mesh->Nodes, tmp_data);
1144 // this is then interpolated onto arg:
1145 interpolateOnDomain(arg, tmp_data);
1146 }
1147 checkFinleyError();
1148 }
1149
1150 //
1151 // return the normal vectors at the location of data points as a Data object
1152 //
1153 void MeshAdapter::setToNormal(escript::Data& normal) const
1154 {
1155 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1156 if (normalDomain!=*this)
1157 throw FinleyAdapterException("Error - Illegal domain of normal locations");
1158 Mesh* mesh=m_finleyMesh.get();
1159 switch(normal.getFunctionSpace().getTypeCode()) {
1160 case Nodes:
1161 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1162 break;
1163 case ReducedNodes:
1164 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1165 break;
1166 case Elements:
1167 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1168 break;
1169 case ReducedElements:
1170 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1171 break;
1172 case FaceElements:
1173 case ReducedFaceElements:
1174 Assemble_getNormal(mesh->Nodes, mesh->FaceElements, normal);
1175 break;
1176 case Points:
1177 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1178 break;
1179 case ContactElementsOne:
1180 case ContactElementsZero:
1181 case ReducedContactElementsOne:
1182 case ReducedContactElementsZero:
1183 Assemble_getNormal(mesh->Nodes, mesh->ContactElements, normal);
1184 break;
1185 case DegreesOfFreedom:
1186 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1187 break;
1188 case ReducedDegreesOfFreedom:
1189 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1190 break;
1191 default:
1192 stringstream temp;
1193 temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1194 throw FinleyAdapterException(temp.str());
1195 break;
1196 }
1197 checkFinleyError();
1198 }
1199
1200 //
1201 // interpolates data to other domain
1202 //
1203 void MeshAdapter::interpolateACross(escript::Data& target, const escript::Data& source) const
1204 {
1205 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains.");
1206 }
1207
1208 //
1209 // calculates the integral of a function defined on arg
1210 //
1211 void MeshAdapter::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
1212 {
1213 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1214 if (argDomain!=*this)
1215 throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1216
1217 double blocktimer_start = blocktimer_time();
1218 Mesh* mesh=m_finleyMesh.get();
1219 switch(arg.getFunctionSpace().getTypeCode()) {
1220 case Nodes:
1221 {
1222 escript::Data temp(arg, escript::function(*this));
1223 Assemble_integrate(mesh->Nodes, mesh->Elements, temp, &integrals[0]);
1224 }
1225 break;
1226 case ReducedNodes:
1227 {
1228 escript::Data temp(arg, escript::function(*this));
1229 Assemble_integrate(mesh->Nodes, mesh->Elements, temp, &integrals[0]);
1230 }
1231 break;
1232 case Elements:
1233 case ReducedElements:
1234 Assemble_integrate(mesh->Nodes, mesh->Elements, arg, &integrals[0]);
1235 break;
1236 case FaceElements:
1237 case ReducedFaceElements:
1238 Assemble_integrate(mesh->Nodes, mesh->FaceElements, arg, &integrals[0]);
1239 break;
1240 case Points:
1241 throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1242 break;
1243 case ContactElementsZero:
1244 case ReducedContactElementsZero:
1245 case ContactElementsOne:
1246 case ReducedContactElementsOne:
1247 Assemble_integrate(mesh->Nodes, mesh->ContactElements, arg, &integrals[0]);
1248 break;
1249 case DegreesOfFreedom:
1250 case ReducedDegreesOfFreedom:
1251 {
1252 escript::Data temp(arg, escript::function(*this));
1253 Assemble_integrate(mesh->Nodes, mesh->Elements, temp, &integrals[0]);
1254 }
1255 break;
1256 default:
1257 stringstream temp;
1258 temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1259 throw FinleyAdapterException(temp.str());
1260 break;
1261 }
1262 checkFinleyError();
1263 blocktimer_increment("integrate()", blocktimer_start);
1264 }
1265
1266 //
1267 // calculates the gradient of arg
1268 //
1269 void MeshAdapter::setToGradient(escript::Data& grad, const escript::Data& arg) const
1270 {
1271 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1272 if (argDomain!=*this)
1273 throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1274 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1275 if (gradDomain!=*this)
1276 throw FinleyAdapterException("Error - Illegal domain of gradient");
1277
1278 Mesh* mesh=m_finleyMesh.get();
1279 escript::Data nodeData;
1280 if (getMPISize()>1) {
1281 if (arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom) {
1282 nodeData=escript::Data(arg, continuousFunction(*this));
1283 } else if(arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom) {
1284 nodeData=escript::Data(arg, reducedContinuousFunction(*this));
1285 } else {
1286 nodeData = arg;
1287 }
1288 } else {
1289 nodeData = arg;
1290 }
1291 switch(grad.getFunctionSpace().getTypeCode()) {
1292 case Nodes:
1293 throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1294 break;
1295 case ReducedNodes:
1296 throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1297 break;
1298 case Elements:
1299 case ReducedElements:
1300 Assemble_gradient(mesh->Nodes, mesh->Elements, grad, nodeData);
1301 break;
1302 case FaceElements:
1303 case ReducedFaceElements:
1304 Assemble_gradient(mesh->Nodes, mesh->FaceElements, grad, nodeData);
1305 break;
1306 case Points:
1307 throw FinleyAdapterException("Error - Gradient at points is not supported.");
1308 break;
1309 case ContactElementsZero:
1310 case ReducedContactElementsZero:
1311 case ContactElementsOne:
1312 case ReducedContactElementsOne:
1313 Assemble_gradient(mesh->Nodes, mesh->ContactElements, grad, nodeData);
1314 break;
1315 case DegreesOfFreedom:
1316 throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1317 break;
1318 case ReducedDegreesOfFreedom:
1319 throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1320 break;
1321 default:
1322 stringstream temp;
1323 temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1324 throw FinleyAdapterException(temp.str());
1325 break;
1326 }
1327 checkFinleyError();
1328 }
1329
1330 //
1331 // returns the size of elements
1332 //
1333 void MeshAdapter::setToSize(escript::Data& size) const
1334 {
1335 Mesh* mesh=m_finleyMesh.get();
1336 switch(size.getFunctionSpace().getTypeCode()) {
1337 case Nodes:
1338 throw FinleyAdapterException("Error - Size of nodes is not supported.");
1339 break;
1340 case ReducedNodes:
1341 throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1342 break;
1343 case Elements:
1344 case ReducedElements:
1345 Assemble_getSize(mesh->Nodes, mesh->Elements, size);
1346 break;
1347 case FaceElements:
1348 case ReducedFaceElements:
1349 Assemble_getSize(mesh->Nodes, mesh->FaceElements, size);
1350 break;
1351 case Points:
1352 throw FinleyAdapterException("Error - Size of point elements is not supported.");
1353 break;
1354 case ContactElementsZero:
1355 case ContactElementsOne:
1356 case ReducedContactElementsZero:
1357 case ReducedContactElementsOne:
1358 Assemble_getSize(mesh->Nodes,mesh->ContactElements,size);
1359 break;
1360 case DegreesOfFreedom:
1361 throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1362 break;
1363 case ReducedDegreesOfFreedom:
1364 throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1365 break;
1366 default:
1367 stringstream temp;
1368 temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1369 throw FinleyAdapterException(temp.str());
1370 break;
1371 }
1372 checkFinleyError();
1373 }
1374
1375 //
1376 // sets the location of nodes
1377 //
1378 void MeshAdapter::setNewX(const escript::Data& new_x)
1379 {
1380 Mesh* mesh=m_finleyMesh.get();
1381 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1382 if (newDomain!=*this)
1383 throw FinleyAdapterException("Error - Illegal domain of new point locations");
1384 if (new_x.getFunctionSpace() == continuousFunction(*this)) {
1385 mesh->setCoordinates(new_x);
1386 } else {
1387 throw FinleyAdapterException("As of escript version 3.3 SetX() only accepts ContinuousFunction arguments. Please interpolate.");
1388 }
1389 checkFinleyError();
1390 }
1391
1392 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1393 {
1394 if (getMPISize() > 1) {
1395 #ifdef ESYS_MPI
1396 int myFirstNode=0, myLastNode=0, k=0;
1397 int* globalNodeIndex=0;
1398 Mesh* mesh_p=m_finleyMesh.get();
1399 /*
1400 * this method is only used by saveDataCSV which would use the returned
1401 * values for reduced nodes wrongly so this case is disabled for now
1402 if (fs_code == FINLEY_REDUCED_NODES) {
1403 myFirstNode = NodeFile_getFirstReducedNode(mesh_p->Nodes);
1404 myLastNode = NodeFile_getLastReducedNode(mesh_p->Nodes);
1405 globalNodeIndex = NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1406 } else
1407 */
1408 if (fs_code == FINLEY_NODES) {
1409 myFirstNode = mesh_p->Nodes->getFirstNode();
1410 myLastNode = mesh_p->Nodes->getLastNode();
1411 globalNodeIndex = mesh_p->Nodes->borrowGlobalNodesIndex();
1412 } else {
1413 throw FinleyAdapterException("Unsupported function space type for ownSample()");
1414 }
1415
1416 k=globalNodeIndex[id];
1417 return ((myFirstNode <= k) && (k < myLastNode));
1418 #endif
1419 }
1420 return true;
1421 }
1422
1423
1424 //
1425 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1426 //
1427 escript::ASM_ptr MeshAdapter::newSystemMatrix(const int row_blocksize,
1428 const escript::FunctionSpace& row_functionspace,
1429 const int column_blocksize,
1430 const escript::FunctionSpace& column_functionspace,
1431 const int type) const
1432 {
1433 // is the domain right?
1434 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1435 if (row_domain!=*this)
1436 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1437 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1438 if (col_domain!=*this)
1439 throw FinleyAdapterException("Error - domain of column function space does not match the domain of matrix generator.");
1440
1441 int reduceRowOrder=0;
1442 int reduceColOrder=0;
1443 // is the function space type right?
1444 if (row_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1445 reduceRowOrder=1;
1446 } else if (row_functionspace.getTypeCode() != DegreesOfFreedom) {
1447 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1448 }
1449 if (column_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1450 reduceColOrder=1;
1451 } else if (column_functionspace.getTypeCode() != DegreesOfFreedom) {
1452 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1453 }
1454
1455 // generate matrix:
1456 Paso_SystemMatrixPattern* fsystemMatrixPattern=
1457 getFinley_Mesh()->getPattern(reduceRowOrder, reduceColOrder);
1458 checkFinleyError();
1459 Paso_SystemMatrix* fsystemMatrix;
1460 const int trilinos = 0;
1461 if (trilinos) {
1462 #ifdef TRILINOS
1463 // FIXME: Allocation Epetra_VrbMatrix here...
1464 #endif
1465 } else {
1466 fsystemMatrix=Paso_SystemMatrix_alloc(type, fsystemMatrixPattern,
1467 row_blocksize, column_blocksize, FALSE);
1468 }
1469 checkPasoError();
1470 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1471 SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1472 return escript::ASM_ptr(sma);
1473 }
1474
1475 //
1476 // creates a TransportProblemAdapter
1477 //
1478 escript::ATP_ptr MeshAdapter::newTransportProblem(const int blocksize,
1479 const escript::FunctionSpace& functionspace, const int type) const
1480 {
1481 // is the domain right?
1482 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1483 if (domain!=*this)
1484 throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1485
1486 // is the function space type right?
1487 int reduceOrder=0;
1488 if (functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1489 reduceOrder=1;
1490 } else if (functionspace.getTypeCode() != DegreesOfFreedom) {
1491 throw FinleyAdapterException("Error - illegal function space type for transport problem.");
1492 }
1493
1494 // generate transport problem:
1495 Paso_SystemMatrixPattern* fsystemMatrixPattern=
1496 getFinley_Mesh()->getPattern(reduceOrder, reduceOrder);
1497 checkFinleyError();
1498 Paso_TransportProblem* transportProblem;
1499 transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern, blocksize);
1500 checkPasoError();
1501 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
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 void MeshAdapter::addDiracPoints(const vector<double>& points,
2166 const vector<int>& tags) const
2167 {
2168 // points will be flattened
2169 const int dim = getDim();
2170 int numPoints=points.size()/dim;
2171 int numTags=tags.size();
2172 Mesh* mesh=m_finleyMesh.get();
2173
2174 if ( points.size() % dim != 0 ) {
2175 throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2176 }
2177
2178 if ((numTags > 0) && (numPoints != numTags))
2179 throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2180
2181 mesh->addPoints(numPoints, &points[0], &tags[0]);
2182 checkFinleyError();
2183 }
2184
2185 // void MeshAdapter::addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2186 // {
2187 // const int dim = getDim();
2188 // int numPoints=boost::python::extract<int>(points.attr("__len__")());
2189 // int numTags=boost::python::extract<int>(tags.attr("__len__")());
2190 // Mesh* mesh=m_finleyMesh.get();
2191 //
2192 // if ( (numTags > 0) && ( numPoints != numTags ) )
2193 // throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2194 //
2195 // double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2196 // int* tags_ptr= TMPMEMALLOC(numPoints, int);
2197 //
2198 // for (int i=0;i<numPoints;++i) {
2199 // int tag_id=-1;
2200 // int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2201 // if ( numComps != dim ) {
2202 // stringstream temp;
2203 // temp << "Error - illegal number of components " << numComps << " for point " << i;
2204 // throw FinleyAdapterException(temp.str());
2205 // }
2206 // points_ptr[ i * dim ] = boost::python::extract<double>(points[i][0]);
2207 // if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2208 // if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2209 //
2210 // if ( numTags > 0) {
2211 // boost::python::extract<string> ex_str(tags[i]);
2212 // if ( ex_str.check() ) {
2213 // tag_id=getTag( ex_str());
2214 // } else {
2215 // boost::python::extract<int> ex_int(tags[i]);
2216 // if ( ex_int.check() ) {
2217 // tag_id=ex_int();
2218 // } else {
2219 // stringstream temp;
2220 // temp << "Error - unable to extract tag for point " << i;
2221 // throw FinleyAdapterException(temp.str());
2222 // }
2223 // }
2224 // }
2225 // tags_ptr[i]=tag_id;
2226 // }
2227 //
2228 // Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2229 // checkPasoError();
2230 //
2231 // TMPMEMFREE(points_ptr);
2232 // TMPMEMFREE(tags_ptr);
2233 // }
2234
2235 /*
2236 void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2237 {
2238 boost::python::list points = boost::python::list();
2239 boost::python::list tags = boost::python::list();
2240 points.append(point);
2241 tags.append(tag);
2242 addDiracPoints(points, tags);
2243 }
2244 */
2245
2246 /*
2247 void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const string& tag) const
2248 {
2249 boost::python::list points = boost::python::list();
2250 boost::python::list tags = boost::python::list();
2251 points.append(point);
2252 tags.append(tag);
2253 addDiracPoints(points, tags);
2254 }
2255 */
2256 } // end of namespace
2257

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26