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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4437 - (show annotations)
Tue Jun 4 05:37:18 2013 UTC (6 years ago) by caltinay
File size: 94260 byte(s)
finley:
-Removed obsolete reference files (OpenDX)
-Converted tag map to a std::map
-loosened tests to ignore tag order in file
-added some decorators for test skips

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26