/[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 3355 - (show annotations)
Tue Nov 16 06:35:06 2010 UTC (8 years, 2 months ago) by caltinay
File size: 85795 byte(s)
Fixed building with boost < 1.40 by
- removing the saveVTK method from dudley's MeshAdapter completely
  (no point introducing it now when we are trying to get rid of it soon)
- adding a helper function to weipa's python layer which can be called from C++
  with older boost versions.


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) const
809 {
810 escriptDataC _mat=mat.getDataC();
811 escriptDataC _D=D.getDataC();
812 escriptDataC _d=d.getDataC();
813
814 Finley_Mesh* mesh=m_finleyMesh.get();
815
816 Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
817 checkFinleyError();
818
819 Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
820 checkFinleyError();
821
822 }
823
824
825 //
826 // adds linear PDE of second order into the right hand side only
827 //
828 void MeshAdapter::addPDEToRHS( escript::Data& rhs, const escript::Data& X,const escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
829 {
830 Finley_Mesh* mesh=m_finleyMesh.get();
831
832 escriptDataC _rhs=rhs.getDataC();
833 escriptDataC _X=X.getDataC();
834 escriptDataC _Y=Y.getDataC();
835 escriptDataC _y=y.getDataC();
836 escriptDataC _y_contact=y_contact.getDataC();
837
838 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
839 checkFinleyError();
840
841 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
842 checkFinleyError();
843
844 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
845 checkFinleyError();
846 }
847 //
848 // adds PDE of second order into a transport problem
849 //
850 void MeshAdapter::addPDEToTransportProblem(
851 AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
852 const escript::Data& A, const escript::Data& B, const escript::Data& C,
853 const escript::Data& D,const escript::Data& X,const escript::Data& Y,
854 const escript::Data& d, const escript::Data& y,
855 const escript::Data& d_contact,const escript::Data& y_contact) const
856 {
857 TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
858 if (tpa==0)
859 {
860 throw FinleyAdapterException("finley only supports its own transport problems.");
861 }
862
863
864 DataTypes::ShapeType shape;
865 source.expand();
866 escriptDataC _source=source.getDataC();
867 escriptDataC _M=M.getDataC();
868 escriptDataC _A=A.getDataC();
869 escriptDataC _B=B.getDataC();
870 escriptDataC _C=C.getDataC();
871 escriptDataC _D=D.getDataC();
872 escriptDataC _X=X.getDataC();
873 escriptDataC _Y=Y.getDataC();
874 escriptDataC _d=d.getDataC();
875 escriptDataC _y=y.getDataC();
876 escriptDataC _d_contact=d_contact.getDataC();
877 escriptDataC _y_contact=y_contact.getDataC();
878
879 Finley_Mesh* mesh=m_finleyMesh.get();
880 Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
881
882 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
883 checkFinleyError();
884
885 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
886 checkFinleyError();
887
888 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
889 checkFinleyError();
890
891 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
892 checkFinleyError();
893 }
894
895 //
896 // interpolates data between different function spaces:
897 //
898 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
899 {
900 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
901 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
902 if (inDomain!=*this)
903 throw FinleyAdapterException("Error - Illegal domain of interpolant.");
904 if (targetDomain!=*this)
905 throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
906
907 Finley_Mesh* mesh=m_finleyMesh.get();
908 escriptDataC _target=target.getDataC();
909 escriptDataC _in=in.getDataC();
910 switch(in.getFunctionSpace().getTypeCode()) {
911 case(Nodes):
912 switch(target.getFunctionSpace().getTypeCode()) {
913 case(Nodes):
914 case(ReducedNodes):
915 case(DegreesOfFreedom):
916 case(ReducedDegreesOfFreedom):
917 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
918 break;
919 case(Elements):
920 case(ReducedElements):
921 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
922 break;
923 case(FaceElements):
924 case(ReducedFaceElements):
925 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
926 break;
927 case(Points):
928 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
929 break;
930 case(ContactElementsZero):
931 case(ReducedContactElementsZero):
932 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
933 break;
934 case(ContactElementsOne):
935 case(ReducedContactElementsOne):
936 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
937 break;
938 default:
939 stringstream temp;
940 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
941 throw FinleyAdapterException(temp.str());
942 break;
943 }
944 break;
945 case(ReducedNodes):
946 switch(target.getFunctionSpace().getTypeCode()) {
947 case(Nodes):
948 case(ReducedNodes):
949 case(DegreesOfFreedom):
950 case(ReducedDegreesOfFreedom):
951 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
952 break;
953 case(Elements):
954 case(ReducedElements):
955 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
956 break;
957 case(FaceElements):
958 case(ReducedFaceElements):
959 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
960 break;
961 case(Points):
962 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
963 break;
964 case(ContactElementsZero):
965 case(ReducedContactElementsZero):
966 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
967 break;
968 case(ContactElementsOne):
969 case(ReducedContactElementsOne):
970 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
971 break;
972 default:
973 stringstream temp;
974 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
975 throw FinleyAdapterException(temp.str());
976 break;
977 }
978 break;
979 case(Elements):
980 if (target.getFunctionSpace().getTypeCode()==Elements) {
981 Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
982 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
983 Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
984 } else {
985 throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
986 }
987 break;
988 case(ReducedElements):
989 if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
990 Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
991 } else {
992 throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
993 }
994 break;
995 case(FaceElements):
996 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
997 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
998 } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
999 Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1000 } else {
1001 throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1002 }
1003 break;
1004 case(ReducedFaceElements):
1005 if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1006 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1007 } else {
1008 throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
1009 }
1010 break;
1011 case(Points):
1012 if (target.getFunctionSpace().getTypeCode()==Points) {
1013 Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1014 } else {
1015 throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1016 }
1017 break;
1018 case(ContactElementsZero):
1019 case(ContactElementsOne):
1020 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1021 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1022 } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1023 Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1024 } else {
1025 throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1026 }
1027 break;
1028 case(ReducedContactElementsZero):
1029 case(ReducedContactElementsOne):
1030 if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1031 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1032 } else {
1033 throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1034 }
1035 break;
1036 case(DegreesOfFreedom):
1037 switch(target.getFunctionSpace().getTypeCode()) {
1038 case(ReducedDegreesOfFreedom):
1039 case(DegreesOfFreedom):
1040 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1041 break;
1042
1043 case(Nodes):
1044 case(ReducedNodes):
1045 if (getMPISize()>1) {
1046 escript::Data temp=escript::Data(in);
1047 temp.expand();
1048 escriptDataC _in2 = temp.getDataC();
1049 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1050 } else {
1051 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1052 }
1053 break;
1054 case(Elements):
1055 case(ReducedElements):
1056 if (getMPISize()>1) {
1057 escript::Data temp=escript::Data( in, continuousFunction(*this) );
1058 escriptDataC _in2 = temp.getDataC();
1059 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1060 } else {
1061 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1062 }
1063 break;
1064 case(FaceElements):
1065 case(ReducedFaceElements):
1066 if (getMPISize()>1) {
1067 escript::Data temp=escript::Data( in, continuousFunction(*this) );
1068 escriptDataC _in2 = temp.getDataC();
1069 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1070
1071 } else {
1072 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1073 }
1074 break;
1075 case(Points):
1076 if (getMPISize()>1) {
1077 escript::Data temp=escript::Data( in, continuousFunction(*this) );
1078 escriptDataC _in2 = temp.getDataC();
1079 } else {
1080 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1081 }
1082 break;
1083 case(ContactElementsZero):
1084 case(ContactElementsOne):
1085 case(ReducedContactElementsZero):
1086 case(ReducedContactElementsOne):
1087 if (getMPISize()>1) {
1088 escript::Data temp=escript::Data( in, continuousFunction(*this) );
1089 escriptDataC _in2 = temp.getDataC();
1090 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1091 } else {
1092 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1093 }
1094 break;
1095 default:
1096 stringstream temp;
1097 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1098 throw FinleyAdapterException(temp.str());
1099 break;
1100 }
1101 break;
1102 case(ReducedDegreesOfFreedom):
1103 switch(target.getFunctionSpace().getTypeCode()) {
1104 case(Nodes):
1105 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1106 break;
1107 case(ReducedNodes):
1108 if (getMPISize()>1) {
1109 escript::Data temp=escript::Data(in);
1110 temp.expand();
1111 escriptDataC _in2 = temp.getDataC();
1112 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1113 } else {
1114 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1115 }
1116 break;
1117 case(DegreesOfFreedom):
1118 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1119 break;
1120 case(ReducedDegreesOfFreedom):
1121 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1122 break;
1123 case(Elements):
1124 case(ReducedElements):
1125 if (getMPISize()>1) {
1126 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
1127 escriptDataC _in2 = temp.getDataC();
1128 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1129 } else {
1130 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1131 }
1132 break;
1133 case(FaceElements):
1134 case(ReducedFaceElements):
1135 if (getMPISize()>1) {
1136 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
1137 escriptDataC _in2 = temp.getDataC();
1138 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1139 } else {
1140 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1141 }
1142 break;
1143 case(Points):
1144 if (getMPISize()>1) {
1145 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
1146 escriptDataC _in2 = temp.getDataC();
1147 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1148 } else {
1149 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1150 }
1151 break;
1152 case(ContactElementsZero):
1153 case(ContactElementsOne):
1154 case(ReducedContactElementsZero):
1155 case(ReducedContactElementsOne):
1156 if (getMPISize()>1) {
1157 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
1158 escriptDataC _in2 = temp.getDataC();
1159 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1160 } else {
1161 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1162 }
1163 break;
1164 default:
1165 stringstream temp;
1166 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1167 throw FinleyAdapterException(temp.str());
1168 break;
1169 }
1170 break;
1171 default:
1172 stringstream temp;
1173 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1174 throw FinleyAdapterException(temp.str());
1175 break;
1176 }
1177 checkFinleyError();
1178 }
1179
1180 //
1181 // copies the locations of sample points into x:
1182 //
1183 void MeshAdapter::setToX(escript::Data& arg) const
1184 {
1185 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1186 if (argDomain!=*this)
1187 throw FinleyAdapterException("Error - Illegal domain of data point locations");
1188 Finley_Mesh* mesh=m_finleyMesh.get();
1189 // in case of values node coordinates we can do the job directly:
1190 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1191 escriptDataC _arg=arg.getDataC();
1192 Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1193 } else {
1194 escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1195 escriptDataC _tmp_data=tmp_data.getDataC();
1196 Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1197 // this is then interpolated onto arg:
1198 interpolateOnDomain(arg,tmp_data);
1199 }
1200 checkFinleyError();
1201 }
1202
1203 //
1204 // return the normal vectors at the location of data points as a Data object:
1205 //
1206 void MeshAdapter::setToNormal(escript::Data& normal) const
1207 {
1208 /* const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1209 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1210 if (normalDomain!=*this)
1211 throw FinleyAdapterException("Error - Illegal domain of normal locations");
1212 Finley_Mesh* mesh=m_finleyMesh.get();
1213 escriptDataC _normal=normal.getDataC();
1214 switch(normal.getFunctionSpace().getTypeCode()) {
1215 case(Nodes):
1216 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1217 break;
1218 case(ReducedNodes):
1219 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1220 break;
1221 case(Elements):
1222 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1223 break;
1224 case(ReducedElements):
1225 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1226 break;
1227 case (FaceElements):
1228 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1229 break;
1230 case (ReducedFaceElements):
1231 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1232 break;
1233 case(Points):
1234 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1235 break;
1236 case (ContactElementsOne):
1237 case (ContactElementsZero):
1238 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1239 break;
1240 case (ReducedContactElementsOne):
1241 case (ReducedContactElementsZero):
1242 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1243 break;
1244 case(DegreesOfFreedom):
1245 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1246 break;
1247 case(ReducedDegreesOfFreedom):
1248 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1249 break;
1250 default:
1251 stringstream temp;
1252 temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1253 throw FinleyAdapterException(temp.str());
1254 break;
1255 }
1256 checkFinleyError();
1257 }
1258
1259 //
1260 // interpolates data to other domain:
1261 //
1262 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1263 {
1264 const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1265 const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1266 if (targetDomain!=this)
1267 throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1268
1269 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
1270 }
1271
1272 //
1273 // calculates the integral of a function defined of arg:
1274 //
1275 void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1276 {
1277 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1278 if (argDomain!=*this)
1279 throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1280
1281 double blocktimer_start = blocktimer_time();
1282 Finley_Mesh* mesh=m_finleyMesh.get();
1283 escriptDataC _temp;
1284 escript::Data temp;
1285 escriptDataC _arg=arg.getDataC();
1286 switch(arg.getFunctionSpace().getTypeCode()) {
1287 case(Nodes):
1288 temp=escript::Data( arg, escript::function(*this) );
1289 _temp=temp.getDataC();
1290 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1291 break;
1292 case(ReducedNodes):
1293 temp=escript::Data( arg, escript::function(*this) );
1294 _temp=temp.getDataC();
1295 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1296 break;
1297 case(Elements):
1298 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1299 break;
1300 case(ReducedElements):
1301 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1302 break;
1303 case(FaceElements):
1304 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1305 break;
1306 case(ReducedFaceElements):
1307 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1308 break;
1309 case(Points):
1310 throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1311 break;
1312 case(ContactElementsZero):
1313 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1314 break;
1315 case(ReducedContactElementsZero):
1316 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1317 break;
1318 case(ContactElementsOne):
1319 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1320 break;
1321 case(ReducedContactElementsOne):
1322 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1323 break;
1324 case(DegreesOfFreedom):
1325 temp=escript::Data( arg, escript::function(*this) );
1326 _temp=temp.getDataC();
1327 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1328 break;
1329 case(ReducedDegreesOfFreedom):
1330 temp=escript::Data( arg, escript::function(*this) );
1331 _temp=temp.getDataC();
1332 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1333 break;
1334 default:
1335 stringstream temp;
1336 temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1337 throw FinleyAdapterException(temp.str());
1338 break;
1339 }
1340 checkFinleyError();
1341 blocktimer_increment("integrate()", blocktimer_start);
1342 }
1343
1344 //
1345 // calculates the gradient of arg:
1346 //
1347 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1348 {
1349 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1350 if (argDomain!=*this)
1351 throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1352 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1353 if (gradDomain!=*this)
1354 throw FinleyAdapterException("Error - Illegal domain of gradient");
1355
1356 Finley_Mesh* mesh=m_finleyMesh.get();
1357 escriptDataC _grad=grad.getDataC();
1358 escriptDataC nodeDataC;
1359 escript::Data temp;
1360 if (getMPISize()>1) {
1361 if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1362 temp=escript::Data( arg, continuousFunction(*this) );
1363 nodeDataC = temp.getDataC();
1364 } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1365 temp=escript::Data( arg, reducedContinuousFunction(*this) );
1366 nodeDataC = temp.getDataC();
1367 } else {
1368 nodeDataC = arg.getDataC();
1369 }
1370 } else {
1371 nodeDataC = arg.getDataC();
1372 }
1373 switch(grad.getFunctionSpace().getTypeCode()) {
1374 case(Nodes):
1375 throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1376 break;
1377 case(ReducedNodes):
1378 throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1379 break;
1380 case(Elements):
1381 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1382 break;
1383 case(ReducedElements):
1384 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1385 break;
1386 case(FaceElements):
1387 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1388 break;
1389 case(ReducedFaceElements):
1390 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1391 break;
1392 case(Points):
1393 throw FinleyAdapterException("Error - Gradient at points is not supported.");
1394 break;
1395 case(ContactElementsZero):
1396 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1397 break;
1398 case(ReducedContactElementsZero):
1399 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1400 break;
1401 case(ContactElementsOne):
1402 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1403 break;
1404 case(ReducedContactElementsOne):
1405 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1406 break;
1407 case(DegreesOfFreedom):
1408 throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1409 break;
1410 case(ReducedDegreesOfFreedom):
1411 throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1412 break;
1413 default:
1414 stringstream temp;
1415 temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1416 throw FinleyAdapterException(temp.str());
1417 break;
1418 }
1419 checkFinleyError();
1420 }
1421
1422 //
1423 // returns the size of elements:
1424 //
1425 void MeshAdapter::setToSize(escript::Data& size) const
1426 {
1427 Finley_Mesh* mesh=m_finleyMesh.get();
1428 escriptDataC tmp=size.getDataC();
1429 switch(size.getFunctionSpace().getTypeCode()) {
1430 case(Nodes):
1431 throw FinleyAdapterException("Error - Size of nodes is not supported.");
1432 break;
1433 case(ReducedNodes):
1434 throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1435 break;
1436 case(Elements):
1437 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1438 break;
1439 case(ReducedElements):
1440 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1441 break;
1442 case(FaceElements):
1443 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1444 break;
1445 case(ReducedFaceElements):
1446 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1447 break;
1448 case(Points):
1449 throw FinleyAdapterException("Error - Size of point elements is not supported.");
1450 break;
1451 case(ContactElementsZero):
1452 case(ContactElementsOne):
1453 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1454 break;
1455 case(ReducedContactElementsZero):
1456 case(ReducedContactElementsOne):
1457 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1458 break;
1459 case(DegreesOfFreedom):
1460 throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1461 break;
1462 case(ReducedDegreesOfFreedom):
1463 throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1464 break;
1465 default:
1466 stringstream temp;
1467 temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1468 throw FinleyAdapterException(temp.str());
1469 break;
1470 }
1471 checkFinleyError();
1472 }
1473
1474 //
1475 // sets the location of nodes
1476 //
1477 void MeshAdapter::setNewX(const escript::Data& new_x)
1478 {
1479 Finley_Mesh* mesh=m_finleyMesh.get();
1480 escriptDataC tmp;
1481 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1482 if (newDomain!=*this)
1483 throw FinleyAdapterException("Error - Illegal domain of new point locations");
1484 if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1485 tmp = new_x.getDataC();
1486 Finley_Mesh_setCoordinates(mesh,&tmp);
1487 } else {
1488 escript::Data new_x_inter=escript::Data( new_x, continuousFunction(*this) );
1489 tmp = new_x_inter.getDataC();
1490 Finley_Mesh_setCoordinates(mesh,&tmp);
1491 }
1492 checkFinleyError();
1493 }
1494
1495 //
1496 // Helper for the save* methods. Extracts optional data variable names and the
1497 // corresponding pointers from python dictionary. Caller must free arrays.
1498 //
1499 void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1500 {
1501 numData = boost::python::extract<int>(arg.attr("__len__")());
1502 /* win32 refactor */
1503 names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1504 data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1505 dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1506
1507 boost::python::list keys=arg.keys();
1508 for (int i=0; i<numData; ++i) {
1509 string n=boost::python::extract<string>(keys[i]);
1510 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1511 if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1512 throw FinleyAdapterException("Error: Data must be defined on same Domain");
1513 data[i] = d.getDataC();
1514 dataPtr[i] = &(data[i]);
1515 names[i] = TMPMEMALLOC(n.length()+1, char);
1516 strcpy(names[i], n.c_str());
1517 }
1518 }
1519
1520 //
1521 // saves mesh and optionally data arrays in openDX format
1522 //
1523 void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1524 {
1525 int num_data;
1526 char **names;
1527 escriptDataC *data;
1528 escriptDataC **ptr_data;
1529
1530 extractArgsFromDict(arg, num_data, names, data, ptr_data);
1531 Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1532 checkFinleyError();
1533
1534 /* win32 refactor */
1535 TMPMEMFREE(data);
1536 TMPMEMFREE(ptr_data);
1537 for(int i=0; i<num_data; i++)
1538 {
1539 TMPMEMFREE(names[i]);
1540 }
1541 TMPMEMFREE(names);
1542
1543 return;
1544 }
1545
1546 //
1547 // saves mesh and optionally data arrays in VTK format
1548 //
1549 void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg, const string& metadata, const string& metadata_schema) const
1550 {
1551 boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1552 pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1553 metadata, metadata_schema, arg);
1554 }
1555
1556 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1557 {
1558 #ifdef ESYS_MPI
1559 index_t myFirstNode=0, myLastNode=0, k=0;
1560 index_t* globalNodeIndex=0;
1561 Finley_Mesh* mesh_p=m_finleyMesh.get();
1562 if (fs_code == FINLEY_REDUCED_NODES)
1563 {
1564 myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1565 myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1566 globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1567 }
1568 else
1569 {
1570 myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1571 myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1572 globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1573 }
1574 k=globalNodeIndex[id];
1575 return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1576 #endif
1577 return true;
1578 }
1579
1580
1581
1582 //
1583 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1584 //
1585 ASM_ptr MeshAdapter::newSystemMatrix(
1586 const int row_blocksize,
1587 const escript::FunctionSpace& row_functionspace,
1588 const int column_blocksize,
1589 const escript::FunctionSpace& column_functionspace,
1590 const int type) const
1591 {
1592 int reduceRowOrder=0;
1593 int reduceColOrder=0;
1594 // is the domain right?
1595 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1596 if (row_domain!=*this)
1597 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1598 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1599 if (col_domain!=*this)
1600 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1601 // is the function space type right
1602 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1603 reduceRowOrder=0;
1604 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1605 reduceRowOrder=1;
1606 } else {
1607 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1608 }
1609 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1610 reduceColOrder=0;
1611 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1612 reduceColOrder=1;
1613 } else {
1614 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1615 }
1616 // generate matrix:
1617
1618 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1619 checkFinleyError();
1620 Paso_SystemMatrix* fsystemMatrix;
1621 int trilinos = 0;
1622 if (trilinos) {
1623 #ifdef TRILINOS
1624 /* Allocation Epetra_VrbMatrix here */
1625 #endif
1626 }
1627 else {
1628 fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1629 }
1630 checkPasoError();
1631 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1632 SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1633 return ASM_ptr(sma);
1634 // return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1635 }
1636
1637 //
1638 // creates a TransportProblemAdapter
1639 //
1640 ATP_ptr MeshAdapter::newTransportProblem(
1641 const bool useBackwardEuler,
1642 const int blocksize,
1643 const escript::FunctionSpace& functionspace,
1644 const int type) const
1645 {
1646 int reduceOrder=0;
1647 // is the domain right?
1648 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1649 if (domain!=*this)
1650 throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1651 // is the function space type right
1652 if (functionspace.getTypeCode()==DegreesOfFreedom) {
1653 reduceOrder=0;
1654 } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1655 reduceOrder=1;
1656 } else {
1657 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1658 }
1659 // generate matrix:
1660
1661 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1662 checkFinleyError();
1663 Paso_TransportProblem* transportProblem;
1664 transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1665 checkPasoError();
1666 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1667 TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1668 return ATP_ptr(tpa);
1669 // return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1670 }
1671
1672 //
1673 // vtkObject MeshAdapter::createVtkObject() const
1674 // TODO:
1675 //
1676
1677 //
1678 // returns true if data at the atom_type is considered as being cell centered:
1679 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1680 {
1681 switch(functionSpaceCode) {
1682 case(Nodes):
1683 case(DegreesOfFreedom):
1684 case(ReducedDegreesOfFreedom):
1685 return false;
1686 break;
1687 case(Elements):
1688 case(FaceElements):
1689 case(Points):
1690 case(ContactElementsZero):
1691 case(ContactElementsOne):
1692 case(ReducedElements):
1693 case(ReducedFaceElements):
1694 case(ReducedContactElementsZero):
1695 case(ReducedContactElementsOne):
1696 return true;
1697 break;
1698 default:
1699 stringstream temp;
1700 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1701 throw FinleyAdapterException(temp.str());
1702 break;
1703 }
1704 checkFinleyError();
1705 return false;
1706 }
1707
1708 bool
1709 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1710 {
1711 /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1712 class 1: DOF <-> Nodes
1713 class 2: ReducedDOF <-> ReducedNodes
1714 class 3: Points
1715 class 4: Elements
1716 class 5: ReducedElements
1717 class 6: FaceElements
1718 class 7: ReducedFaceElements
1719 class 8: ContactElementZero <-> ContactElementOne
1720 class 9: ReducedContactElementZero <-> ReducedContactElementOne
1721
1722 There is also a set of lines. Interpolation is possible down a line but not between lines.
1723 class 1 and 2 belong to all lines so aren't considered.
1724 line 0: class 3
1725 line 1: class 4,5
1726 line 2: class 6,7
1727 line 3: class 8,9
1728
1729 For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1730 eg hasnodes is true if we have at least one instance of Nodes.
1731 */
1732 if (fs.empty())
1733 {
1734 return false;
1735 }
1736 vector<int> hasclass(10);
1737 vector<int> hasline(4);
1738 bool hasnodes=false;
1739 bool hasrednodes=false;
1740 bool hascez=false;
1741 bool hasrcez=false;
1742 for (int i=0;i<fs.size();++i)
1743 {
1744 switch(fs[i])
1745 {
1746 case(Nodes): hasnodes=true; // no break is deliberate
1747 case(DegreesOfFreedom):
1748 hasclass[1]=1;
1749 break;
1750 case(ReducedNodes): hasrednodes=true; // no break is deliberate
1751 case(ReducedDegreesOfFreedom):
1752 hasclass[2]=1;
1753 break;
1754 case(Points):
1755 hasline[0]=1;
1756 hasclass[3]=1;
1757 break;
1758 case(Elements):
1759 hasclass[4]=1;
1760 hasline[1]=1;
1761 break;
1762 case(ReducedElements):
1763 hasclass[5]=1;
1764 hasline[1]=1;
1765 break;
1766 case(FaceElements):
1767 hasclass[6]=1;
1768 hasline[2]=1;
1769 break;
1770 case(ReducedFaceElements):
1771 hasclass[7]=1;
1772 hasline[2]=1;
1773 break;
1774 case(ContactElementsZero): hascez=true; // no break is deliberate
1775 case(ContactElementsOne):
1776 hasclass[8]=1;
1777 hasline[3]=1;
1778 break;
1779 case(ReducedContactElementsZero): hasrcez=true; // no break is deliberate
1780 case(ReducedContactElementsOne):
1781 hasclass[9]=1;
1782 hasline[3]=1;
1783 break;
1784 default:
1785 return false;
1786 }
1787 }
1788 int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1789 // fail if we have more than one leaf group
1790
1791 if (totlines>1)
1792 {
1793 return false; // there are at least two branches we can't interpolate between
1794 }
1795 else if (totlines==1)
1796 {
1797 if (hasline[0]==1) // we have points
1798 {
1799 resultcode=Points;
1800 }
1801 else if (hasline[1]==1)
1802 {
1803 if (hasclass[5]==1)
1804 {
1805 resultcode=ReducedElements;
1806 }
1807 else
1808 {
1809 resultcode=Elements;
1810 }
1811 }
1812 else if (hasline[2]==1)
1813 {
1814 if (hasclass[7]==1)
1815 {
1816 resultcode=ReducedFaceElements;
1817 }
1818 else
1819 {
1820 resultcode=FaceElements;
1821 }
1822 }
1823 else // so we must be in line3
1824 {
1825 if (hasclass[9]==1)
1826 {
1827 // need something from class 9
1828 resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1829 }
1830 else
1831 {
1832 // something from class 8
1833 resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1834 }
1835 }
1836 }
1837 else // totlines==0
1838 {
1839 if (hasclass[2]==1)
1840 {
1841 // something from class 2
1842 resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1843 }
1844 else
1845 { // something from class 1
1846 resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1847 }
1848 }
1849 return true;
1850 }
1851
1852 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1853 {
1854 switch(functionSpaceType_source) {
1855 case(Nodes):
1856 switch(functionSpaceType_target) {
1857 case(Nodes):
1858 case(ReducedNodes):
1859 case(ReducedDegreesOfFreedom):
1860 case(DegreesOfFreedom):
1861 case(Elements):
1862 case(ReducedElements):
1863 case(FaceElements):
1864 case(ReducedFaceElements):
1865 case(Points):
1866 case(ContactElementsZero):
1867 case(ReducedContactElementsZero):
1868 case(ContactElementsOne):
1869 case(ReducedContactElementsOne):
1870 return true;
1871 default:
1872 stringstream temp;
1873 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1874 throw FinleyAdapterException(temp.str());
1875 }
1876 break;
1877 case(ReducedNodes):
1878 switch(functionSpaceType_target) {
1879 case(ReducedNodes):
1880 case(ReducedDegreesOfFreedom):
1881 case(Elements):
1882 case(ReducedElements):
1883 case(FaceElements):
1884 case(ReducedFaceElements):
1885 case(Points):
1886 case(ContactElementsZero):
1887 case(ReducedContactElementsZero):
1888 case(ContactElementsOne):
1889 case(ReducedContactElementsOne):
1890 return true;
1891 case(Nodes):
1892 case(DegreesOfFreedom):
1893 return false;
1894 default:
1895 stringstream temp;
1896 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1897 throw FinleyAdapterException(temp.str());
1898 }
1899 break;
1900 case(Elements):
1901 if (functionSpaceType_target==Elements) {
1902 return true;
1903 } else if (functionSpaceType_target==ReducedElements) {
1904 return true;
1905 } else {
1906 return false;
1907 }
1908 case(ReducedElements):
1909 if (functionSpaceType_target==ReducedElements) {
1910 return true;
1911 } else {
1912 return false;
1913 }
1914 case(FaceElements):
1915 if (functionSpaceType_target==FaceElements) {
1916 return true;
1917 } else if (functionSpaceType_target==ReducedFaceElements) {
1918 return true;
1919 } else {
1920 return false;
1921 }
1922 case(ReducedFaceElements):
1923 if (functionSpaceType_target==ReducedFaceElements) {
1924 return true;
1925 } else {
1926 return false;
1927 }
1928 case(Points):
1929 if (functionSpaceType_target==Points) {
1930 return true;
1931 } else {
1932 return false;
1933 }
1934 case(ContactElementsZero):
1935 case(ContactElementsOne):
1936 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1937 return true;
1938 } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1939 return true;
1940 } else {
1941 return false;
1942 }
1943 case(ReducedContactElementsZero):
1944 case(ReducedContactElementsOne):
1945 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1946 return true;
1947 } else {
1948 return false;
1949 }
1950 case(DegreesOfFreedom):
1951 switch(functionSpaceType_target) {
1952 case(ReducedDegreesOfFreedom):
1953 case(DegreesOfFreedom):
1954 case(Nodes):
1955 case(ReducedNodes):
1956 case(Elements):
1957 case(ReducedElements):
1958 case(Points):
1959 case(FaceElements):
1960 case(ReducedFaceElements):
1961 case(ContactElementsZero):
1962 case(ReducedContactElementsZero):
1963 case(ContactElementsOne):
1964 case(ReducedContactElementsOne):
1965 return true;
1966 default:
1967 stringstream temp;
1968 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1969 throw FinleyAdapterException(temp.str());
1970 }
1971 break;
1972 case(ReducedDegreesOfFreedom):
1973 switch(functionSpaceType_target) {
1974 case(ReducedDegreesOfFreedom):
1975 case(ReducedNodes):
1976 case(Elements):
1977 case(ReducedElements):
1978 case(FaceElements):
1979 case(ReducedFaceElements):
1980 case(Points):
1981 case(ContactElementsZero):
1982 case(ReducedContactElementsZero):
1983 case(ContactElementsOne):
1984 case(ReducedContactElementsOne):
1985 return true;
1986 case(Nodes):
1987 case(DegreesOfFreedom):
1988 return false;
1989 default:
1990 stringstream temp;
1991 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1992 throw FinleyAdapterException(temp.str());
1993 }
1994 break;
1995 default:
1996 stringstream temp;
1997 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1998 throw FinleyAdapterException(temp.str());
1999 break;
2000 }
2001 checkFinleyError();
2002 return false;
2003 }
2004
2005 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
2006 {
2007 return false;
2008 }
2009
2010 bool MeshAdapter::operator==(const AbstractDomain& other) const
2011 {
2012 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2013 if (temp!=0) {
2014 return (m_finleyMesh==temp->m_finleyMesh);
2015 } else {
2016 return false;
2017 }
2018 }
2019
2020 bool MeshAdapter::operator!=(const AbstractDomain& other) const
2021 {
2022 return !(operator==(other));
2023 }
2024
2025 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2026 {
2027 Finley_Mesh* mesh=m_finleyMesh.get();
2028 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2029 checkPasoError();
2030 return out;
2031 }
2032 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2033 {
2034 Finley_Mesh* mesh=m_finleyMesh.get();
2035 int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2036 checkPasoError();
2037 return out;
2038 }
2039
2040 escript::Data MeshAdapter::getX() const
2041 {
2042 return continuousFunction(*this).getX();
2043 }
2044
2045 escript::Data MeshAdapter::getNormal() const
2046 {
2047 return functionOnBoundary(*this).getNormal();
2048 }
2049
2050 escript::Data MeshAdapter::getSize() const
2051 {
2052 return escript::function(*this).getSize();
2053 }
2054
2055 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2056 {
2057 int *out = NULL;
2058 Finley_Mesh* mesh=m_finleyMesh.get();
2059 switch (functionSpaceType) {
2060 case(Nodes):
2061 out=mesh->Nodes->Id;
2062 break;
2063 case(ReducedNodes):
2064 out=mesh->Nodes->reducedNodesId;
2065 break;
2066 case(Elements):
2067 out=mesh->Elements->Id;
2068 break;
2069 case(ReducedElements):
2070 out=mesh->Elements->Id;
2071 break;
2072 case(FaceElements):
2073 out=mesh->FaceElements->Id;
2074 break;
2075 case(ReducedFaceElements):
2076 out=mesh->FaceElements->Id;
2077 break;
2078 case(Points):
2079 out=mesh->Points->Id;
2080 break;
2081 case(ContactElementsZero):
2082 out=mesh->ContactElements->Id;
2083 break;
2084 case(ReducedContactElementsZero):
2085 out=mesh->ContactElements->Id;
2086 break;
2087 case(ContactElementsOne):
2088 out=mesh->ContactElements->Id;
2089 break;
2090 case(ReducedContactElementsOne):
2091 out=mesh->ContactElements->Id;
2092 break;
2093 case(DegreesOfFreedom):
2094 out=mesh->Nodes->degreesOfFreedomId;
2095 break;
2096 case(ReducedDegreesOfFreedom):
2097 out=mesh->Nodes->reducedDegreesOfFreedomId;
2098 break;
2099 default:
2100 stringstream temp;
2101 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2102 throw FinleyAdapterException(temp.str());
2103 break;
2104 }
2105 return out;
2106 }
2107 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2108 {
2109 int out=0;
2110 Finley_Mesh* mesh=m_finleyMesh.get();
2111 switch (functionSpaceType) {
2112 case(Nodes):
2113 out=mesh->Nodes->Tag[sampleNo];
2114 break;
2115 case(ReducedNodes):
2116 throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2117 break;
2118 case(Elements):
2119 out=mesh->Elements->Tag[sampleNo];
2120 break;
2121 case(ReducedElements):
2122 out=mesh->Elements->Tag[sampleNo];
2123 break;
2124 case(FaceElements):
2125 out=mesh->FaceElements->Tag[sampleNo];
2126 break;
2127 case(ReducedFaceElements):
2128 out=mesh->FaceElements->Tag[sampleNo];
2129 break;
2130 case(Points):
2131 out=mesh->Points->Tag[sampleNo];
2132 break;
2133 case(ContactElementsZero):
2134 out=mesh->ContactElements->Tag[sampleNo];
2135 break;
2136 case(ReducedContactElementsZero):
2137 out=mesh->ContactElements->Tag[sampleNo];
2138 break;
2139 case(ContactElementsOne):
2140 out=mesh->ContactElements->Tag[sampleNo];
2141 break;
2142 case(ReducedContactElementsOne):
2143 out=mesh->ContactElements->Tag[sampleNo];
2144 break;
2145 case(DegreesOfFreedom):
2146 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2147 break;
2148 case(ReducedDegreesOfFreedom):
2149 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2150 break;
2151 default:
2152 stringstream temp;
2153 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2154 throw FinleyAdapterException(temp.str());
2155 break;
2156 }
2157 return out;
2158 }
2159
2160
2161 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2162 {
2163 Finley_Mesh* mesh=m_finleyMesh.get();
2164 escriptDataC tmp=mask.getDataC();
2165 switch(functionSpaceType) {
2166 case(Nodes):
2167 Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2168 break;
2169 case(ReducedNodes):
2170 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2171 break;
2172 case(DegreesOfFreedom):
2173 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2174 break;
2175 case(ReducedDegreesOfFreedom):
2176 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2177 break;
2178 case(Elements):
2179 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2180 break;
2181 case(ReducedElements):
2182 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2183 break;
2184 case(FaceElements):
2185 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2186 break;
2187 case(ReducedFaceElements):
2188 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2189 break;
2190 case(Points):
2191 Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2192 break;
2193 case(ContactElementsZero):
2194 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2195 break;
2196 case(ReducedContactElementsZero):
2197 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2198 break;
2199 case(ContactElementsOne):
2200 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2201 break;
2202 case(ReducedContactElementsOne):
2203 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2204 break;
2205 default:
2206 stringstream temp;
2207 temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2208 throw FinleyAdapterException(temp.str());
2209 }
2210 checkFinleyError();
2211 return;
2212 }
2213
2214 void MeshAdapter::setTagMap(const string& name, int tag)
2215 {
2216 Finley_Mesh* mesh=m_finleyMesh.get();
2217 Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2218 checkPasoError();
2219 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2220 }
2221
2222 int MeshAdapter::getTag(const string& name) const
2223 {
2224 Finley_Mesh* mesh=m_finleyMesh.get();
2225 int tag=0;
2226 tag=Finley_Mesh_getTag(mesh, name.c_str());
2227 checkPasoError();
2228 // throwStandardException("MeshAdapter::getTag is not implemented.");
2229 return tag;
2230 }
2231
2232 bool MeshAdapter::isValidTagName(const string& name) const
2233 {
2234 Finley_Mesh* mesh=m_finleyMesh.get();
2235 return Finley_Mesh_isValidTagName(mesh,name.c_str());
2236 }
2237
2238 string MeshAdapter::showTagNames() const
2239 {
2240 stringstream temp;
2241 Finley_Mesh* mesh=m_finleyMesh.get();
2242 Finley_TagMap* tag_map=mesh->TagMap;
2243 while (tag_map) {
2244 temp << tag_map->name;
2245 tag_map=tag_map->next;
2246 if (tag_map) temp << ", ";
2247 }
2248 return temp.str();
2249 }
2250
2251 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2252 {
2253 Finley_Mesh* mesh=m_finleyMesh.get();
2254 dim_t numTags=0;
2255 switch(functionSpaceCode) {
2256 case(Nodes):
2257 numTags=mesh->Nodes->numTagsInUse;
2258 break;
2259 case(ReducedNodes):
2260 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2261 break;
2262 case(DegreesOfFreedom):
2263 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2264 break;
2265 case(ReducedDegreesOfFreedom):
2266 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2267 break;
2268 case(Elements):
2269 case(ReducedElements):
2270 numTags=mesh->Elements->numTagsInUse;
2271 break;
2272 case(FaceElements):
2273 case(ReducedFaceElements):
2274 numTags=mesh->FaceElements->numTagsInUse;
2275 break;
2276 case(Points):
2277 numTags=mesh->Points->numTagsInUse;
2278 break;
2279 case(ContactElementsZero):
2280 case(ReducedContactElementsZero):
2281 case(ContactElementsOne):
2282 case(ReducedContactElementsOne):
2283 numTags=mesh->ContactElements->numTagsInUse;
2284 break;
2285 default:
2286 stringstream temp;
2287 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2288 throw FinleyAdapterException(temp.str());
2289 }
2290 return numTags;
2291 }
2292
2293 const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2294 {
2295 Finley_Mesh* mesh=m_finleyMesh.get();
2296 index_t* tags=NULL;
2297 switch(functionSpaceCode) {
2298 case(Nodes):
2299 tags=mesh->Nodes->tagsInUse;
2300 break;
2301 case(ReducedNodes):
2302 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2303 break;
2304 case(DegreesOfFreedom):
2305 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2306 break;
2307 case(ReducedDegreesOfFreedom):
2308 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2309 break;
2310 case(Elements):
2311 case(ReducedElements):
2312 tags=mesh->Elements->tagsInUse;
2313 break;
2314 case(FaceElements):
2315 case(ReducedFaceElements):
2316 tags=mesh->FaceElements->tagsInUse;
2317 break;
2318 case(Points):
2319 tags=mesh->Points->tagsInUse;
2320 break;
2321 case(ContactElementsZero):
2322 case(ReducedContactElementsZero):
2323 case(ContactElementsOne):
2324 case(ReducedContactElementsOne):
2325 tags=mesh->ContactElements->tagsInUse;
2326 break;
2327 default:
2328 stringstream temp;
2329 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2330 throw FinleyAdapterException(temp.str());
2331 }
2332 return tags;
2333 }
2334
2335
2336 bool MeshAdapter::canTag(int functionSpaceCode) const
2337 {
2338 switch(functionSpaceCode) {
2339 case(Nodes):
2340 case(Elements):
2341 case(ReducedElements):
2342 case(FaceElements):
2343 case(ReducedFaceElements):
2344 case(Points):
2345 case(ContactElementsZero):
2346 case(ReducedContactElementsZero):
2347 case(ContactElementsOne):
2348 case(ReducedContactElementsOne):
2349 return true;
2350 case(ReducedNodes):
2351 case(DegreesOfFreedom):
2352 case(ReducedDegreesOfFreedom):
2353 return false;
2354 default:
2355 return false;
2356 }
2357 }
2358
2359 AbstractDomain::StatusType MeshAdapter::getStatus() const
2360 {
2361 Finley_Mesh* mesh=m_finleyMesh.get();
2362 return Finley_Mesh_getStatus(mesh);
2363 }
2364
2365 int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2366 {
2367
2368 Finley_Mesh* mesh=m_finleyMesh.get();
2369 int order =-1;
2370 switch(functionSpaceCode) {
2371 case(Nodes):
2372 case(DegreesOfFreedom):
2373 order=mesh->approximationOrder;
2374 break;
2375 case(ReducedNodes):
2376 case(ReducedDegreesOfFreedom):
2377 order=mesh->reducedApproximationOrder;
2378 break;
2379 case(Elements):
2380 case(FaceElements):
2381 case(Points):
2382 case(ContactElementsZero):
2383 case(ContactElementsOne):
2384 order=mesh->integrationOrder;
2385 break;
2386 case(ReducedElements):
2387 case(ReducedFaceElements):
2388 case(ReducedContactElementsZero):
2389 case(ReducedContactElementsOne):
2390 order=mesh->reducedIntegrationOrder;
2391 break;
2392 default:
2393 stringstream temp;
2394 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2395 throw FinleyAdapterException(temp.str());
2396 }
2397 return order;
2398 }
2399
2400 bool MeshAdapter::supportsContactElements() const
2401 {
2402 return true;
2403 }
2404
2405 ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2406 {
2407 m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2408 }
2409
2410 ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2411 {
2412 Finley_ReferenceElementSet_dealloc(m_refSet);
2413 }
2414
2415 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26