/[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 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years ago) by jfenwick
File size: 85719 byte(s)
Merging dudley and scons updates from branches

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26