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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26