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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4255 - (show annotations)
Wed Feb 27 03:06:21 2013 UTC (6 years, 1 month ago) by jfenwick
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 89528 byte(s)
Hopefully, this will address the interpolation problems.

New canInterpolate() function exposed to python which calls probeInterpolation

AbstractDomain now has an additional virtual.
preferredIntrpolationOnDomain()

This will return 0 if interpolation is impossible, 1 if possible and preferred.
It will return -1 if interpolation is possible and preferred in the 
oposite direction.

A value of -1 does not say that the proposed interpolation is possible or not.
Rather it indicates "please use the other way".
If you really _need_ to test it that way, use probeInterpolationOnDomain 


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 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 signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1985 {
1986 if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1987 {
1988 return 1;
1989 }
1990 else if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1991 {
1992 return -1;
1993 }
1994 return 0;
1995 }
1996
1997 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1998 {
1999 return false;
2000 }
2001
2002 bool MeshAdapter::operator==(const AbstractDomain& other) const
2003 {
2004 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2005 if (temp!=0) {
2006 return (m_finleyMesh==temp->m_finleyMesh);
2007 } else {
2008 return false;
2009 }
2010 }
2011
2012 bool MeshAdapter::operator!=(const AbstractDomain& other) const
2013 {
2014 return !(operator==(other));
2015 }
2016
2017 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2018 {
2019 Finley_Mesh* mesh=m_finleyMesh.get();
2020 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2021 package, symmetry, mesh->MPIInfo);
2022 }
2023
2024 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2025 {
2026 Finley_Mesh* mesh=m_finleyMesh.get();
2027 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2028 package, symmetry, mesh->MPIInfo);
2029 }
2030
2031 escript::Data MeshAdapter::getX() const
2032 {
2033 return continuousFunction(*this).getX();
2034 }
2035
2036 escript::Data MeshAdapter::getNormal() const
2037 {
2038 return functionOnBoundary(*this).getNormal();
2039 }
2040
2041 escript::Data MeshAdapter::getSize() const
2042 {
2043 return escript::function(*this).getSize();
2044 }
2045
2046 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2047 {
2048 int *out = NULL;
2049 Finley_Mesh* mesh=m_finleyMesh.get();
2050 switch (functionSpaceType) {
2051 case(Nodes):
2052 out=mesh->Nodes->Id;
2053 break;
2054 case(ReducedNodes):
2055 out=mesh->Nodes->reducedNodesId;
2056 break;
2057 case(Elements):
2058 out=mesh->Elements->Id;
2059 break;
2060 case(ReducedElements):
2061 out=mesh->Elements->Id;
2062 break;
2063 case(FaceElements):
2064 out=mesh->FaceElements->Id;
2065 break;
2066 case(ReducedFaceElements):
2067 out=mesh->FaceElements->Id;
2068 break;
2069 case(Points):
2070 out=mesh->Points->Id;
2071 break;
2072 case(ContactElementsZero):
2073 out=mesh->ContactElements->Id;
2074 break;
2075 case(ReducedContactElementsZero):
2076 out=mesh->ContactElements->Id;
2077 break;
2078 case(ContactElementsOne):
2079 out=mesh->ContactElements->Id;
2080 break;
2081 case(ReducedContactElementsOne):
2082 out=mesh->ContactElements->Id;
2083 break;
2084 case(DegreesOfFreedom):
2085 out=mesh->Nodes->degreesOfFreedomId;
2086 break;
2087 case(ReducedDegreesOfFreedom):
2088 out=mesh->Nodes->reducedDegreesOfFreedomId;
2089 break;
2090 default:
2091 stringstream temp;
2092 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2093 throw FinleyAdapterException(temp.str());
2094 break;
2095 }
2096 return out;
2097 }
2098 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2099 {
2100 int out=0;
2101 Finley_Mesh* mesh=m_finleyMesh.get();
2102 switch (functionSpaceType) {
2103 case(Nodes):
2104 out=mesh->Nodes->Tag[sampleNo];
2105 break;
2106 case(ReducedNodes):
2107 throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2108 break;
2109 case(Elements):
2110 out=mesh->Elements->Tag[sampleNo];
2111 break;
2112 case(ReducedElements):
2113 out=mesh->Elements->Tag[sampleNo];
2114 break;
2115 case(FaceElements):
2116 out=mesh->FaceElements->Tag[sampleNo];
2117 break;
2118 case(ReducedFaceElements):
2119 out=mesh->FaceElements->Tag[sampleNo];
2120 break;
2121 case(Points):
2122 out=mesh->Points->Tag[sampleNo];
2123 break;
2124 case(ContactElementsZero):
2125 out=mesh->ContactElements->Tag[sampleNo];
2126 break;
2127 case(ReducedContactElementsZero):
2128 out=mesh->ContactElements->Tag[sampleNo];
2129 break;
2130 case(ContactElementsOne):
2131 out=mesh->ContactElements->Tag[sampleNo];
2132 break;
2133 case(ReducedContactElementsOne):
2134 out=mesh->ContactElements->Tag[sampleNo];
2135 break;
2136 case(DegreesOfFreedom):
2137 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2138 break;
2139 case(ReducedDegreesOfFreedom):
2140 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2141 break;
2142 default:
2143 stringstream temp;
2144 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2145 throw FinleyAdapterException(temp.str());
2146 break;
2147 }
2148 return out;
2149 }
2150
2151
2152 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2153 {
2154 Finley_Mesh* mesh=m_finleyMesh.get();
2155 escriptDataC tmp=mask.getDataC();
2156 switch(functionSpaceType) {
2157 case(Nodes):
2158 Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2159 break;
2160 case(ReducedNodes):
2161 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2162 break;
2163 case(DegreesOfFreedom):
2164 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2165 break;
2166 case(ReducedDegreesOfFreedom):
2167 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2168 break;
2169 case(Elements):
2170 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2171 break;
2172 case(ReducedElements):
2173 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2174 break;
2175 case(FaceElements):
2176 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2177 break;
2178 case(ReducedFaceElements):
2179 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2180 break;
2181 case(Points):
2182 Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2183 break;
2184 case(ContactElementsZero):
2185 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2186 break;
2187 case(ReducedContactElementsZero):
2188 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2189 break;
2190 case(ContactElementsOne):
2191 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2192 break;
2193 case(ReducedContactElementsOne):
2194 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2195 break;
2196 default:
2197 stringstream temp;
2198 temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2199 throw FinleyAdapterException(temp.str());
2200 }
2201 checkFinleyError();
2202 return;
2203 }
2204
2205 void MeshAdapter::setTagMap(const string& name, int tag)
2206 {
2207 Finley_Mesh* mesh=m_finleyMesh.get();
2208 Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2209 checkFinleyError();
2210 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2211 }
2212
2213 int MeshAdapter::getTag(const string& name) const
2214 {
2215 Finley_Mesh* mesh=m_finleyMesh.get();
2216 int tag=0;
2217 tag=Finley_Mesh_getTag(mesh, name.c_str());
2218 checkFinleyError();
2219 // throwStandardException("MeshAdapter::getTag is not implemented.");
2220 return tag;
2221 }
2222
2223 bool MeshAdapter::isValidTagName(const string& name) const
2224 {
2225 Finley_Mesh* mesh=m_finleyMesh.get();
2226 return Finley_Mesh_isValidTagName(mesh,name.c_str());
2227 }
2228
2229 string MeshAdapter::showTagNames() const
2230 {
2231 stringstream temp;
2232 Finley_Mesh* mesh=m_finleyMesh.get();
2233 Finley_TagMap* tag_map=mesh->TagMap;
2234 while (tag_map) {
2235 temp << tag_map->name;
2236 tag_map=tag_map->next;
2237 if (tag_map) temp << ", ";
2238 }
2239 return temp.str();
2240 }
2241
2242 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2243 {
2244 Finley_Mesh* mesh=m_finleyMesh.get();
2245 dim_t numTags=0;
2246 switch(functionSpaceCode) {
2247 case(Nodes):
2248 numTags=mesh->Nodes->numTagsInUse;
2249 break;
2250 case(ReducedNodes):
2251 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2252 break;
2253 case(DegreesOfFreedom):
2254 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2255 break;
2256 case(ReducedDegreesOfFreedom):
2257 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2258 break;
2259 case(Elements):
2260 case(ReducedElements):
2261 numTags=mesh->Elements->numTagsInUse;
2262 break;
2263 case(FaceElements):
2264 case(ReducedFaceElements):
2265 numTags=mesh->FaceElements->numTagsInUse;
2266 break;
2267 case(Points):
2268 numTags=mesh->Points->numTagsInUse;
2269 break;
2270 case(ContactElementsZero):
2271 case(ReducedContactElementsZero):
2272 case(ContactElementsOne):
2273 case(ReducedContactElementsOne):
2274 numTags=mesh->ContactElements->numTagsInUse;
2275 break;
2276 default:
2277 stringstream temp;
2278 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2279 throw FinleyAdapterException(temp.str());
2280 }
2281 return numTags;
2282 }
2283
2284 const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2285 {
2286 Finley_Mesh* mesh=m_finleyMesh.get();
2287 index_t* tags=NULL;
2288 switch(functionSpaceCode) {
2289 case(Nodes):
2290 tags=mesh->Nodes->tagsInUse;
2291 break;
2292 case(ReducedNodes):
2293 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2294 break;
2295 case(DegreesOfFreedom):
2296 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2297 break;
2298 case(ReducedDegreesOfFreedom):
2299 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2300 break;
2301 case(Elements):
2302 case(ReducedElements):
2303 tags=mesh->Elements->tagsInUse;
2304 break;
2305 case(FaceElements):
2306 case(ReducedFaceElements):
2307 tags=mesh->FaceElements->tagsInUse;
2308 break;
2309 case(Points):
2310 tags=mesh->Points->tagsInUse;
2311 break;
2312 case(ContactElementsZero):
2313 case(ReducedContactElementsZero):
2314 case(ContactElementsOne):
2315 case(ReducedContactElementsOne):
2316 tags=mesh->ContactElements->tagsInUse;
2317 break;
2318 default:
2319 stringstream temp;
2320 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2321 throw FinleyAdapterException(temp.str());
2322 }
2323 return tags;
2324 }
2325
2326
2327 bool MeshAdapter::canTag(int functionSpaceCode) const
2328 {
2329 switch(functionSpaceCode) {
2330 case(Nodes):
2331 case(Elements):
2332 case(ReducedElements):
2333 case(FaceElements):
2334 case(ReducedFaceElements):
2335 case(Points):
2336 case(ContactElementsZero):
2337 case(ReducedContactElementsZero):
2338 case(ContactElementsOne):
2339 case(ReducedContactElementsOne):
2340 return true;
2341 case(ReducedNodes):
2342 case(DegreesOfFreedom):
2343 case(ReducedDegreesOfFreedom):
2344 return false;
2345 default:
2346 return false;
2347 }
2348 }
2349
2350 AbstractDomain::StatusType MeshAdapter::getStatus() const
2351 {
2352 Finley_Mesh* mesh=m_finleyMesh.get();
2353 return Finley_Mesh_getStatus(mesh);
2354 }
2355
2356 int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2357 {
2358
2359 Finley_Mesh* mesh=m_finleyMesh.get();
2360 int order =-1;
2361 switch(functionSpaceCode) {
2362 case(Nodes):
2363 case(DegreesOfFreedom):
2364 order=mesh->approximationOrder;
2365 break;
2366 case(ReducedNodes):
2367 case(ReducedDegreesOfFreedom):
2368 order=mesh->reducedApproximationOrder;
2369 break;
2370 case(Elements):
2371 case(FaceElements):
2372 case(Points):
2373 case(ContactElementsZero):
2374 case(ContactElementsOne):
2375 order=mesh->integrationOrder;
2376 break;
2377 case(ReducedElements):
2378 case(ReducedFaceElements):
2379 case(ReducedContactElementsZero):
2380 case(ReducedContactElementsOne):
2381 order=mesh->reducedIntegrationOrder;
2382 break;
2383 default:
2384 stringstream temp;
2385 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2386 throw FinleyAdapterException(temp.str());
2387 }
2388 return order;
2389 }
2390
2391 bool MeshAdapter::supportsContactElements() const
2392 {
2393 return true;
2394 }
2395
2396 ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2397 {
2398 m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2399 }
2400
2401 ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2402 {
2403 Finley_ReferenceElementSet_dealloc(m_refSet);
2404 }
2405
2406 // points will be flattened
2407 void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2408 {
2409 const int dim = getDim();
2410 int numPoints=points.size()/dim;
2411 int numTags=tags.size();
2412 Finley_Mesh* mesh=m_finleyMesh.get();
2413
2414 if ( points.size() % dim != 0 )
2415 {
2416 throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2417 }
2418
2419 if ( (numTags > 0) && ( numPoints != numTags ) )
2420 throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2421
2422 double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2423 int* tags_ptr= TMPMEMALLOC(numPoints, int);
2424
2425 for (int i=0;i<numPoints;++i) {
2426 points_ptr[ i * dim ] = points[i * dim ];
2427 if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2428 if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];
2429 tags_ptr[i]=tags[i];
2430 }
2431
2432 Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2433 checkFinleyError();
2434
2435 TMPMEMFREE(points_ptr);
2436 TMPMEMFREE(tags_ptr);
2437 }
2438
2439
2440 // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2441 // {
2442 // const int dim = getDim();
2443 // int numPoints=boost::python::extract<int>(points.attr("__len__")());
2444 // int numTags=boost::python::extract<int>(tags.attr("__len__")());
2445 // Finley_Mesh* mesh=m_finleyMesh.get();
2446 //
2447 // if ( (numTags > 0) && ( numPoints != numTags ) )
2448 // throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2449 //
2450 // double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2451 // int* tags_ptr= TMPMEMALLOC(numPoints, int);
2452 //
2453 // for (int i=0;i<numPoints;++i) {
2454 // int tag_id=-1;
2455 // int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2456 // if ( numComps != dim ) {
2457 // stringstream temp;
2458 // temp << "Error - illegal number of components " << numComps << " for point " << i;
2459 // throw FinleyAdapterException(temp.str());
2460 // }
2461 // points_ptr[ i * dim ] = boost::python::extract<double>(points[i][0]);
2462 // if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2463 // if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2464 //
2465 // if ( numTags > 0) {
2466 // boost::python::extract<string> ex_str(tags[i]);
2467 // if ( ex_str.check() ) {
2468 // tag_id=getTag( ex_str());
2469 // } else {
2470 // boost::python::extract<int> ex_int(tags[i]);
2471 // if ( ex_int.check() ) {
2472 // tag_id=ex_int();
2473 // } else {
2474 // stringstream temp;
2475 // temp << "Error - unable to extract tag for point " << i;
2476 // throw FinleyAdapterException(temp.str());
2477 // }
2478 // }
2479 // }
2480 // tags_ptr[i]=tag_id;
2481 // }
2482 //
2483 // Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2484 // checkPasoError();
2485 //
2486 // TMPMEMFREE(points_ptr);
2487 // TMPMEMFREE(tags_ptr);
2488 // }
2489
2490 /*
2491 void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2492 {
2493 boost::python::list points = boost::python::list();
2494 boost::python::list tags = boost::python::list();
2495 points.append(point);
2496 tags.append(tag);
2497 addDiracPoints(points, tags);
2498 }
2499 */
2500
2501 /*
2502 void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2503 {
2504 boost::python::list points = boost::python::list();
2505 boost::python::list tags = boost::python::list();
2506 points.append(point);
2507 tags.append(tag);
2508 addDiracPoints(points, tags);
2509 }
2510 */
2511 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26