/[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 3490 - (show annotations)
Wed Mar 30 02:24:33 2011 UTC (8 years, 5 months ago) by caltinay
File size: 85874 byte(s)
More gcc-4.6 fixes (mostly initialized-but-unused-var warnings)

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
27 using namespace std;
28 using namespace escript;
29
30 namespace finley {
31
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(Finley_Mesh* finleyMesh)
50 {
51 setFunctionSpaceTypeNames();
52 //
53 // need to use a null_deleter as Finley_Mesh_free deletes the pointer
54 // for us.
55 m_finleyMesh.reset(finleyMesh,null_deleter());
56 }
57
58 //
59 // The copy constructor should just increment the use count
60 MeshAdapter::MeshAdapter(const MeshAdapter& in):
61 m_finleyMesh(in.m_finleyMesh)
62 {
63 setFunctionSpaceTypeNames();
64 }
65
66 MeshAdapter::~MeshAdapter()
67 {
68 //
69 // I hope the case for the pointer being zero has been taken care of.
70 // cout << "In MeshAdapter destructor." << endl;
71 if (m_finleyMesh.unique()) {
72 Finley_Mesh_free(m_finleyMesh.get());
73 }
74 }
75
76 int MeshAdapter::getMPISize() const
77 {
78 return m_finleyMesh.get()->MPIInfo->size;
79 }
80 int MeshAdapter::getMPIRank() const
81 {
82 return m_finleyMesh.get()->MPIInfo->rank;
83 }
84 void MeshAdapter::MPIBarrier() const
85 {
86 #ifdef ESYS_MPI
87 MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
88 #endif
89 return;
90 }
91 bool MeshAdapter::onMasterProcessor() const
92 {
93 return m_finleyMesh.get()->MPIInfo->rank == 0;
94 }
95
96
97 #ifdef ESYS_MPI
98 MPI_Comm
99 #else
100 unsigned int
101 #endif
102 MeshAdapter::getMPIComm() const
103 {
104 #ifdef ESYS_MPI
105 return m_finleyMesh->MPIInfo->comm;
106 #else
107 return 0;
108 #endif
109 }
110
111
112 Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
113 return m_finleyMesh.get();
114 }
115
116 void MeshAdapter::write(const string& fileName) const
117 {
118 char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
119 strcpy(fName,fileName.c_str());
120 Finley_Mesh_write(m_finleyMesh.get(),fName);
121 checkFinleyError();
122 TMPMEMFREE(fName);
123 }
124
125 void MeshAdapter::Print_Mesh_Info(const bool full) const
126 {
127 Finley_PrintMesh_Info(m_finleyMesh.get(), full);
128 }
129
130 void MeshAdapter::dump(const string& fileName) const
131 {
132 #ifdef USE_NETCDF
133 const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
134 NcVar *ids;
135 int *int_ptr;
136 Finley_Mesh *mesh = m_finleyMesh.get();
137 Finley_TagMap* tag_map;
138 int num_Tags = 0;
139 int mpi_size = mesh->MPIInfo->size;
140 int mpi_rank = mesh->MPIInfo->rank;
141 int numDim = mesh->Nodes->numDim;
142 int numNodes = mesh->Nodes->numNodes;
143 int num_Elements = mesh->Elements->numElements;
144 int num_FaceElements = mesh->FaceElements->numElements;
145 int num_ContactElements = mesh->ContactElements->numElements;
146 int num_Points = mesh->Points->numElements;
147 int num_Elements_numNodes = mesh->Elements->numNodes;
148 int num_FaceElements_numNodes = mesh->FaceElements->numNodes;
149 int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
150 #ifdef ESYS_MPI
151 MPI_Status status;
152 #endif
153
154 /* Incoming token indicates it's my turn to write */
155 #ifdef ESYS_MPI
156 if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
157 #endif
158
159 char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
160 mpi_size, mpi_rank);
161
162 /* Figure out how much storage is required for tags */
163 tag_map = mesh->TagMap;
164 num_Tags = 0;
165 while (tag_map) {
166 num_Tags++;
167 tag_map=tag_map->next;
168 }
169
170 // NetCDF error handler
171 NcError err(NcError::verbose_nonfatal);
172 // Create the file.
173 NcFile dataFile(newFileName, NcFile::Replace);
174 string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
175 // check if writing was successful
176 if (!dataFile.is_valid())
177 throw DataException(msgPrefix+"Open file for output");
178
179 // Define dimensions (num_Elements and dim_Elements are identical,
180 // dim_Elements only appears if > 0)
181 if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
182 throw DataException(msgPrefix+"add_dim(numNodes)");
183 if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
184 throw DataException(msgPrefix+"add_dim(numDim)");
185 if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )
186 throw DataException(msgPrefix+"add_dim(mpi_size)");
187 if (num_Elements>0)
188 if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
189 throw DataException(msgPrefix+"add_dim(dim_Elements)");
190 if (num_FaceElements>0)
191 if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
192 throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
193 if (num_ContactElements>0)
194 if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
195 throw DataException(msgPrefix+"add_dim(dim_ContactElements)");
196 if (num_Points>0)
197 if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
198 throw DataException(msgPrefix+"add_dim(dim_Points)");
199 if (num_Elements>0)
200 if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
201 throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
202 if (num_FaceElements>0)
203 if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
204 throw DataException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
205 if (num_ContactElements>0)
206 if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
207 throw DataException(msgPrefix+"add_dim(dim_ContactElements_numNodes)");
208 if (num_Tags>0)
209 if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
210 throw DataException(msgPrefix+"add_dim(dim_Tags)");
211
212 // Attributes: MPI size, MPI rank, Name, order, reduced_order
213 if (!dataFile.add_att("mpi_size", mpi_size) )
214 throw DataException(msgPrefix+"add_att(mpi_size)");
215 if (!dataFile.add_att("mpi_rank", mpi_rank) )
216 throw DataException(msgPrefix+"add_att(mpi_rank)");
217 if (!dataFile.add_att("Name",mesh->Name) )
218 throw DataException(msgPrefix+"add_att(Name)");
219 if (!dataFile.add_att("numDim",numDim) )
220 throw DataException(msgPrefix+"add_att(order)");
221 if (!dataFile.add_att("order",mesh->integrationOrder) )
222 throw DataException(msgPrefix+"add_att(order)");
223 if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
224 throw DataException(msgPrefix+"add_att(reduced_order)");
225 if (!dataFile.add_att("numNodes",numNodes) )
226 throw DataException(msgPrefix+"add_att(numNodes)");
227 if (!dataFile.add_att("num_Elements",num_Elements) )
228 throw DataException(msgPrefix+"add_att(num_Elements)");
229 if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
230 throw DataException(msgPrefix+"add_att(num_FaceElements)");
231 if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
232 throw DataException(msgPrefix+"add_att(num_ContactElements)");
233 if (!dataFile.add_att("num_Points",num_Points) )
234 throw DataException(msgPrefix+"add_att(num_Points)");
235 if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
236 throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
237 if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
238 throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
239 if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
240 throw DataException(msgPrefix+"add_att(num_ContactElements_numNodes)");
241 if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
242 throw DataException(msgPrefix+"add_att(Elements_TypeId)");
243 if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
244 throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
245 if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
246 throw DataException(msgPrefix+"add_att(ContactElements_TypeId)");
247 if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
248 throw DataException(msgPrefix+"add_att(Points_TypeId)");
249 if (!dataFile.add_att("num_Tags", num_Tags) )
250 throw DataException(msgPrefix+"add_att(num_Tags)");
251
252 // // // // // Nodes // // // // //
253
254 // Nodes nodeDistribution
255 if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
256 throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
257 int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
258 if (! (ids->put(int_ptr, mpi_size+1)) )
259 throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
260
261 // Nodes degreesOfFreedomDistribution
262 if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
263 throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
264 int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
265 if (! (ids->put(int_ptr, mpi_size+1)) )
266 throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
267
268 // Only write nodes if non-empty because NetCDF doesn't like empty arrays
269 // (it treats them as NC_UNLIMITED)
270 if (numNodes>0) {
271
272 // Nodes Id
273 if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
274 throw DataException(msgPrefix+"add_var(Nodes_Id)");
275 int_ptr = &mesh->Nodes->Id[0];
276 if (! (ids->put(int_ptr, numNodes)) )
277 throw DataException(msgPrefix+"put(Nodes_Id)");
278
279 // Nodes Tag
280 if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
281 throw DataException(msgPrefix+"add_var(Nodes_Tag)");
282 int_ptr = &mesh->Nodes->Tag[0];
283 if (! (ids->put(int_ptr, numNodes)) )
284 throw DataException(msgPrefix+"put(Nodes_Tag)");
285
286 // Nodes gDOF
287 if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
288 throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
289 int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
290 if (! (ids->put(int_ptr, numNodes)) )
291 throw DataException(msgPrefix+"put(Nodes_gDOF)");
292
293 // Nodes global node index
294 if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
295 throw DataException(msgPrefix+"add_var(Nodes_gNI)");
296 int_ptr = &mesh->Nodes->globalNodesIndex[0];
297 if (! (ids->put(int_ptr, numNodes)) )
298 throw DataException(msgPrefix+"put(Nodes_gNI)");
299
300 // Nodes grDof
301 if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
302 throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
303 int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
304 if (! (ids->put(int_ptr, numNodes)) )
305 throw DataException(msgPrefix+"put(Nodes_grDfI)");
306
307 // Nodes grNI
308 if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
309 throw DataException(msgPrefix+"add_var(Nodes_grNI)");
310 int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
311 if (! (ids->put(int_ptr, numNodes)) )
312 throw DataException(msgPrefix+"put(Nodes_grNI)");
313
314 // Nodes Coordinates
315 if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
316 throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
317 if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
318 throw DataException(msgPrefix+"put(Nodes_Coordinates)");
319
320 }
321
322 // // // // // Elements // // // // //
323
324 if (num_Elements>0) {
325
326 // Elements_Id
327 if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
328 throw DataException(msgPrefix+"add_var(Elements_Id)");
329 int_ptr = &mesh->Elements->Id[0];
330 if (! (ids->put(int_ptr, num_Elements)) )
331 throw DataException(msgPrefix+"put(Elements_Id)");
332
333 // Elements_Tag
334 if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
335 throw DataException(msgPrefix+"add_var(Elements_Tag)");
336 int_ptr = &mesh->Elements->Tag[0];
337 if (! (ids->put(int_ptr, num_Elements)) )
338 throw DataException(msgPrefix+"put(Elements_Tag)");
339
340 // Elements_Owner
341 if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
342 throw DataException(msgPrefix+"add_var(Elements_Owner)");
343 int_ptr = &mesh->Elements->Owner[0];
344 if (! (ids->put(int_ptr, num_Elements)) )
345 throw DataException(msgPrefix+"put(Elements_Owner)");
346
347 // Elements_Color
348 if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
349 throw DataException(msgPrefix+"add_var(Elements_Color)");
350 int_ptr = &mesh->Elements->Color[0];
351 if (! (ids->put(int_ptr, num_Elements)) )
352 throw DataException(msgPrefix+"put(Elements_Color)");
353
354 // Elements_Nodes
355 if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
356 throw DataException(msgPrefix+"add_var(Elements_Nodes)");
357 if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
358 throw DataException(msgPrefix+"put(Elements_Nodes)");
359
360 }
361
362 // // // // // Face_Elements // // // // //
363
364 if (num_FaceElements>0) {
365
366 // FaceElements_Id
367 if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
368 throw DataException(msgPrefix+"add_var(FaceElements_Id)");
369 int_ptr = &mesh->FaceElements->Id[0];
370 if (! (ids->put(int_ptr, num_FaceElements)) )
371 throw DataException(msgPrefix+"put(FaceElements_Id)");
372
373 // FaceElements_Tag
374 if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
375 throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
376 int_ptr = &mesh->FaceElements->Tag[0];
377 if (! (ids->put(int_ptr, num_FaceElements)) )
378 throw DataException(msgPrefix+"put(FaceElements_Tag)");
379
380 // FaceElements_Owner
381 if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
382 throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
383 int_ptr = &mesh->FaceElements->Owner[0];
384 if (! (ids->put(int_ptr, num_FaceElements)) )
385 throw DataException(msgPrefix+"put(FaceElements_Owner)");
386
387 // FaceElements_Color
388 if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
389 throw DataException(msgPrefix+"add_var(FaceElements_Color)");
390 int_ptr = &mesh->FaceElements->Color[0];
391 if (! (ids->put(int_ptr, num_FaceElements)) )
392 throw DataException(msgPrefix+"put(FaceElements_Color)");
393
394 // FaceElements_Nodes
395 if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
396 throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
397 if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
398 throw DataException(msgPrefix+"put(FaceElements_Nodes)");
399
400 }
401
402 // // // // // Contact_Elements // // // // //
403
404 if (num_ContactElements>0) {
405
406 // ContactElements_Id
407 if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
408 throw DataException(msgPrefix+"add_var(ContactElements_Id)");
409 int_ptr = &mesh->ContactElements->Id[0];
410 if (! (ids->put(int_ptr, num_ContactElements)) )
411 throw DataException(msgPrefix+"put(ContactElements_Id)");
412
413 // ContactElements_Tag
414 if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
415 throw DataException(msgPrefix+"add_var(ContactElements_Tag)");
416 int_ptr = &mesh->ContactElements->Tag[0];
417 if (! (ids->put(int_ptr, num_ContactElements)) )
418 throw DataException(msgPrefix+"put(ContactElements_Tag)");
419
420 // ContactElements_Owner
421 if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
422 throw DataException(msgPrefix+"add_var(ContactElements_Owner)");
423 int_ptr = &mesh->ContactElements->Owner[0];
424 if (! (ids->put(int_ptr, num_ContactElements)) )
425 throw DataException(msgPrefix+"put(ContactElements_Owner)");
426
427 // ContactElements_Color
428 if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
429 throw DataException(msgPrefix+"add_var(ContactElements_Color)");
430 int_ptr = &mesh->ContactElements->Color[0];
431 if (! (ids->put(int_ptr, num_ContactElements)) )
432 throw DataException(msgPrefix+"put(ContactElements_Color)");
433
434 // ContactElements_Nodes
435 if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
436 throw DataException(msgPrefix+"add_var(ContactElements_Nodes)");
437 if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
438 throw DataException(msgPrefix+"put(ContactElements_Nodes)");
439
440 }
441
442 // // // // // Points // // // // //
443
444 if (num_Points>0) {
445
446 fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");
447
448 // Points_Id
449 if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
450 throw DataException(msgPrefix+"add_var(Points_Id)");
451 int_ptr = &mesh->Points->Id[0];
452 if (! (ids->put(int_ptr, num_Points)) )
453 throw DataException(msgPrefix+"put(Points_Id)");
454
455 // Points_Tag
456 if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
457 throw DataException(msgPrefix+"add_var(Points_Tag)");
458 int_ptr = &mesh->Points->Tag[0];
459 if (! (ids->put(int_ptr, num_Points)) )
460 throw DataException(msgPrefix+"put(Points_Tag)");
461
462 // Points_Owner
463 if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
464 throw DataException(msgPrefix+"add_var(Points_Owner)");
465 int_ptr = &mesh->Points->Owner[0];
466 if (! (ids->put(int_ptr, num_Points)) )
467 throw DataException(msgPrefix+"put(Points_Owner)");
468
469 // Points_Color
470 if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
471 throw DataException(msgPrefix+"add_var(Points_Color)");
472 int_ptr = &mesh->Points->Color[0];
473 if (! (ids->put(int_ptr, num_Points)) )
474 throw DataException(msgPrefix+"put(Points_Color)");
475
476 // Points_Nodes
477 // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
478 if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
479 throw DataException(msgPrefix+"add_var(Points_Nodes)");
480 if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
481 throw DataException(msgPrefix+"put(Points_Nodes)");
482
483 }
484
485 // // // // // TagMap // // // // //
486
487 if (num_Tags>0) {
488
489 // Temp storage to gather node IDs
490 int *Tags_keys = TMPMEMALLOC(num_Tags, int);
491 char name_temp[4096];
492
493 /* Copy tag data into temp arrays */
494 tag_map = mesh->TagMap;
495 if (tag_map) {
496 int i = 0;
497 while (tag_map) {
498 Tags_keys[i++] = tag_map->tag_key;
499 tag_map=tag_map->next;
500 }
501 }
502
503 // Tags_keys
504 if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
505 throw DataException(msgPrefix+"add_var(Tags_keys)");
506 int_ptr = &Tags_keys[0];
507 if (! (ids->put(int_ptr, num_Tags)) )
508 throw DataException(msgPrefix+"put(Tags_keys)");
509
510 // Tags_names_*
511 // This is an array of strings, it should be stored as an array but
512 // instead I have hacked in one attribute per string because the NetCDF
513 // manual doesn't tell how to do an array of strings
514 tag_map = mesh->TagMap;
515 if (tag_map) {
516 int i = 0;
517 while (tag_map) {
518 sprintf(name_temp, "Tags_name_%d", i);
519 if (!dataFile.add_att(name_temp, tag_map->name) )
520 throw DataException(msgPrefix+"add_att(Tags_names_XX)");
521 tag_map=tag_map->next;
522 i++;
523 }
524 }
525
526 TMPMEMFREE(Tags_keys);
527 }
528
529 /* Send token to next MPI process so he can take his turn */
530 #ifdef ESYS_MPI
531 if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
532 #endif
533
534 // NetCDF file is closed by destructor of NcFile object
535
536 #else
537 Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
538 #endif /* USE_NETCDF */
539 checkFinleyError();
540 }
541
542 string MeshAdapter::getDescription() const
543 {
544 return "FinleyMesh";
545 }
546
547 string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
548 {
549 FunctionSpaceNamesMapType::iterator loc;
550 loc=m_functionSpaceTypeNames.find(functionSpaceType);
551 if (loc==m_functionSpaceTypeNames.end()) {
552 return "Invalid function space type code.";
553 } else {
554 return loc->second;
555 }
556 }
557
558 bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
559 {
560 FunctionSpaceNamesMapType::iterator loc;
561 loc=m_functionSpaceTypeNames.find(functionSpaceType);
562 return (loc!=m_functionSpaceTypeNames.end());
563 }
564
565 void MeshAdapter::setFunctionSpaceTypeNames()
566 {
567 m_functionSpaceTypeNames.insert
568 (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));
569 m_functionSpaceTypeNames.insert
570 (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));
571 m_functionSpaceTypeNames.insert
572 (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
573 m_functionSpaceTypeNames.insert
574 (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));
575 m_functionSpaceTypeNames.insert
576 (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
577 m_functionSpaceTypeNames.insert
578 (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
579 m_functionSpaceTypeNames.insert
580 (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
581 m_functionSpaceTypeNames.insert
582 (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
583 m_functionSpaceTypeNames.insert
584 (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
585 m_functionSpaceTypeNames.insert
586 (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
587 m_functionSpaceTypeNames.insert
588 (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
589 m_functionSpaceTypeNames.insert
590 (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
591 m_functionSpaceTypeNames.insert
592 (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
593 }
594
595 int MeshAdapter::getContinuousFunctionCode() const
596 {
597 return Nodes;
598 }
599 int MeshAdapter::getReducedContinuousFunctionCode() const
600 {
601 return ReducedNodes;
602 }
603
604 int MeshAdapter::getFunctionCode() const
605 {
606 return Elements;
607 }
608 int MeshAdapter::getReducedFunctionCode() const
609 {
610 return ReducedElements;
611 }
612
613 int MeshAdapter::getFunctionOnBoundaryCode() const
614 {
615 return FaceElements;
616 }
617 int MeshAdapter::getReducedFunctionOnBoundaryCode() const
618 {
619 return ReducedFaceElements;
620 }
621
622 int MeshAdapter::getFunctionOnContactZeroCode() const
623 {
624 return ContactElementsZero;
625 }
626 int MeshAdapter::getReducedFunctionOnContactZeroCode() const
627 {
628 return ReducedContactElementsZero;
629 }
630
631 int MeshAdapter::getFunctionOnContactOneCode() const
632 {
633 return ContactElementsOne;
634 }
635 int MeshAdapter::getReducedFunctionOnContactOneCode() const
636 {
637 return ReducedContactElementsOne;
638 }
639
640 int MeshAdapter::getSolutionCode() const
641 {
642 return DegreesOfFreedom;
643 }
644
645 int MeshAdapter::getReducedSolutionCode() const
646 {
647 return ReducedDegreesOfFreedom;
648 }
649
650 int MeshAdapter::getDiracDeltaFunctionCode() const
651 {
652 return Points;
653 }
654
655 //
656 // return the spatial dimension of the Mesh:
657 //
658 int MeshAdapter::getDim() const
659 {
660 Finley_Mesh* mesh=m_finleyMesh.get();
661 int numDim=Finley_Mesh_getDim(mesh);
662 checkFinleyError();
663 return numDim;
664 }
665
666 //
667 // Return the number of data points summed across all MPI processes
668 //
669 int MeshAdapter::getNumDataPointsGlobal() const
670 {
671 return Finley_NodeFile_getGlobalNumNodes(m_finleyMesh.get()->Nodes);
672 }
673
674 //
675 // return the number of data points per sample and the number of samples
676 // needed to represent data on a parts of the mesh.
677 //
678 pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
679 {
680 int numDataPointsPerSample=0;
681 int numSamples=0;
682 Finley_Mesh* mesh=m_finleyMesh.get();
683 switch (functionSpaceCode) {
684 case(Nodes):
685 numDataPointsPerSample=1;
686 numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);
687 break;
688 case(ReducedNodes):
689 numDataPointsPerSample=1;
690 numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
691 break;
692 case(Elements):
693 if (mesh->Elements!=NULL) {
694 numSamples=mesh->Elements->numElements;
695 numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
696 }
697 break;
698 case(ReducedElements):
699 if (mesh->Elements!=NULL) {
700 numSamples=mesh->Elements->numElements;
701 numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
702 }
703 break;
704 case(FaceElements):
705 if (mesh->FaceElements!=NULL) {
706 numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
707 numSamples=mesh->FaceElements->numElements;
708 }
709 break;
710 case(ReducedFaceElements):
711 if (mesh->FaceElements!=NULL) {
712 numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
713 numSamples=mesh->FaceElements->numElements;
714 }
715 break;
716 case(Points):
717 if (mesh->Points!=NULL) {
718 numDataPointsPerSample=1;
719 numSamples=mesh->Points->numElements;
720 }
721 break;
722 case(ContactElementsZero):
723 if (mesh->ContactElements!=NULL) {
724 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
725 numSamples=mesh->ContactElements->numElements;
726 }
727 break;
728 case(ReducedContactElementsZero):
729 if (mesh->ContactElements!=NULL) {
730 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
731 numSamples=mesh->ContactElements->numElements;
732 }
733 break;
734 case(ContactElementsOne):
735 if (mesh->ContactElements!=NULL) {
736 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
737 numSamples=mesh->ContactElements->numElements;
738 }
739 break;
740 case(ReducedContactElementsOne):
741 if (mesh->ContactElements!=NULL) {
742 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
743 numSamples=mesh->ContactElements->numElements;
744 }
745 break;
746 case(DegreesOfFreedom):
747 if (mesh->Nodes!=NULL) {
748 numDataPointsPerSample=1;
749 numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
750 }
751 break;
752 case(ReducedDegreesOfFreedom):
753 if (mesh->Nodes!=NULL) {
754 numDataPointsPerSample=1;
755 numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
756 }
757 break;
758 default:
759 stringstream temp;
760 temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
761 throw FinleyAdapterException(temp.str());
762 break;
763 }
764 return pair<int,int>(numDataPointsPerSample,numSamples);
765 }
766
767 //
768 // adds linear PDE of second order into a given stiffness matrix and right hand side:
769 //
770 void MeshAdapter::addPDEToSystem(
771 AbstractSystemMatrix& mat, escript::Data& rhs,
772 const escript::Data& A, const escript::Data& B, const escript::Data& C,const escript::Data& D,const escript::Data& X,const escript::Data& Y,
773 const escript::Data& d, const escript::Data& y,
774 const escript::Data& d_contact,const escript::Data& y_contact) const
775 {
776 SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
777 if (smat==0)
778 {
779 throw FinleyAdapterException("finley only supports its own system matrices.");
780 }
781 escriptDataC _rhs=rhs.getDataC();
782 escriptDataC _A =A.getDataC();
783 escriptDataC _B=B.getDataC();
784 escriptDataC _C=C.getDataC();
785 escriptDataC _D=D.getDataC();
786 escriptDataC _X=X.getDataC();
787 escriptDataC _Y=Y.getDataC();
788 escriptDataC _d=d.getDataC();
789 escriptDataC _y=y.getDataC();
790 escriptDataC _d_contact=d_contact.getDataC();
791 escriptDataC _y_contact=y_contact.getDataC();
792
793 Finley_Mesh* mesh=m_finleyMesh.get();
794
795 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
796 checkFinleyError();
797
798 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
799 checkFinleyError();
800
801 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
802 checkFinleyError();
803 }
804
805 void MeshAdapter::addPDEToLumpedSystem(
806 escript::Data& mat,
807 const escript::Data& D,
808 const escript::Data& d,
809 const bool useHRZ) 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, useHRZ);
818 checkFinleyError();
819
820 Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
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(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1554 metadata, metadata_schema, arg);
1555 }
1556
1557 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1558 {
1559 #ifdef ESYS_MPI
1560 index_t myFirstNode=0, myLastNode=0, k=0;
1561 index_t* globalNodeIndex=0;
1562 Finley_Mesh* mesh_p=m_finleyMesh.get();
1563 if (fs_code == FINLEY_REDUCED_NODES)
1564 {
1565 myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1566 myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1567 globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1568 }
1569 else
1570 {
1571 myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1572 myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1573 globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1574 }
1575 k=globalNodeIndex[id];
1576 return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1577 #endif
1578 return true;
1579 }
1580
1581
1582
1583 //
1584 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1585 //
1586 ASM_ptr MeshAdapter::newSystemMatrix(
1587 const int row_blocksize,
1588 const escript::FunctionSpace& row_functionspace,
1589 const int column_blocksize,
1590 const escript::FunctionSpace& column_functionspace,
1591 const int type) const
1592 {
1593 int reduceRowOrder=0;
1594 int reduceColOrder=0;
1595 // is the domain right?
1596 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1597 if (row_domain!=*this)
1598 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1599 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1600 if (col_domain!=*this)
1601 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1602 // is the function space type right
1603 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1604 reduceRowOrder=0;
1605 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1606 reduceRowOrder=1;
1607 } else {
1608 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1609 }
1610 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1611 reduceColOrder=0;
1612 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1613 reduceColOrder=1;
1614 } else {
1615 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1616 }
1617 // generate matrix:
1618
1619 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1620 checkFinleyError();
1621 Paso_SystemMatrix* fsystemMatrix;
1622 int trilinos = 0;
1623 if (trilinos) {
1624 #ifdef TRILINOS
1625 /* Allocation Epetra_VrbMatrix here */
1626 #endif
1627 }
1628 else {
1629 fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1630 }
1631 checkPasoError();
1632 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1633 SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1634 return ASM_ptr(sma);
1635 // return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1636 }
1637
1638 //
1639 // creates a TransportProblemAdapter
1640 //
1641 ATP_ptr MeshAdapter::newTransportProblem(
1642 const bool useBackwardEuler,
1643 const int blocksize,
1644 const escript::FunctionSpace& functionspace,
1645 const int type) const
1646 {
1647 int reduceOrder=0;
1648 // is the domain right?
1649 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1650 if (domain!=*this)
1651 throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1652 // is the function space type right
1653 if (functionspace.getTypeCode()==DegreesOfFreedom) {
1654 reduceOrder=0;
1655 } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1656 reduceOrder=1;
1657 } else {
1658 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1659 }
1660 // generate matrix:
1661
1662 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1663 checkFinleyError();
1664 Paso_TransportProblem* transportProblem;
1665 transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1666 checkPasoError();
1667 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1668 TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1669 return ATP_ptr(tpa);
1670 // return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1671 }
1672
1673 //
1674 // vtkObject MeshAdapter::createVtkObject() const
1675 // TODO:
1676 //
1677
1678 //
1679 // returns true if data at the atom_type is considered as being cell centered:
1680 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1681 {
1682 switch(functionSpaceCode) {
1683 case(Nodes):
1684 case(DegreesOfFreedom):
1685 case(ReducedDegreesOfFreedom):
1686 return false;
1687 break;
1688 case(Elements):
1689 case(FaceElements):
1690 case(Points):
1691 case(ContactElementsZero):
1692 case(ContactElementsOne):
1693 case(ReducedElements):
1694 case(ReducedFaceElements):
1695 case(ReducedContactElementsZero):
1696 case(ReducedContactElementsOne):
1697 return true;
1698 break;
1699 default:
1700 stringstream temp;
1701 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1702 throw FinleyAdapterException(temp.str());
1703 break;
1704 }
1705 checkFinleyError();
1706 return false;
1707 }
1708
1709 bool
1710 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1711 {
1712 /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1713 class 1: DOF <-> Nodes
1714 class 2: ReducedDOF <-> ReducedNodes
1715 class 3: Points
1716 class 4: Elements
1717 class 5: ReducedElements
1718 class 6: FaceElements
1719 class 7: ReducedFaceElements
1720 class 8: ContactElementZero <-> ContactElementOne
1721 class 9: ReducedContactElementZero <-> ReducedContactElementOne
1722
1723 There is also a set of lines. Interpolation is possible down a line but not between lines.
1724 class 1 and 2 belong to all lines so aren't considered.
1725 line 0: class 3
1726 line 1: class 4,5
1727 line 2: class 6,7
1728 line 3: class 8,9
1729
1730 For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1731 eg hasnodes is true if we have at least one instance of Nodes.
1732 */
1733 if (fs.empty())
1734 {
1735 return false;
1736 }
1737 vector<int> hasclass(10);
1738 vector<int> hasline(4);
1739 bool hasnodes=false;
1740 bool hasrednodes=false;
1741 bool hascez=false;
1742 bool hasrcez=false;
1743 for (int i=0;i<fs.size();++i)
1744 {
1745 switch(fs[i])
1746 {
1747 case(Nodes): hasnodes=true; // no break is deliberate
1748 case(DegreesOfFreedom):
1749 hasclass[1]=1;
1750 break;
1751 case(ReducedNodes): hasrednodes=true; // no break is deliberate
1752 case(ReducedDegreesOfFreedom):
1753 hasclass[2]=1;
1754 break;
1755 case(Points):
1756 hasline[0]=1;
1757 hasclass[3]=1;
1758 break;
1759 case(Elements):
1760 hasclass[4]=1;
1761 hasline[1]=1;
1762 break;
1763 case(ReducedElements):
1764 hasclass[5]=1;
1765 hasline[1]=1;
1766 break;
1767 case(FaceElements):
1768 hasclass[6]=1;
1769 hasline[2]=1;
1770 break;
1771 case(ReducedFaceElements):
1772 hasclass[7]=1;
1773 hasline[2]=1;
1774 break;
1775 case(ContactElementsZero): hascez=true; // no break is deliberate
1776 case(ContactElementsOne):
1777 hasclass[8]=1;
1778 hasline[3]=1;
1779 break;
1780 case(ReducedContactElementsZero): hasrcez=true; // no break is deliberate
1781 case(ReducedContactElementsOne):
1782 hasclass[9]=1;
1783 hasline[3]=1;
1784 break;
1785 default:
1786 return false;
1787 }
1788 }
1789 int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1790 // fail if we have more than one leaf group
1791
1792 if (totlines>1)
1793 {
1794 return false; // there are at least two branches we can't interpolate between
1795 }
1796 else if (totlines==1)
1797 {
1798 if (hasline[0]==1) // we have points
1799 {
1800 resultcode=Points;
1801 }
1802 else if (hasline[1]==1)
1803 {
1804 if (hasclass[5]==1)
1805 {
1806 resultcode=ReducedElements;
1807 }
1808 else
1809 {
1810 resultcode=Elements;
1811 }
1812 }
1813 else if (hasline[2]==1)
1814 {
1815 if (hasclass[7]==1)
1816 {
1817 resultcode=ReducedFaceElements;
1818 }
1819 else
1820 {
1821 resultcode=FaceElements;
1822 }
1823 }
1824 else // so we must be in line3
1825 {
1826 if (hasclass[9]==1)
1827 {
1828 // need something from class 9
1829 resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1830 }
1831 else
1832 {
1833 // something from class 8
1834 resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1835 }
1836 }
1837 }
1838 else // totlines==0
1839 {
1840 if (hasclass[2]==1)
1841 {
1842 // something from class 2
1843 resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1844 }
1845 else
1846 { // something from class 1
1847 resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1848 }
1849 }
1850 return true;
1851 }
1852
1853 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1854 {
1855 switch(functionSpaceType_source) {
1856 case(Nodes):
1857 switch(functionSpaceType_target) {
1858 case(Nodes):
1859 case(ReducedNodes):
1860 case(ReducedDegreesOfFreedom):
1861 case(DegreesOfFreedom):
1862 case(Elements):
1863 case(ReducedElements):
1864 case(FaceElements):
1865 case(ReducedFaceElements):
1866 case(Points):
1867 case(ContactElementsZero):
1868 case(ReducedContactElementsZero):
1869 case(ContactElementsOne):
1870 case(ReducedContactElementsOne):
1871 return true;
1872 default:
1873 stringstream temp;
1874 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1875 throw FinleyAdapterException(temp.str());
1876 }
1877 break;
1878 case(ReducedNodes):
1879 switch(functionSpaceType_target) {
1880 case(ReducedNodes):
1881 case(ReducedDegreesOfFreedom):
1882 case(Elements):
1883 case(ReducedElements):
1884 case(FaceElements):
1885 case(ReducedFaceElements):
1886 case(Points):
1887 case(ContactElementsZero):
1888 case(ReducedContactElementsZero):
1889 case(ContactElementsOne):
1890 case(ReducedContactElementsOne):
1891 return true;
1892 case(Nodes):
1893 case(DegreesOfFreedom):
1894 return false;
1895 default:
1896 stringstream temp;
1897 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1898 throw FinleyAdapterException(temp.str());
1899 }
1900 break;
1901 case(Elements):
1902 if (functionSpaceType_target==Elements) {
1903 return true;
1904 } else if (functionSpaceType_target==ReducedElements) {
1905 return true;
1906 } else {
1907 return false;
1908 }
1909 case(ReducedElements):
1910 if (functionSpaceType_target==ReducedElements) {
1911 return true;
1912 } else {
1913 return false;
1914 }
1915 case(FaceElements):
1916 if (functionSpaceType_target==FaceElements) {
1917 return true;
1918 } else if (functionSpaceType_target==ReducedFaceElements) {
1919 return true;
1920 } else {
1921 return false;
1922 }
1923 case(ReducedFaceElements):
1924 if (functionSpaceType_target==ReducedFaceElements) {
1925 return true;
1926 } else {
1927 return false;
1928 }
1929 case(Points):
1930 if (functionSpaceType_target==Points) {
1931 return true;
1932 } else {
1933 return false;
1934 }
1935 case(ContactElementsZero):
1936 case(ContactElementsOne):
1937 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1938 return true;
1939 } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1940 return true;
1941 } else {
1942 return false;
1943 }
1944 case(ReducedContactElementsZero):
1945 case(ReducedContactElementsOne):
1946 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1947 return true;
1948 } else {
1949 return false;
1950 }
1951 case(DegreesOfFreedom):
1952 switch(functionSpaceType_target) {
1953 case(ReducedDegreesOfFreedom):
1954 case(DegreesOfFreedom):
1955 case(Nodes):
1956 case(ReducedNodes):
1957 case(Elements):
1958 case(ReducedElements):
1959 case(Points):
1960 case(FaceElements):
1961 case(ReducedFaceElements):
1962 case(ContactElementsZero):
1963 case(ReducedContactElementsZero):
1964 case(ContactElementsOne):
1965 case(ReducedContactElementsOne):
1966 return true;
1967 default:
1968 stringstream temp;
1969 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1970 throw FinleyAdapterException(temp.str());
1971 }
1972 break;
1973 case(ReducedDegreesOfFreedom):
1974 switch(functionSpaceType_target) {
1975 case(ReducedDegreesOfFreedom):
1976 case(ReducedNodes):
1977 case(Elements):
1978 case(ReducedElements):
1979 case(FaceElements):
1980 case(ReducedFaceElements):
1981 case(Points):
1982 case(ContactElementsZero):
1983 case(ReducedContactElementsZero):
1984 case(ContactElementsOne):
1985 case(ReducedContactElementsOne):
1986 return true;
1987 case(Nodes):
1988 case(DegreesOfFreedom):
1989 return false;
1990 default:
1991 stringstream temp;
1992 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1993 throw FinleyAdapterException(temp.str());
1994 }
1995 break;
1996 default:
1997 stringstream temp;
1998 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1999 throw FinleyAdapterException(temp.str());
2000 break;
2001 }
2002 checkFinleyError();
2003 return false;
2004 }
2005
2006 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
2007 {
2008 return false;
2009 }
2010
2011 bool MeshAdapter::operator==(const AbstractDomain& other) const
2012 {
2013 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2014 if (temp!=0) {
2015 return (m_finleyMesh==temp->m_finleyMesh);
2016 } else {
2017 return false;
2018 }
2019 }
2020
2021 bool MeshAdapter::operator!=(const AbstractDomain& other) const
2022 {
2023 return !(operator==(other));
2024 }
2025
2026 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2027 {
2028 Finley_Mesh* mesh=m_finleyMesh.get();
2029 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2030 checkPasoError();
2031 return out;
2032 }
2033 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2034 {
2035 Finley_Mesh* mesh=m_finleyMesh.get();
2036 int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2037 checkPasoError();
2038 return out;
2039 }
2040
2041 escript::Data MeshAdapter::getX() const
2042 {
2043 return continuousFunction(*this).getX();
2044 }
2045
2046 escript::Data MeshAdapter::getNormal() const
2047 {
2048 return functionOnBoundary(*this).getNormal();
2049 }
2050
2051 escript::Data MeshAdapter::getSize() const
2052 {
2053 return escript::function(*this).getSize();
2054 }
2055
2056 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2057 {
2058 int *out = NULL;
2059 Finley_Mesh* mesh=m_finleyMesh.get();
2060 switch (functionSpaceType) {
2061 case(Nodes):
2062 out=mesh->Nodes->Id;
2063 break;
2064 case(ReducedNodes):
2065 out=mesh->Nodes->reducedNodesId;
2066 break;
2067 case(Elements):
2068 out=mesh->Elements->Id;
2069 break;
2070 case(ReducedElements):
2071 out=mesh->Elements->Id;
2072 break;
2073 case(FaceElements):
2074 out=mesh->FaceElements->Id;
2075 break;
2076 case(ReducedFaceElements):
2077 out=mesh->FaceElements->Id;
2078 break;
2079 case(Points):
2080 out=mesh->Points->Id;
2081 break;
2082 case(ContactElementsZero):
2083 out=mesh->ContactElements->Id;
2084 break;
2085 case(ReducedContactElementsZero):
2086 out=mesh->ContactElements->Id;
2087 break;
2088 case(ContactElementsOne):
2089 out=mesh->ContactElements->Id;
2090 break;
2091 case(ReducedContactElementsOne):
2092 out=mesh->ContactElements->Id;
2093 break;
2094 case(DegreesOfFreedom):
2095 out=mesh->Nodes->degreesOfFreedomId;
2096 break;
2097 case(ReducedDegreesOfFreedom):
2098 out=mesh->Nodes->reducedDegreesOfFreedomId;
2099 break;
2100 default:
2101 stringstream temp;
2102 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2103 throw FinleyAdapterException(temp.str());
2104 break;
2105 }
2106 return out;
2107 }
2108 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2109 {
2110 int out=0;
2111 Finley_Mesh* mesh=m_finleyMesh.get();
2112 switch (functionSpaceType) {
2113 case(Nodes):
2114 out=mesh->Nodes->Tag[sampleNo];
2115 break;
2116 case(ReducedNodes):
2117 throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2118 break;
2119 case(Elements):
2120 out=mesh->Elements->Tag[sampleNo];
2121 break;
2122 case(ReducedElements):
2123 out=mesh->Elements->Tag[sampleNo];
2124 break;
2125 case(FaceElements):
2126 out=mesh->FaceElements->Tag[sampleNo];
2127 break;
2128 case(ReducedFaceElements):
2129 out=mesh->FaceElements->Tag[sampleNo];
2130 break;
2131 case(Points):
2132 out=mesh->Points->Tag[sampleNo];
2133 break;
2134 case(ContactElementsZero):
2135 out=mesh->ContactElements->Tag[sampleNo];
2136 break;
2137 case(ReducedContactElementsZero):
2138 out=mesh->ContactElements->Tag[sampleNo];
2139 break;
2140 case(ContactElementsOne):
2141 out=mesh->ContactElements->Tag[sampleNo];
2142 break;
2143 case(ReducedContactElementsOne):
2144 out=mesh->ContactElements->Tag[sampleNo];
2145 break;
2146 case(DegreesOfFreedom):
2147 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2148 break;
2149 case(ReducedDegreesOfFreedom):
2150 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2151 break;
2152 default:
2153 stringstream temp;
2154 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2155 throw FinleyAdapterException(temp.str());
2156 break;
2157 }
2158 return out;
2159 }
2160
2161
2162 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2163 {
2164 Finley_Mesh* mesh=m_finleyMesh.get();
2165 escriptDataC tmp=mask.getDataC();
2166 switch(functionSpaceType) {
2167 case(Nodes):
2168 Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2169 break;
2170 case(ReducedNodes):
2171 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2172 break;
2173 case(DegreesOfFreedom):
2174 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2175 break;
2176 case(ReducedDegreesOfFreedom):
2177 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2178 break;
2179 case(Elements):
2180 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2181 break;
2182 case(ReducedElements):
2183 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2184 break;
2185 case(FaceElements):
2186 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2187 break;
2188 case(ReducedFaceElements):
2189 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2190 break;
2191 case(Points):
2192 Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2193 break;
2194 case(ContactElementsZero):
2195 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2196 break;
2197 case(ReducedContactElementsZero):
2198 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2199 break;
2200 case(ContactElementsOne):
2201 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2202 break;
2203 case(ReducedContactElementsOne):
2204 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2205 break;
2206 default:
2207 stringstream temp;
2208 temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2209 throw FinleyAdapterException(temp.str());
2210 }
2211 checkFinleyError();
2212 return;
2213 }
2214
2215 void MeshAdapter::setTagMap(const string& name, int tag)
2216 {
2217 Finley_Mesh* mesh=m_finleyMesh.get();
2218 Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2219 checkPasoError();
2220 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2221 }
2222
2223 int MeshAdapter::getTag(const string& name) const
2224 {
2225 Finley_Mesh* mesh=m_finleyMesh.get();
2226 int tag=0;
2227 tag=Finley_Mesh_getTag(mesh, name.c_str());
2228 checkPasoError();
2229 // throwStandardException("MeshAdapter::getTag is not implemented.");
2230 return tag;
2231 }
2232
2233 bool MeshAdapter::isValidTagName(const string& name) const
2234 {
2235 Finley_Mesh* mesh=m_finleyMesh.get();
2236 return Finley_Mesh_isValidTagName(mesh,name.c_str());
2237 }
2238
2239 string MeshAdapter::showTagNames() const
2240 {
2241 stringstream temp;
2242 Finley_Mesh* mesh=m_finleyMesh.get();
2243 Finley_TagMap* tag_map=mesh->TagMap;
2244 while (tag_map) {
2245 temp << tag_map->name;
2246 tag_map=tag_map->next;
2247 if (tag_map) temp << ", ";
2248 }
2249 return temp.str();
2250 }
2251
2252 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2253 {
2254 Finley_Mesh* mesh=m_finleyMesh.get();
2255 dim_t numTags=0;
2256 switch(functionSpaceCode) {
2257 case(Nodes):
2258 numTags=mesh->Nodes->numTagsInUse;
2259 break;
2260 case(ReducedNodes):
2261 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2262 break;
2263 case(DegreesOfFreedom):
2264 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2265 break;
2266 case(ReducedDegreesOfFreedom):
2267 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2268 break;
2269 case(Elements):
2270 case(ReducedElements):
2271 numTags=mesh->Elements->numTagsInUse;
2272 break;
2273 case(FaceElements):
2274 case(ReducedFaceElements):
2275 numTags=mesh->FaceElements->numTagsInUse;
2276 break;
2277 case(Points):
2278 numTags=mesh->Points->numTagsInUse;
2279 break;
2280 case(ContactElementsZero):
2281 case(ReducedContactElementsZero):
2282 case(ContactElementsOne):
2283 case(ReducedContactElementsOne):
2284 numTags=mesh->ContactElements->numTagsInUse;
2285 break;
2286 default:
2287 stringstream temp;
2288 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2289 throw FinleyAdapterException(temp.str());
2290 }
2291 return numTags;
2292 }
2293
2294 const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2295 {
2296 Finley_Mesh* mesh=m_finleyMesh.get();
2297 index_t* tags=NULL;
2298 switch(functionSpaceCode) {
2299 case(Nodes):
2300 tags=mesh->Nodes->tagsInUse;
2301 break;
2302 case(ReducedNodes):
2303 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2304 break;
2305 case(DegreesOfFreedom):
2306 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2307 break;
2308 case(ReducedDegreesOfFreedom):
2309 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2310 break;
2311 case(Elements):
2312 case(ReducedElements):
2313 tags=mesh->Elements->tagsInUse;
2314 break;
2315 case(FaceElements):
2316 case(ReducedFaceElements):
2317 tags=mesh->FaceElements->tagsInUse;
2318 break;
2319 case(Points):
2320 tags=mesh->Points->tagsInUse;
2321 break;
2322 case(ContactElementsZero):
2323 case(ReducedContactElementsZero):
2324 case(ContactElementsOne):
2325 case(ReducedContactElementsOne):
2326 tags=mesh->ContactElements->tagsInUse;
2327 break;
2328 default:
2329 stringstream temp;
2330 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2331 throw FinleyAdapterException(temp.str());
2332 }
2333 return tags;
2334 }
2335
2336
2337 bool MeshAdapter::canTag(int functionSpaceCode) const
2338 {
2339 switch(functionSpaceCode) {
2340 case(Nodes):
2341 case(Elements):
2342 case(ReducedElements):
2343 case(FaceElements):
2344 case(ReducedFaceElements):
2345 case(Points):
2346 case(ContactElementsZero):
2347 case(ReducedContactElementsZero):
2348 case(ContactElementsOne):
2349 case(ReducedContactElementsOne):
2350 return true;
2351 case(ReducedNodes):
2352 case(DegreesOfFreedom):
2353 case(ReducedDegreesOfFreedom):
2354 return false;
2355 default:
2356 return false;
2357 }
2358 }
2359
2360 AbstractDomain::StatusType MeshAdapter::getStatus() const
2361 {
2362 Finley_Mesh* mesh=m_finleyMesh.get();
2363 return Finley_Mesh_getStatus(mesh);
2364 }
2365
2366 int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2367 {
2368
2369 Finley_Mesh* mesh=m_finleyMesh.get();
2370 int order =-1;
2371 switch(functionSpaceCode) {
2372 case(Nodes):
2373 case(DegreesOfFreedom):
2374 order=mesh->approximationOrder;
2375 break;
2376 case(ReducedNodes):
2377 case(ReducedDegreesOfFreedom):
2378 order=mesh->reducedApproximationOrder;
2379 break;
2380 case(Elements):
2381 case(FaceElements):
2382 case(Points):
2383 case(ContactElementsZero):
2384 case(ContactElementsOne):
2385 order=mesh->integrationOrder;
2386 break;
2387 case(ReducedElements):
2388 case(ReducedFaceElements):
2389 case(ReducedContactElementsZero):
2390 case(ReducedContactElementsOne):
2391 order=mesh->reducedIntegrationOrder;
2392 break;
2393 default:
2394 stringstream temp;
2395 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2396 throw FinleyAdapterException(temp.str());
2397 }
2398 return order;
2399 }
2400
2401 bool MeshAdapter::supportsContactElements() const
2402 {
2403 return true;
2404 }
2405
2406 ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2407 {
2408 m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2409 }
2410
2411 ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2412 {
2413 Finley_ReferenceElementSet_dealloc(m_refSet);
2414 }
2415
2416 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26