/[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 3269 - (show annotations)
Wed Oct 13 03:21:50 2010 UTC (8 years, 11 months ago) by jfenwick
File size: 86093 byte(s)
Fixed some intel compiler warnings.
Put linearPDEs back the way it was and present a common interface for dudley and finley (as per Lutz)

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26