/[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 4114 - (show annotations)
Fri Dec 14 04:24:46 2012 UTC (6 years, 9 months ago) by caltinay
File size: 89142 byte(s)
Time to remove deprecated saveVTK/DX methods from Data and Domain.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26