/[escript]/branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Contents of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2989 - (show annotations)
Thu Mar 18 01:45:55 2010 UTC (9 years, 1 month ago) by gross
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 85811 byte(s)
- bug in error handling in the case of lumped mass matrix fixed.
- switched back to rowsum lumping as there are problems with HRZ lumping for order 2 


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 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 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(asAbstractContinuousDomain()) );
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(asAbstractContinuousDomain()) );
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(asAbstractContinuousDomain()) );
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(asAbstractContinuousDomain()) );
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(asAbstractContinuousDomain()) );
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(asAbstractContinuousDomain()) );
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(asAbstractContinuousDomain()) );
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(asAbstractContinuousDomain()) );
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(asAbstractContinuousDomain()),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(asAbstractContinuousDomain()) );
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(asAbstractContinuousDomain()) );
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(asAbstractContinuousDomain()) );
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(asAbstractContinuousDomain()) );
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(asAbstractContinuousDomain()) );
1353 nodeDataC = temp.getDataC();
1354 } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1355 temp=escript::Data( arg, reducedContinuousFunction(asAbstractContinuousDomain()) );
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(asAbstractContinuousDomain()) ) {
1475 tmp = new_x.getDataC();
1476 Finley_Mesh_setCoordinates(mesh,&tmp);
1477 } else {
1478 escript::Data new_x_inter=escript::Data( new_x, continuousFunction(asAbstractContinuousDomain()) );
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 PASO_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 SystemMatrixAdapter 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 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1637 }
1638
1639 //
1640 // creates a TransportProblemAdapter
1641 //
1642 TransportProblemAdapter MeshAdapter::newTransportProblem(
1643 const bool useBackwardEuler,
1644 const int blocksize,
1645 const escript::FunctionSpace& functionspace,
1646 const int type) const
1647 {
1648 int reduceOrder=0;
1649 // is the domain right?
1650 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1651 if (domain!=*this)
1652 throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1653 // is the function space type right
1654 if (functionspace.getTypeCode()==DegreesOfFreedom) {
1655 reduceOrder=0;
1656 } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1657 reduceOrder=1;
1658 } else {
1659 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1660 }
1661 // generate matrix:
1662
1663 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1664 checkFinleyError();
1665 Paso_TransportProblem* transportProblem;
1666 transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1667 checkPasoError();
1668 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1669 return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1670 }
1671
1672 //
1673 // vtkObject MeshAdapter::createVtkObject() const
1674 // TODO:
1675 //
1676
1677 //
1678 // returns true if data at the atom_type is considered as being cell centered:
1679 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1680 {
1681 switch(functionSpaceCode) {
1682 case(Nodes):
1683 case(DegreesOfFreedom):
1684 case(ReducedDegreesOfFreedom):
1685 return false;
1686 break;
1687 case(Elements):
1688 case(FaceElements):
1689 case(Points):
1690 case(ContactElementsZero):
1691 case(ContactElementsOne):
1692 case(ReducedElements):
1693 case(ReducedFaceElements):
1694 case(ReducedContactElementsZero):
1695 case(ReducedContactElementsOne):
1696 return true;
1697 break;
1698 default:
1699 stringstream temp;
1700 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1701 throw FinleyAdapterException(temp.str());
1702 break;
1703 }
1704 checkFinleyError();
1705 return false;
1706 }
1707
1708 bool
1709 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1710 {
1711 /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1712 class 1: DOF <-> Nodes
1713 class 2: ReducedDOF <-> ReducedNodes
1714 class 3: Points
1715 class 4: Elements
1716 class 5: ReducedElements
1717 class 6: FaceElements
1718 class 7: ReducedFaceElements
1719 class 8: ContactElementZero <-> ContactElementOne
1720 class 9: ReducedContactElementZero <-> ReducedContactElementOne
1721
1722 There is also a set of lines. Interpolation is possible down a line but not between lines.
1723 class 1 and 2 belong to all lines so aren't considered.
1724 line 0: class 3
1725 line 1: class 4,5
1726 line 2: class 6,7
1727 line 3: class 8,9
1728
1729 For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1730 eg hasnodes is true if we have at least one instance of Nodes.
1731 */
1732 if (fs.empty())
1733 {
1734 return false;
1735 }
1736 vector<int> hasclass(10);
1737 vector<int> hasline(4);
1738 bool hasnodes=false;
1739 bool hasrednodes=false;
1740 bool hascez=false;
1741 bool hasrcez=false;
1742 for (int i=0;i<fs.size();++i)
1743 {
1744 switch(fs[i])
1745 {
1746 case(Nodes): hasnodes=true; // no break is deliberate
1747 case(DegreesOfFreedom):
1748 hasclass[1]=1;
1749 break;
1750 case(ReducedNodes): hasrednodes=true; // no break is deliberate
1751 case(ReducedDegreesOfFreedom):
1752 hasclass[2]=1;
1753 break;
1754 case(Points):
1755 hasline[0]=1;
1756 hasclass[3]=1;
1757 break;
1758 case(Elements):
1759 hasclass[4]=1;
1760 hasline[1]=1;
1761 break;
1762 case(ReducedElements):
1763 hasclass[5]=1;
1764 hasline[1]=1;
1765 break;
1766 case(FaceElements):
1767 hasclass[6]=1;
1768 hasline[2]=1;
1769 break;
1770 case(ReducedFaceElements):
1771 hasclass[7]=1;
1772 hasline[2]=1;
1773 break;
1774 case(ContactElementsZero): hascez=true; // no break is deliberate
1775 case(ContactElementsOne):
1776 hasclass[8]=1;
1777 hasline[3]=1;
1778 break;
1779 case(ReducedContactElementsZero): hasrcez=true; // no break is deliberate
1780 case(ReducedContactElementsOne):
1781 hasclass[9]=1;
1782 hasline[3]=1;
1783 break;
1784 default:
1785 return false;
1786 }
1787 }
1788 int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1789 // fail if we have more than one leaf group
1790
1791 if (totlines>1)
1792 {
1793 return false; // there are at least two branches we can't interpolate between
1794 }
1795 else if (totlines==1)
1796 {
1797 if (hasline[0]==1) // we have points
1798 {
1799 resultcode=Points;
1800 }
1801 else if (hasline[1]==1)
1802 {
1803 if (hasclass[5]==1)
1804 {
1805 resultcode=ReducedElements;
1806 }
1807 else
1808 {
1809 resultcode=Elements;
1810 }
1811 }
1812 else if (hasline[2]==1)
1813 {
1814 if (hasclass[7]==1)
1815 {
1816 resultcode=ReducedFaceElements;
1817 }
1818 else
1819 {
1820 resultcode=FaceElements;
1821 }
1822 }
1823 else // so we must be in line3
1824 {
1825 if (hasclass[9]==1)
1826 {
1827 // need something from class 9
1828 resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1829 }
1830 else
1831 {
1832 // something from class 8
1833 resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1834 }
1835 }
1836 }
1837 else // totlines==0
1838 {
1839 if (hasclass[2]==1)
1840 {
1841 // something from class 2
1842 resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1843 }
1844 else
1845 { // something from class 1
1846 resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1847 }
1848 }
1849 return true;
1850 }
1851
1852 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1853 {
1854 switch(functionSpaceType_source) {
1855 case(Nodes):
1856 switch(functionSpaceType_target) {
1857 case(Nodes):
1858 case(ReducedNodes):
1859 case(ReducedDegreesOfFreedom):
1860 case(DegreesOfFreedom):
1861 case(Elements):
1862 case(ReducedElements):
1863 case(FaceElements):
1864 case(ReducedFaceElements):
1865 case(Points):
1866 case(ContactElementsZero):
1867 case(ReducedContactElementsZero):
1868 case(ContactElementsOne):
1869 case(ReducedContactElementsOne):
1870 return true;
1871 default:
1872 stringstream temp;
1873 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1874 throw FinleyAdapterException(temp.str());
1875 }
1876 break;
1877 case(ReducedNodes):
1878 switch(functionSpaceType_target) {
1879 case(ReducedNodes):
1880 case(ReducedDegreesOfFreedom):
1881 case(Elements):
1882 case(ReducedElements):
1883 case(FaceElements):
1884 case(ReducedFaceElements):
1885 case(Points):
1886 case(ContactElementsZero):
1887 case(ReducedContactElementsZero):
1888 case(ContactElementsOne):
1889 case(ReducedContactElementsOne):
1890 return true;
1891 case(Nodes):
1892 case(DegreesOfFreedom):
1893 return false;
1894 default:
1895 stringstream temp;
1896 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1897 throw FinleyAdapterException(temp.str());
1898 }
1899 break;
1900 case(Elements):
1901 if (functionSpaceType_target==Elements) {
1902 return true;
1903 } else if (functionSpaceType_target==ReducedElements) {
1904 return true;
1905 } else {
1906 return false;
1907 }
1908 case(ReducedElements):
1909 if (functionSpaceType_target==ReducedElements) {
1910 return true;
1911 } else {
1912 return false;
1913 }
1914 case(FaceElements):
1915 if (functionSpaceType_target==FaceElements) {
1916 return true;
1917 } else if (functionSpaceType_target==ReducedFaceElements) {
1918 return true;
1919 } else {
1920 return false;
1921 }
1922 case(ReducedFaceElements):
1923 if (functionSpaceType_target==ReducedFaceElements) {
1924 return true;
1925 } else {
1926 return false;
1927 }
1928 case(Points):
1929 if (functionSpaceType_target==Points) {
1930 return true;
1931 } else {
1932 return false;
1933 }
1934 case(ContactElementsZero):
1935 case(ContactElementsOne):
1936 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1937 return true;
1938 } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1939 return true;
1940 } else {
1941 return false;
1942 }
1943 case(ReducedContactElementsZero):
1944 case(ReducedContactElementsOne):
1945 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1946 return true;
1947 } else {
1948 return false;
1949 }
1950 case(DegreesOfFreedom):
1951 switch(functionSpaceType_target) {
1952 case(ReducedDegreesOfFreedom):
1953 case(DegreesOfFreedom):
1954 case(Nodes):
1955 case(ReducedNodes):
1956 case(Elements):
1957 case(ReducedElements):
1958 case(Points):
1959 case(FaceElements):
1960 case(ReducedFaceElements):
1961 case(ContactElementsZero):
1962 case(ReducedContactElementsZero):
1963 case(ContactElementsOne):
1964 case(ReducedContactElementsOne):
1965 return true;
1966 default:
1967 stringstream temp;
1968 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1969 throw FinleyAdapterException(temp.str());
1970 }
1971 break;
1972 case(ReducedDegreesOfFreedom):
1973 switch(functionSpaceType_target) {
1974 case(ReducedDegreesOfFreedom):
1975 case(ReducedNodes):
1976 case(Elements):
1977 case(ReducedElements):
1978 case(FaceElements):
1979 case(ReducedFaceElements):
1980 case(Points):
1981 case(ContactElementsZero):
1982 case(ReducedContactElementsZero):
1983 case(ContactElementsOne):
1984 case(ReducedContactElementsOne):
1985 return true;
1986 case(Nodes):
1987 case(DegreesOfFreedom):
1988 return false;
1989 default:
1990 stringstream temp;
1991 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1992 throw FinleyAdapterException(temp.str());
1993 }
1994 break;
1995 default:
1996 stringstream temp;
1997 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1998 throw FinleyAdapterException(temp.str());
1999 break;
2000 }
2001 checkFinleyError();
2002 return false;
2003 }
2004
2005 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
2006 {
2007 return false;
2008 }
2009
2010 bool MeshAdapter::operator==(const AbstractDomain& other) const
2011 {
2012 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2013 if (temp!=0) {
2014 return (m_finleyMesh==temp->m_finleyMesh);
2015 } else {
2016 return false;
2017 }
2018 }
2019
2020 bool MeshAdapter::operator!=(const AbstractDomain& other) const
2021 {
2022 return !(operator==(other));
2023 }
2024
2025 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2026 {
2027 Finley_Mesh* mesh=m_finleyMesh.get();
2028 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2029 checkPasoError();
2030 return out;
2031 }
2032 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2033 {
2034 Finley_Mesh* mesh=m_finleyMesh.get();
2035 int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2036 checkPasoError();
2037 return out;
2038 }
2039
2040 escript::Data MeshAdapter::getX() const
2041 {
2042 return continuousFunction(asAbstractContinuousDomain()).getX();
2043 }
2044
2045 escript::Data MeshAdapter::getNormal() const
2046 {
2047 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
2048 }
2049
2050 escript::Data MeshAdapter::getSize() const
2051 {
2052 return escript::function(asAbstractContinuousDomain()).getSize();
2053 }
2054
2055 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2056 {
2057 int *out = NULL;
2058 Finley_Mesh* mesh=m_finleyMesh.get();
2059 switch (functionSpaceType) {
2060 case(Nodes):
2061 out=mesh->Nodes->Id;
2062 break;
2063 case(ReducedNodes):
2064 out=mesh->Nodes->reducedNodesId;
2065 break;
2066 case(Elements):
2067 out=mesh->Elements->Id;
2068 break;
2069 case(ReducedElements):
2070 out=mesh->Elements->Id;
2071 break;
2072 case(FaceElements):
2073 out=mesh->FaceElements->Id;
2074 break;
2075 case(ReducedFaceElements):
2076 out=mesh->FaceElements->Id;
2077 break;
2078 case(Points):
2079 out=mesh->Points->Id;
2080 break;
2081 case(ContactElementsZero):
2082 out=mesh->ContactElements->Id;
2083 break;
2084 case(ReducedContactElementsZero):
2085 out=mesh->ContactElements->Id;
2086 break;
2087 case(ContactElementsOne):
2088 out=mesh->ContactElements->Id;
2089 break;
2090 case(ReducedContactElementsOne):
2091 out=mesh->ContactElements->Id;
2092 break;
2093 case(DegreesOfFreedom):
2094 out=mesh->Nodes->degreesOfFreedomId;
2095 break;
2096 case(ReducedDegreesOfFreedom):
2097 out=mesh->Nodes->reducedDegreesOfFreedomId;
2098 break;
2099 default:
2100 stringstream temp;
2101 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2102 throw FinleyAdapterException(temp.str());
2103 break;
2104 }
2105 return out;
2106 }
2107 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2108 {
2109 int out=0;
2110 Finley_Mesh* mesh=m_finleyMesh.get();
2111 switch (functionSpaceType) {
2112 case(Nodes):
2113 out=mesh->Nodes->Tag[sampleNo];
2114 break;
2115 case(ReducedNodes):
2116 throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2117 break;
2118 case(Elements):
2119 out=mesh->Elements->Tag[sampleNo];
2120 break;
2121 case(ReducedElements):
2122 out=mesh->Elements->Tag[sampleNo];
2123 break;
2124 case(FaceElements):
2125 out=mesh->FaceElements->Tag[sampleNo];
2126 break;
2127 case(ReducedFaceElements):
2128 out=mesh->FaceElements->Tag[sampleNo];
2129 break;
2130 case(Points):
2131 out=mesh->Points->Tag[sampleNo];
2132 break;
2133 case(ContactElementsZero):
2134 out=mesh->ContactElements->Tag[sampleNo];
2135 break;
2136 case(ReducedContactElementsZero):
2137 out=mesh->ContactElements->Tag[sampleNo];
2138 break;
2139 case(ContactElementsOne):
2140 out=mesh->ContactElements->Tag[sampleNo];
2141 break;
2142 case(ReducedContactElementsOne):
2143 out=mesh->ContactElements->Tag[sampleNo];
2144 break;
2145 case(DegreesOfFreedom):
2146 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2147 break;
2148 case(ReducedDegreesOfFreedom):
2149 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2150 break;
2151 default:
2152 stringstream temp;
2153 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2154 throw FinleyAdapterException(temp.str());
2155 break;
2156 }
2157 return out;
2158 }
2159
2160
2161 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2162 {
2163 Finley_Mesh* mesh=m_finleyMesh.get();
2164 escriptDataC tmp=mask.getDataC();
2165 switch(functionSpaceType) {
2166 case(Nodes):
2167 Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2168 break;
2169 case(ReducedNodes):
2170 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2171 break;
2172 case(DegreesOfFreedom):
2173 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2174 break;
2175 case(ReducedDegreesOfFreedom):
2176 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2177 break;
2178 case(Elements):
2179 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2180 break;
2181 case(ReducedElements):
2182 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2183 break;
2184 case(FaceElements):
2185 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2186 break;
2187 case(ReducedFaceElements):
2188 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2189 break;
2190 case(Points):
2191 Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2192 break;
2193 case(ContactElementsZero):
2194 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2195 break;
2196 case(ReducedContactElementsZero):
2197 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2198 break;
2199 case(ContactElementsOne):
2200 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2201 break;
2202 case(ReducedContactElementsOne):
2203 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2204 break;
2205 default:
2206 stringstream temp;
2207 temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2208 throw FinleyAdapterException(temp.str());
2209 }
2210 checkFinleyError();
2211 return;
2212 }
2213
2214 void MeshAdapter::setTagMap(const string& name, int tag)
2215 {
2216 Finley_Mesh* mesh=m_finleyMesh.get();
2217 Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2218 checkPasoError();
2219 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2220 }
2221
2222 int MeshAdapter::getTag(const string& name) const
2223 {
2224 Finley_Mesh* mesh=m_finleyMesh.get();
2225 int tag=0;
2226 tag=Finley_Mesh_getTag(mesh, name.c_str());
2227 checkPasoError();
2228 // throwStandardException("MeshAdapter::getTag is not implemented.");
2229 return tag;
2230 }
2231
2232 bool MeshAdapter::isValidTagName(const string& name) const
2233 {
2234 Finley_Mesh* mesh=m_finleyMesh.get();
2235 return Finley_Mesh_isValidTagName(mesh,name.c_str());
2236 }
2237
2238 string MeshAdapter::showTagNames() const
2239 {
2240 stringstream temp;
2241 Finley_Mesh* mesh=m_finleyMesh.get();
2242 Finley_TagMap* tag_map=mesh->TagMap;
2243 while (tag_map) {
2244 temp << tag_map->name;
2245 tag_map=tag_map->next;
2246 if (tag_map) temp << ", ";
2247 }
2248 return temp.str();
2249 }
2250
2251 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2252 {
2253 Finley_Mesh* mesh=m_finleyMesh.get();
2254 dim_t numTags=0;
2255 switch(functionSpaceCode) {
2256 case(Nodes):
2257 numTags=mesh->Nodes->numTagsInUse;
2258 break;
2259 case(ReducedNodes):
2260 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2261 break;
2262 case(DegreesOfFreedom):
2263 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2264 break;
2265 case(ReducedDegreesOfFreedom):
2266 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2267 break;
2268 case(Elements):
2269 case(ReducedElements):
2270 numTags=mesh->Elements->numTagsInUse;
2271 break;
2272 case(FaceElements):
2273 case(ReducedFaceElements):
2274 numTags=mesh->FaceElements->numTagsInUse;
2275 break;
2276 case(Points):
2277 numTags=mesh->Points->numTagsInUse;
2278 break;
2279 case(ContactElementsZero):
2280 case(ReducedContactElementsZero):
2281 case(ContactElementsOne):
2282 case(ReducedContactElementsOne):
2283 numTags=mesh->ContactElements->numTagsInUse;
2284 break;
2285 default:
2286 stringstream temp;
2287 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2288 throw FinleyAdapterException(temp.str());
2289 }
2290 return numTags;
2291 }
2292
2293 const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2294 {
2295 Finley_Mesh* mesh=m_finleyMesh.get();
2296 index_t* tags=NULL;
2297 switch(functionSpaceCode) {
2298 case(Nodes):
2299 tags=mesh->Nodes->tagsInUse;
2300 break;
2301 case(ReducedNodes):
2302 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2303 break;
2304 case(DegreesOfFreedom):
2305 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2306 break;
2307 case(ReducedDegreesOfFreedom):
2308 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2309 break;
2310 case(Elements):
2311 case(ReducedElements):
2312 tags=mesh->Elements->tagsInUse;
2313 break;
2314 case(FaceElements):
2315 case(ReducedFaceElements):
2316 tags=mesh->FaceElements->tagsInUse;
2317 break;
2318 case(Points):
2319 tags=mesh->Points->tagsInUse;
2320 break;
2321 case(ContactElementsZero):
2322 case(ReducedContactElementsZero):
2323 case(ContactElementsOne):
2324 case(ReducedContactElementsOne):
2325 tags=mesh->ContactElements->tagsInUse;
2326 break;
2327 default:
2328 stringstream temp;
2329 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2330 throw FinleyAdapterException(temp.str());
2331 }
2332 return tags;
2333 }
2334
2335
2336 bool MeshAdapter::canTag(int functionSpaceCode) const
2337 {
2338 switch(functionSpaceCode) {
2339 case(Nodes):
2340 case(Elements):
2341 case(ReducedElements):
2342 case(FaceElements):
2343 case(ReducedFaceElements):
2344 case(Points):
2345 case(ContactElementsZero):
2346 case(ReducedContactElementsZero):
2347 case(ContactElementsOne):
2348 case(ReducedContactElementsOne):
2349 return true;
2350 case(ReducedNodes):
2351 case(DegreesOfFreedom):
2352 case(ReducedDegreesOfFreedom):
2353 return false;
2354 default:
2355 return false;
2356 }
2357 }
2358
2359 AbstractDomain::StatusType MeshAdapter::getStatus() const
2360 {
2361 Finley_Mesh* mesh=m_finleyMesh.get();
2362 return Finley_Mesh_getStatus(mesh);
2363 }
2364
2365 int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2366 {
2367
2368 Finley_Mesh* mesh=m_finleyMesh.get();
2369 int order =-1;
2370 switch(functionSpaceCode) {
2371 case(Nodes):
2372 case(DegreesOfFreedom):
2373 order=mesh->approximationOrder;
2374 break;
2375 case(ReducedNodes):
2376 case(ReducedDegreesOfFreedom):
2377 order=mesh->reducedApproximationOrder;
2378 break;
2379 case(Elements):
2380 case(FaceElements):
2381 case(Points):
2382 case(ContactElementsZero):
2383 case(ContactElementsOne):
2384 order=mesh->integrationOrder;
2385 break;
2386 case(ReducedElements):
2387 case(ReducedFaceElements):
2388 case(ReducedContactElementsZero):
2389 case(ReducedContactElementsOne):
2390 order=mesh->reducedIntegrationOrder;
2391 break;
2392 default:
2393 stringstream temp;
2394 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2395 throw FinleyAdapterException(temp.str());
2396 }
2397 return order;
2398 }
2399
2400 ReferenceElementSetWrapper::ReferenceElementSetWrapper(ElementTypeId id, index_t order, index_t reducedOrder)
2401 {
2402 m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2403 }
2404
2405 ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2406 {
2407 Finley_ReferenceElementSet_dealloc(m_refSet);
2408 }
2409
2410 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26