/[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 4410 - (show annotations)
Tue May 14 11:29:23 2013 UTC (5 years, 11 months ago) by caltinay
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 95049 byte(s)
Oops, I saw this coming.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26