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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26