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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26