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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26