/[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 3998 - (show annotations)
Thu Sep 27 01:17:28 2012 UTC (6 years, 11 months ago) by caltinay
File size: 91311 byte(s)
Finley now throws in ownSample() for reduced function with mpisize>1.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26