/[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 3317 - (show annotations)
Thu Oct 28 00:50:41 2010 UTC (8 years, 7 months ago) by caltinay
File size: 86021 byte(s)
Removed bogus mpi.h includes, replaced others by our esysUtils wrapper
and rearranged so that the wrapper always comes before netcdf which fixes
linking problems when disabling mpi on systems with netCDF 4.x.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26