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

Contents of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26