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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3344 - (show annotations)
Thu Nov 11 23:26:52 2010 UTC (8 years, 9 months ago) by caltinay
File size: 85874 byte(s)
Phew!
-escript, finley, and dudley now uses weipa's saveVTK implementation
-moved tests from finley to weipa accordingly; dudley still to do
-rebaselined all test files
-fixed a few issues in weipa.saveVTK, e.g. saving metadata without schema
-added a deprecation warning to esys.escript.util.saveVTK
-todo: change doco, tests and other places to use weipa.saveVTK


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26