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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26