/[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 3675 - (show annotations)
Thu Nov 17 00:53:38 2011 UTC (7 years, 5 months ago) by jfenwick
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 90988 byte(s)
pasowrap joins the trunk.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26