/[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 3993 - (show annotations)
Wed Sep 26 07:00:55 2012 UTC (7 years ago) by caltinay
File size: 91124 byte(s)
Addressing mantis 651. Test for finley will fail under MPI since there is
sth wrong with the way ReducedContinuousFunction samples are interrogated.

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 //
1526 // Helper for the save* methods. Extracts optional data variable names and the
1527 // corresponding pointers from python dictionary. Caller must free arrays.
1528 //
1529 void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1530 {
1531 numData = boost::python::extract<int>(arg.attr("__len__")());
1532 /* win32 refactor */
1533 names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1534 data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1535 dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1536
1537 boost::python::list keys=arg.keys();
1538 for (int i=0; i<numData; ++i) {
1539 string n=boost::python::extract<string>(keys[i]);
1540 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1541 if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1542 throw FinleyAdapterException("Error: Data must be defined on same Domain");
1543 data[i] = d.getDataC();
1544 dataPtr[i] = &(data[i]);
1545 names[i] = TMPMEMALLOC(n.length()+1, char);
1546 strcpy(names[i], n.c_str());
1547 }
1548 }
1549
1550 //
1551 // saves mesh and optionally data arrays in openDX format
1552 //
1553 void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1554 {
1555 int num_data;
1556 char **names;
1557 escriptDataC *data;
1558 escriptDataC **ptr_data;
1559
1560 extractArgsFromDict(arg, num_data, names, data, ptr_data);
1561 Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1562 checkFinleyError();
1563
1564 /* win32 refactor */
1565 TMPMEMFREE(data);
1566 TMPMEMFREE(ptr_data);
1567 for(int i=0; i<num_data; i++)
1568 {
1569 TMPMEMFREE(names[i]);
1570 }
1571 TMPMEMFREE(names);
1572 }
1573
1574 //
1575 // saves mesh and optionally data arrays in VTK format
1576 //
1577 void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg, const string& metadata, const string& metadata_schema) const
1578 {
1579 boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1580 pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1581 metadata, metadata_schema, arg);
1582 }
1583
1584 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1585 {
1586 if (getMPISize() > 1) {
1587 #ifdef ESYS_MPI
1588 index_t myFirstNode=0, myLastNode=0, k=0;
1589 index_t* globalNodeIndex=0;
1590 Finley_Mesh* mesh_p=m_finleyMesh.get();
1591 if (fs_code == FINLEY_REDUCED_NODES)
1592 {
1593 myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1594 myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1595 globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1596 }
1597 else if (fs_code == FINLEY_NODES)
1598 {
1599 myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1600 myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1601 globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1602 }
1603 else
1604 {
1605 throw FinleyAdapterException("Unsupported function space type for ownSample()");
1606 }
1607
1608 k=globalNodeIndex[id];
1609 return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1610 #endif
1611 }
1612 return true;
1613 }
1614
1615
1616 //
1617 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1618 //
1619 ASM_ptr MeshAdapter::newSystemMatrix(
1620 const int row_blocksize,
1621 const escript::FunctionSpace& row_functionspace,
1622 const int column_blocksize,
1623 const escript::FunctionSpace& column_functionspace,
1624 const int type) const
1625 {
1626 int reduceRowOrder=0;
1627 int reduceColOrder=0;
1628 // is the domain right?
1629 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1630 if (row_domain!=*this)
1631 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1632 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1633 if (col_domain!=*this)
1634 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1635 // is the function space type right
1636 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1637 reduceRowOrder=0;
1638 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1639 reduceRowOrder=1;
1640 } else {
1641 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1642 }
1643 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1644 reduceColOrder=0;
1645 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1646 reduceColOrder=1;
1647 } else {
1648 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1649 }
1650 // generate matrix:
1651
1652 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1653 checkFinleyError();
1654 Paso_SystemMatrix* fsystemMatrix;
1655 int trilinos = 0;
1656 if (trilinos) {
1657 #ifdef TRILINOS
1658 /* Allocation Epetra_VrbMatrix here */
1659 #endif
1660 }
1661 else {
1662 fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1663 }
1664 checkPasoError();
1665 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1666 SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1667 return ASM_ptr(sma);
1668 // return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1669 }
1670
1671 //
1672 // creates a TransportProblemAdapter
1673 //
1674 ATP_ptr MeshAdapter::newTransportProblem(
1675 const int blocksize,
1676 const escript::FunctionSpace& functionspace,
1677 const int type) const
1678 {
1679 int reduceOrder=0;
1680 // is the domain right?
1681 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1682 if (domain!=*this)
1683 throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1684 // is the function space type right
1685 if (functionspace.getTypeCode()==DegreesOfFreedom) {
1686 reduceOrder=0;
1687 } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1688 reduceOrder=1;
1689 } else {
1690 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1691 }
1692 // generate matrix:
1693
1694 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1695 checkFinleyError();
1696 Paso_TransportProblem* transportProblem;
1697 transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1698 checkPasoError();
1699 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1700 TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1701 return ATP_ptr(tpa);
1702 // return TransportProblemAdapter(transportProblem, blocksize,functionspace);
1703 }
1704
1705 //
1706 // vtkObject MeshAdapter::createVtkObject() const
1707 // TODO:
1708 //
1709
1710 //
1711 // returns true if data at the atom_type is considered as being cell centered:
1712 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1713 {
1714 switch(functionSpaceCode) {
1715 case(Nodes):
1716 case(DegreesOfFreedom):
1717 case(ReducedDegreesOfFreedom):
1718 return false;
1719 break;
1720 case(Elements):
1721 case(FaceElements):
1722 case(Points):
1723 case(ContactElementsZero):
1724 case(ContactElementsOne):
1725 case(ReducedElements):
1726 case(ReducedFaceElements):
1727 case(ReducedContactElementsZero):
1728 case(ReducedContactElementsOne):
1729 return true;
1730 break;
1731 default:
1732 stringstream temp;
1733 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1734 throw FinleyAdapterException(temp.str());
1735 break;
1736 }
1737 checkFinleyError();
1738 return false;
1739 }
1740
1741 bool
1742 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1743 {
1744 /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1745 class 1: DOF <-> Nodes
1746 class 2: ReducedDOF <-> ReducedNodes
1747 class 3: Points
1748 class 4: Elements
1749 class 5: ReducedElements
1750 class 6: FaceElements
1751 class 7: ReducedFaceElements
1752 class 8: ContactElementZero <-> ContactElementOne
1753 class 9: ReducedContactElementZero <-> ReducedContactElementOne
1754
1755 There is also a set of lines. Interpolation is possible down a line but not between lines.
1756 class 1 and 2 belong to all lines so aren't considered.
1757 line 0: class 3
1758 line 1: class 4,5
1759 line 2: class 6,7
1760 line 3: class 8,9
1761
1762 For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1763 eg hasnodes is true if we have at least one instance of Nodes.
1764 */
1765 if (fs.empty())
1766 {
1767 return false;
1768 }
1769 vector<int> hasclass(10);
1770 vector<int> hasline(4);
1771 bool hasnodes=false;
1772 bool hasrednodes=false;
1773 bool hascez=false;
1774 bool hasrcez=false;
1775 for (int i=0;i<fs.size();++i)
1776 {
1777 switch(fs[i])
1778 {
1779 case(Nodes): hasnodes=true; // no break is deliberate
1780 case(DegreesOfFreedom):
1781 hasclass[1]=1;
1782 break;
1783 case(ReducedNodes): hasrednodes=true; // no break is deliberate
1784 case(ReducedDegreesOfFreedom):
1785 hasclass[2]=1;
1786 break;
1787 case(Points):
1788 hasline[0]=1;
1789 hasclass[3]=1;
1790 break;
1791 case(Elements):
1792 hasclass[4]=1;
1793 hasline[1]=1;
1794 break;
1795 case(ReducedElements):
1796 hasclass[5]=1;
1797 hasline[1]=1;
1798 break;
1799 case(FaceElements):
1800 hasclass[6]=1;
1801 hasline[2]=1;
1802 break;
1803 case(ReducedFaceElements):
1804 hasclass[7]=1;
1805 hasline[2]=1;
1806 break;
1807 case(ContactElementsZero): hascez=true; // no break is deliberate
1808 case(ContactElementsOne):
1809 hasclass[8]=1;
1810 hasline[3]=1;
1811 break;
1812 case(ReducedContactElementsZero): hasrcez=true; // no break is deliberate
1813 case(ReducedContactElementsOne):
1814 hasclass[9]=1;
1815 hasline[3]=1;
1816 break;
1817 default:
1818 return false;
1819 }
1820 }
1821 int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1822 // fail if we have more than one leaf group
1823
1824 if (totlines>1)
1825 {
1826 return false; // there are at least two branches we can't interpolate between
1827 }
1828 else if (totlines==1)
1829 {
1830 if (hasline[0]==1) // we have points
1831 {
1832 resultcode=Points;
1833 }
1834 else if (hasline[1]==1)
1835 {
1836 if (hasclass[5]==1)
1837 {
1838 resultcode=ReducedElements;
1839 }
1840 else
1841 {
1842 resultcode=Elements;
1843 }
1844 }
1845 else if (hasline[2]==1)
1846 {
1847 if (hasclass[7]==1)
1848 {
1849 resultcode=ReducedFaceElements;
1850 }
1851 else
1852 {
1853 resultcode=FaceElements;
1854 }
1855 }
1856 else // so we must be in line3
1857 {
1858 if (hasclass[9]==1)
1859 {
1860 // need something from class 9
1861 resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1862 }
1863 else
1864 {
1865 // something from class 8
1866 resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1867 }
1868 }
1869 }
1870 else // totlines==0
1871 {
1872 if (hasclass[2]==1)
1873 {
1874 // something from class 2
1875 resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1876 }
1877 else
1878 { // something from class 1
1879 resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1880 }
1881 }
1882 return true;
1883 }
1884
1885 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1886 {
1887 switch(functionSpaceType_source) {
1888 case(Nodes):
1889 switch(functionSpaceType_target) {
1890 case(Nodes):
1891 case(ReducedNodes):
1892 case(ReducedDegreesOfFreedom):
1893 case(DegreesOfFreedom):
1894 case(Elements):
1895 case(ReducedElements):
1896 case(FaceElements):
1897 case(ReducedFaceElements):
1898 case(Points):
1899 case(ContactElementsZero):
1900 case(ReducedContactElementsZero):
1901 case(ContactElementsOne):
1902 case(ReducedContactElementsOne):
1903 return true;
1904 default:
1905 stringstream temp;
1906 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1907 throw FinleyAdapterException(temp.str());
1908 }
1909 break;
1910 case(ReducedNodes):
1911 switch(functionSpaceType_target) {
1912 case(ReducedNodes):
1913 case(ReducedDegreesOfFreedom):
1914 case(Elements):
1915 case(ReducedElements):
1916 case(FaceElements):
1917 case(ReducedFaceElements):
1918 case(Points):
1919 case(ContactElementsZero):
1920 case(ReducedContactElementsZero):
1921 case(ContactElementsOne):
1922 case(ReducedContactElementsOne):
1923 return true;
1924 case(Nodes):
1925 case(DegreesOfFreedom):
1926 return false;
1927 default:
1928 stringstream temp;
1929 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1930 throw FinleyAdapterException(temp.str());
1931 }
1932 break;
1933 case(Elements):
1934 if (functionSpaceType_target==Elements) {
1935 return true;
1936 } else if (functionSpaceType_target==ReducedElements) {
1937 return true;
1938 } else {
1939 return false;
1940 }
1941 case(ReducedElements):
1942 if (functionSpaceType_target==ReducedElements) {
1943 return true;
1944 } else {
1945 return false;
1946 }
1947 case(FaceElements):
1948 if (functionSpaceType_target==FaceElements) {
1949 return true;
1950 } else if (functionSpaceType_target==ReducedFaceElements) {
1951 return true;
1952 } else {
1953 return false;
1954 }
1955 case(ReducedFaceElements):
1956 if (functionSpaceType_target==ReducedFaceElements) {
1957 return true;
1958 } else {
1959 return false;
1960 }
1961 case(Points):
1962 if (functionSpaceType_target==Points) {
1963 return true;
1964 } else {
1965 return false;
1966 }
1967 case(ContactElementsZero):
1968 case(ContactElementsOne):
1969 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1970 return true;
1971 } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1972 return true;
1973 } else {
1974 return false;
1975 }
1976 case(ReducedContactElementsZero):
1977 case(ReducedContactElementsOne):
1978 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1979 return true;
1980 } else {
1981 return false;
1982 }
1983 case(DegreesOfFreedom):
1984 switch(functionSpaceType_target) {
1985 case(ReducedDegreesOfFreedom):
1986 case(DegreesOfFreedom):
1987 case(Nodes):
1988 case(ReducedNodes):
1989 case(Elements):
1990 case(ReducedElements):
1991 case(Points):
1992 case(FaceElements):
1993 case(ReducedFaceElements):
1994 case(ContactElementsZero):
1995 case(ReducedContactElementsZero):
1996 case(ContactElementsOne):
1997 case(ReducedContactElementsOne):
1998 return true;
1999 default:
2000 stringstream temp;
2001 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
2002 throw FinleyAdapterException(temp.str());
2003 }
2004 break;
2005 case(ReducedDegreesOfFreedom):
2006 switch(functionSpaceType_target) {
2007 case(ReducedDegreesOfFreedom):
2008 case(ReducedNodes):
2009 case(Elements):
2010 case(ReducedElements):
2011 case(FaceElements):
2012 case(ReducedFaceElements):
2013 case(Points):
2014 case(ContactElementsZero):
2015 case(ReducedContactElementsZero):
2016 case(ContactElementsOne):
2017 case(ReducedContactElementsOne):
2018 return true;
2019 case(Nodes):
2020 case(DegreesOfFreedom):
2021 return false;
2022 default:
2023 stringstream temp;
2024 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
2025 throw FinleyAdapterException(temp.str());
2026 }
2027 break;
2028 default:
2029 stringstream temp;
2030 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
2031 throw FinleyAdapterException(temp.str());
2032 break;
2033 }
2034 checkFinleyError();
2035 return false;
2036 }
2037
2038 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
2039 {
2040 return false;
2041 }
2042
2043 bool MeshAdapter::operator==(const AbstractDomain& other) const
2044 {
2045 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2046 if (temp!=0) {
2047 return (m_finleyMesh==temp->m_finleyMesh);
2048 } else {
2049 return false;
2050 }
2051 }
2052
2053 bool MeshAdapter::operator!=(const AbstractDomain& other) const
2054 {
2055 return !(operator==(other));
2056 }
2057
2058 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2059 {
2060 Finley_Mesh* mesh=m_finleyMesh.get();
2061 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2062 package, symmetry, mesh->MPIInfo);
2063 }
2064
2065 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2066 {
2067 Finley_Mesh* mesh=m_finleyMesh.get();
2068 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2069 package, symmetry, mesh->MPIInfo);
2070 }
2071
2072 escript::Data MeshAdapter::getX() const
2073 {
2074 return continuousFunction(*this).getX();
2075 }
2076
2077 escript::Data MeshAdapter::getNormal() const
2078 {
2079 return functionOnBoundary(*this).getNormal();
2080 }
2081
2082 escript::Data MeshAdapter::getSize() const
2083 {
2084 return escript::function(*this).getSize();
2085 }
2086
2087 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2088 {
2089 int *out = NULL;
2090 Finley_Mesh* mesh=m_finleyMesh.get();
2091 switch (functionSpaceType) {
2092 case(Nodes):
2093 out=mesh->Nodes->Id;
2094 break;
2095 case(ReducedNodes):
2096 out=mesh->Nodes->reducedNodesId;
2097 break;
2098 case(Elements):
2099 out=mesh->Elements->Id;
2100 break;
2101 case(ReducedElements):
2102 out=mesh->Elements->Id;
2103 break;
2104 case(FaceElements):
2105 out=mesh->FaceElements->Id;
2106 break;
2107 case(ReducedFaceElements):
2108 out=mesh->FaceElements->Id;
2109 break;
2110 case(Points):
2111 out=mesh->Points->Id;
2112 break;
2113 case(ContactElementsZero):
2114 out=mesh->ContactElements->Id;
2115 break;
2116 case(ReducedContactElementsZero):
2117 out=mesh->ContactElements->Id;
2118 break;
2119 case(ContactElementsOne):
2120 out=mesh->ContactElements->Id;
2121 break;
2122 case(ReducedContactElementsOne):
2123 out=mesh->ContactElements->Id;
2124 break;
2125 case(DegreesOfFreedom):
2126 out=mesh->Nodes->degreesOfFreedomId;
2127 break;
2128 case(ReducedDegreesOfFreedom):
2129 out=mesh->Nodes->reducedDegreesOfFreedomId;
2130 break;
2131 default:
2132 stringstream temp;
2133 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2134 throw FinleyAdapterException(temp.str());
2135 break;
2136 }
2137 return out;
2138 }
2139 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2140 {
2141 int out=0;
2142 Finley_Mesh* mesh=m_finleyMesh.get();
2143 switch (functionSpaceType) {
2144 case(Nodes):
2145 out=mesh->Nodes->Tag[sampleNo];
2146 break;
2147 case(ReducedNodes):
2148 throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2149 break;
2150 case(Elements):
2151 out=mesh->Elements->Tag[sampleNo];
2152 break;
2153 case(ReducedElements):
2154 out=mesh->Elements->Tag[sampleNo];
2155 break;
2156 case(FaceElements):
2157 out=mesh->FaceElements->Tag[sampleNo];
2158 break;
2159 case(ReducedFaceElements):
2160 out=mesh->FaceElements->Tag[sampleNo];
2161 break;
2162 case(Points):
2163 out=mesh->Points->Tag[sampleNo];
2164 break;
2165 case(ContactElementsZero):
2166 out=mesh->ContactElements->Tag[sampleNo];
2167 break;
2168 case(ReducedContactElementsZero):
2169 out=mesh->ContactElements->Tag[sampleNo];
2170 break;
2171 case(ContactElementsOne):
2172 out=mesh->ContactElements->Tag[sampleNo];
2173 break;
2174 case(ReducedContactElementsOne):
2175 out=mesh->ContactElements->Tag[sampleNo];
2176 break;
2177 case(DegreesOfFreedom):
2178 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2179 break;
2180 case(ReducedDegreesOfFreedom):
2181 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2182 break;
2183 default:
2184 stringstream temp;
2185 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2186 throw FinleyAdapterException(temp.str());
2187 break;
2188 }
2189 return out;
2190 }
2191
2192
2193 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2194 {
2195 Finley_Mesh* mesh=m_finleyMesh.get();
2196 escriptDataC tmp=mask.getDataC();
2197 switch(functionSpaceType) {
2198 case(Nodes):
2199 Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2200 break;
2201 case(ReducedNodes):
2202 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2203 break;
2204 case(DegreesOfFreedom):
2205 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2206 break;
2207 case(ReducedDegreesOfFreedom):
2208 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2209 break;
2210 case(Elements):
2211 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2212 break;
2213 case(ReducedElements):
2214 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2215 break;
2216 case(FaceElements):
2217 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2218 break;
2219 case(ReducedFaceElements):
2220 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2221 break;
2222 case(Points):
2223 Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2224 break;
2225 case(ContactElementsZero):
2226 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2227 break;
2228 case(ReducedContactElementsZero):
2229 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2230 break;
2231 case(ContactElementsOne):
2232 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2233 break;
2234 case(ReducedContactElementsOne):
2235 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2236 break;
2237 default:
2238 stringstream temp;
2239 temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2240 throw FinleyAdapterException(temp.str());
2241 }
2242 checkFinleyError();
2243 return;
2244 }
2245
2246 void MeshAdapter::setTagMap(const string& name, int tag)
2247 {
2248 Finley_Mesh* mesh=m_finleyMesh.get();
2249 Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2250 checkFinleyError();
2251 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2252 }
2253
2254 int MeshAdapter::getTag(const string& name) const
2255 {
2256 Finley_Mesh* mesh=m_finleyMesh.get();
2257 int tag=0;
2258 tag=Finley_Mesh_getTag(mesh, name.c_str());
2259 checkFinleyError();
2260 // throwStandardException("MeshAdapter::getTag is not implemented.");
2261 return tag;
2262 }
2263
2264 bool MeshAdapter::isValidTagName(const string& name) const
2265 {
2266 Finley_Mesh* mesh=m_finleyMesh.get();
2267 return Finley_Mesh_isValidTagName(mesh,name.c_str());
2268 }
2269
2270 string MeshAdapter::showTagNames() const
2271 {
2272 stringstream temp;
2273 Finley_Mesh* mesh=m_finleyMesh.get();
2274 Finley_TagMap* tag_map=mesh->TagMap;
2275 while (tag_map) {
2276 temp << tag_map->name;
2277 tag_map=tag_map->next;
2278 if (tag_map) temp << ", ";
2279 }
2280 return temp.str();
2281 }
2282
2283 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2284 {
2285 Finley_Mesh* mesh=m_finleyMesh.get();
2286 dim_t numTags=0;
2287 switch(functionSpaceCode) {
2288 case(Nodes):
2289 numTags=mesh->Nodes->numTagsInUse;
2290 break;
2291 case(ReducedNodes):
2292 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2293 break;
2294 case(DegreesOfFreedom):
2295 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2296 break;
2297 case(ReducedDegreesOfFreedom):
2298 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2299 break;
2300 case(Elements):
2301 case(ReducedElements):
2302 numTags=mesh->Elements->numTagsInUse;
2303 break;
2304 case(FaceElements):
2305 case(ReducedFaceElements):
2306 numTags=mesh->FaceElements->numTagsInUse;
2307 break;
2308 case(Points):
2309 numTags=mesh->Points->numTagsInUse;
2310 break;
2311 case(ContactElementsZero):
2312 case(ReducedContactElementsZero):
2313 case(ContactElementsOne):
2314 case(ReducedContactElementsOne):
2315 numTags=mesh->ContactElements->numTagsInUse;
2316 break;
2317 default:
2318 stringstream temp;
2319 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2320 throw FinleyAdapterException(temp.str());
2321 }
2322 return numTags;
2323 }
2324
2325 const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2326 {
2327 Finley_Mesh* mesh=m_finleyMesh.get();
2328 index_t* tags=NULL;
2329 switch(functionSpaceCode) {
2330 case(Nodes):
2331 tags=mesh->Nodes->tagsInUse;
2332 break;
2333 case(ReducedNodes):
2334 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2335 break;
2336 case(DegreesOfFreedom):
2337 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2338 break;
2339 case(ReducedDegreesOfFreedom):
2340 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2341 break;
2342 case(Elements):
2343 case(ReducedElements):
2344 tags=mesh->Elements->tagsInUse;
2345 break;
2346 case(FaceElements):
2347 case(ReducedFaceElements):
2348 tags=mesh->FaceElements->tagsInUse;
2349 break;
2350 case(Points):
2351 tags=mesh->Points->tagsInUse;
2352 break;
2353 case(ContactElementsZero):
2354 case(ReducedContactElementsZero):
2355 case(ContactElementsOne):
2356 case(ReducedContactElementsOne):
2357 tags=mesh->ContactElements->tagsInUse;
2358 break;
2359 default:
2360 stringstream temp;
2361 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2362 throw FinleyAdapterException(temp.str());
2363 }
2364 return tags;
2365 }
2366
2367
2368 bool MeshAdapter::canTag(int functionSpaceCode) const
2369 {
2370 switch(functionSpaceCode) {
2371 case(Nodes):
2372 case(Elements):
2373 case(ReducedElements):
2374 case(FaceElements):
2375 case(ReducedFaceElements):
2376 case(Points):
2377 case(ContactElementsZero):
2378 case(ReducedContactElementsZero):
2379 case(ContactElementsOne):
2380 case(ReducedContactElementsOne):
2381 return true;
2382 case(ReducedNodes):
2383 case(DegreesOfFreedom):
2384 case(ReducedDegreesOfFreedom):
2385 return false;
2386 default:
2387 return false;
2388 }
2389 }
2390
2391 AbstractDomain::StatusType MeshAdapter::getStatus() const
2392 {
2393 Finley_Mesh* mesh=m_finleyMesh.get();
2394 return Finley_Mesh_getStatus(mesh);
2395 }
2396
2397 int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2398 {
2399
2400 Finley_Mesh* mesh=m_finleyMesh.get();
2401 int order =-1;
2402 switch(functionSpaceCode) {
2403 case(Nodes):
2404 case(DegreesOfFreedom):
2405 order=mesh->approximationOrder;
2406 break;
2407 case(ReducedNodes):
2408 case(ReducedDegreesOfFreedom):
2409 order=mesh->reducedApproximationOrder;
2410 break;
2411 case(Elements):
2412 case(FaceElements):
2413 case(Points):
2414 case(ContactElementsZero):
2415 case(ContactElementsOne):
2416 order=mesh->integrationOrder;
2417 break;
2418 case(ReducedElements):
2419 case(ReducedFaceElements):
2420 case(ReducedContactElementsZero):
2421 case(ReducedContactElementsOne):
2422 order=mesh->reducedIntegrationOrder;
2423 break;
2424 default:
2425 stringstream temp;
2426 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2427 throw FinleyAdapterException(temp.str());
2428 }
2429 return order;
2430 }
2431
2432 bool MeshAdapter::supportsContactElements() const
2433 {
2434 return true;
2435 }
2436
2437 ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2438 {
2439 m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2440 }
2441
2442 ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2443 {
2444 Finley_ReferenceElementSet_dealloc(m_refSet);
2445 }
2446
2447 // points will be flattened
2448 void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2449 {
2450 const int dim = getDim();
2451 int numPoints=points.size()/dim;
2452 int numTags=tags.size();
2453 Finley_Mesh* mesh=m_finleyMesh.get();
2454
2455 if ( points.size() % dim != 0 )
2456 {
2457 throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2458 }
2459
2460 if ( (numTags > 0) && ( numPoints != numTags ) )
2461 throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2462
2463 double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2464 int* tags_ptr= TMPMEMALLOC(numPoints, int);
2465
2466 for (int i=0;i<numPoints;++i) {
2467 points_ptr[ i * dim ] = points[i * dim ];
2468 if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2469 if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];
2470 tags_ptr[i]=tags[i];
2471 }
2472
2473 Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2474 checkFinleyError();
2475
2476 TMPMEMFREE(points_ptr);
2477 TMPMEMFREE(tags_ptr);
2478 }
2479
2480
2481 // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2482 // {
2483 // const int dim = getDim();
2484 // int numPoints=boost::python::extract<int>(points.attr("__len__")());
2485 // int numTags=boost::python::extract<int>(tags.attr("__len__")());
2486 // Finley_Mesh* mesh=m_finleyMesh.get();
2487 //
2488 // if ( (numTags > 0) && ( numPoints != numTags ) )
2489 // throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2490 //
2491 // double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2492 // int* tags_ptr= TMPMEMALLOC(numPoints, int);
2493 //
2494 // for (int i=0;i<numPoints;++i) {
2495 // int tag_id=-1;
2496 // int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2497 // if ( numComps != dim ) {
2498 // stringstream temp;
2499 // temp << "Error - illegal number of components " << numComps << " for point " << i;
2500 // throw FinleyAdapterException(temp.str());
2501 // }
2502 // points_ptr[ i * dim ] = boost::python::extract<double>(points[i][0]);
2503 // if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2504 // if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2505 //
2506 // if ( numTags > 0) {
2507 // boost::python::extract<string> ex_str(tags[i]);
2508 // if ( ex_str.check() ) {
2509 // tag_id=getTag( ex_str());
2510 // } else {
2511 // boost::python::extract<int> ex_int(tags[i]);
2512 // if ( ex_int.check() ) {
2513 // tag_id=ex_int();
2514 // } else {
2515 // stringstream temp;
2516 // temp << "Error - unable to extract tag for point " << i;
2517 // throw FinleyAdapterException(temp.str());
2518 // }
2519 // }
2520 // }
2521 // tags_ptr[i]=tag_id;
2522 // }
2523 //
2524 // Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2525 // checkPasoError();
2526 //
2527 // TMPMEMFREE(points_ptr);
2528 // TMPMEMFREE(tags_ptr);
2529 // }
2530
2531 /*
2532 void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2533 {
2534 boost::python::list points = boost::python::list();
2535 boost::python::list tags = boost::python::list();
2536 points.append(point);
2537 tags.append(tag);
2538 addDiracPoints(points, tags);
2539 }
2540 */
2541
2542 /*
2543 void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2544 {
2545 boost::python::list points = boost::python::list();
2546 boost::python::list tags = boost::python::list();
2547 points.append(point);
2548 tags.append(tag);
2549 addDiracPoints(points, tags);
2550 }
2551 */
2552 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26