/[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 4428 - (show annotations)
Thu May 30 06:39:10 2013 UTC (5 years, 10 months ago) by caltinay
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 94529 byte(s)
finley's NodeFile is now a class.
Associated changes:
- use of some C++ constructs/functions/types (1st pass)
- removal of obsolete pointer check
- merging of some duplicated code
- ...

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 m_finleyMesh.get()->Nodes->getGlobalNumNodes();
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=mesh->Nodes->getNumNodes();
668 break;
669 case ReducedNodes:
670 numDataPointsPerSample=1;
671 numSamples=mesh->Nodes->getNumReducedNodes();
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=mesh->Nodes->getNumDegreesOfFreedom();
731 }
732 break;
733 case ReducedDegreesOfFreedom:
734 if (mesh->Nodes!=NULL) {
735 numDataPointsPerSample=1;
736 numSamples=mesh->Nodes->getNumReducedDegreesOfFreedom();
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 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1449 if (newDomain!=*this)
1450 throw FinleyAdapterException("Error - Illegal domain of new point locations");
1451 if (new_x.getFunctionSpace() == continuousFunction(*this)) {
1452 Finley_Mesh_setCoordinates(mesh, new_x);
1453 } else {
1454 throw FinleyAdapterException("As of escript version 3.3 SetX() only accepts ContinuousFunction arguments. Please interpolate.");
1455 //escript::Data new_x_inter=escript::Data(new_x, continuousFunction(*this));
1456 //Finley_Mesh_setCoordinates(mesh, new_x_inter);
1457 }
1458 checkFinleyError();
1459 }
1460
1461 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1462 {
1463 if (getMPISize() > 1) {
1464 #ifdef ESYS_MPI
1465 index_t myFirstNode=0, myLastNode=0, k=0;
1466 index_t* globalNodeIndex=0;
1467 Finley_Mesh* mesh_p=m_finleyMesh.get();
1468 /*
1469 * this method is only used by saveDataCSV which would use the returned
1470 * values for reduced nodes wrongly so this case is disabled for now
1471 if (fs_code == FINLEY_REDUCED_NODES) {
1472 myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1473 myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1474 globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1475 } else
1476 */
1477 if (fs_code == FINLEY_NODES) {
1478 myFirstNode = mesh_p->Nodes->getFirstNode();
1479 myLastNode = mesh_p->Nodes->getLastNode();
1480 globalNodeIndex = mesh_p->Nodes->borrowGlobalNodesIndex();
1481 } else {
1482 throw FinleyAdapterException("Unsupported function space type for ownSample()");
1483 }
1484
1485 k=globalNodeIndex[id];
1486 return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1487 #endif
1488 }
1489 return true;
1490 }
1491
1492
1493 //
1494 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1495 //
1496 escript::ASM_ptr MeshAdapter::newSystemMatrix(const int row_blocksize,
1497 const escript::FunctionSpace& row_functionspace,
1498 const int column_blocksize,
1499 const escript::FunctionSpace& column_functionspace,
1500 const int type) const
1501 {
1502 // is the domain right?
1503 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1504 if (row_domain!=*this)
1505 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1506 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1507 if (col_domain!=*this)
1508 throw FinleyAdapterException("Error - domain of column function space does not match the domain of matrix generator.");
1509
1510 int reduceRowOrder=0;
1511 int reduceColOrder=0;
1512 // is the function space type right?
1513 if (row_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1514 reduceRowOrder=1;
1515 } else if (row_functionspace.getTypeCode() != DegreesOfFreedom) {
1516 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1517 }
1518 if (column_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1519 reduceColOrder=1;
1520 } else if (column_functionspace.getTypeCode() != DegreesOfFreedom) {
1521 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1522 }
1523
1524 // generate matrix:
1525 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(
1526 getFinley_Mesh(), reduceRowOrder, reduceColOrder);
1527 checkFinleyError();
1528 Paso_SystemMatrix* fsystemMatrix;
1529 const int trilinos = 0;
1530 if (trilinos) {
1531 #ifdef TRILINOS
1532 // FIXME: Allocation Epetra_VrbMatrix here...
1533 #endif
1534 } else {
1535 fsystemMatrix=Paso_SystemMatrix_alloc(type, fsystemMatrixPattern,
1536 row_blocksize, column_blocksize, FALSE);
1537 }
1538 checkPasoError();
1539 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1540 SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1541 return escript::ASM_ptr(sma);
1542 }
1543
1544 //
1545 // creates a TransportProblemAdapter
1546 //
1547 escript::ATP_ptr MeshAdapter::newTransportProblem(const int blocksize,
1548 const escript::FunctionSpace& functionspace, const int type) const
1549 {
1550 // is the domain right?
1551 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1552 if (domain!=*this)
1553 throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1554
1555 // is the function space type right?
1556 int reduceOrder=0;
1557 if (functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1558 reduceOrder=1;
1559 } else if (functionspace.getTypeCode() != DegreesOfFreedom) {
1560 throw FinleyAdapterException("Error - illegal function space type for transport problem.");
1561 }
1562
1563 // generate transport problem:
1564 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(
1565 getFinley_Mesh(), reduceOrder, reduceOrder);
1566 checkFinleyError();
1567 Paso_TransportProblem* transportProblem;
1568 transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern, blocksize);
1569 checkPasoError();
1570 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1571 TransportProblemAdapter* tpa=new TransportProblemAdapter(
1572 transportProblem, blocksize, functionspace);
1573 return escript::ATP_ptr(tpa);
1574 }
1575
1576 //
1577 // returns true if data on functionSpaceCode is considered as being cell centered
1578 //
1579 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1580 {
1581 switch(functionSpaceCode) {
1582 case Nodes:
1583 case DegreesOfFreedom:
1584 case ReducedDegreesOfFreedom:
1585 return false;
1586 case Elements:
1587 case FaceElements:
1588 case Points:
1589 case ContactElementsZero:
1590 case ContactElementsOne:
1591 case ReducedElements:
1592 case ReducedFaceElements:
1593 case ReducedContactElementsZero:
1594 case ReducedContactElementsOne:
1595 return true;
1596 default:
1597 stringstream temp;
1598 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1599 throw FinleyAdapterException(temp.str());
1600 break;
1601 }
1602 return false;
1603 }
1604
1605 bool
1606 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1607 {
1608 /* The idea is to use equivalence classes. [Types which can be interpolated
1609 back and forth]:
1610 class 1: DOF <-> Nodes
1611 class 2: ReducedDOF <-> ReducedNodes
1612 class 3: Points
1613 class 4: Elements
1614 class 5: ReducedElements
1615 class 6: FaceElements
1616 class 7: ReducedFaceElements
1617 class 8: ContactElementZero <-> ContactElementOne
1618 class 9: ReducedContactElementZero <-> ReducedContactElementOne
1619
1620 There is also a set of lines. Interpolation is possible down a line but
1621 not between lines.
1622 class 1 and 2 belong to all lines so aren't considered.
1623 line 0: class 3
1624 line 1: class 4,5
1625 line 2: class 6,7
1626 line 3: class 8,9
1627
1628 For classes with multiple members (eg class 2) we have vars to record if
1629 there is at least one instance.
1630 e.g. hasnodes is true if we have at least one instance of Nodes.
1631 */
1632 if (fs.empty())
1633 return false;
1634
1635 vector<int> hasclass(10);
1636 vector<int> hasline(4);
1637 bool hasnodes=false;
1638 bool hasrednodes=false;
1639 bool hascez=false;
1640 bool hasrcez=false;
1641 for (int i=0;i<fs.size();++i) {
1642 switch(fs[i]) {
1643 case Nodes:
1644 hasnodes=true; // fall through
1645 case DegreesOfFreedom:
1646 hasclass[1]=1;
1647 break;
1648 case ReducedNodes:
1649 hasrednodes=true; // fall through
1650 case ReducedDegreesOfFreedom:
1651 hasclass[2]=1;
1652 break;
1653 case Points:
1654 hasline[0]=1;
1655 hasclass[3]=1;
1656 break;
1657 case Elements:
1658 hasclass[4]=1;
1659 hasline[1]=1;
1660 break;
1661 case ReducedElements:
1662 hasclass[5]=1;
1663 hasline[1]=1;
1664 break;
1665 case FaceElements:
1666 hasclass[6]=1;
1667 hasline[2]=1;
1668 break;
1669 case ReducedFaceElements:
1670 hasclass[7]=1;
1671 hasline[2]=1;
1672 break;
1673 case ContactElementsZero:
1674 hascez=true; // fall through
1675 case ContactElementsOne:
1676 hasclass[8]=1;
1677 hasline[3]=1;
1678 break;
1679 case ReducedContactElementsZero:
1680 hasrcez=true; // fall through
1681 case ReducedContactElementsOne:
1682 hasclass[9]=1;
1683 hasline[3]=1;
1684 break;
1685 default:
1686 return false;
1687 }
1688 }
1689 int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1690
1691 // fail if we have more than one leaf group
1692 if (totlines>1)
1693 return false; // there are at least two branches we can't interpolate between
1694 else if (totlines==1) {
1695 if (hasline[0]==1) // we have points
1696 resultcode=Points;
1697 else if (hasline[1]==1) {
1698 if (hasclass[5]==1)
1699 resultcode=ReducedElements;
1700 else
1701 resultcode=Elements;
1702 } else if (hasline[2]==1) {
1703 if (hasclass[7]==1)
1704 resultcode=ReducedFaceElements;
1705 else
1706 resultcode=FaceElements;
1707 } else { // so we must be in line3
1708 if (hasclass[9]==1) {
1709 // need something from class 9
1710 resultcode=(hasrcez ? ReducedContactElementsZero : ReducedContactElementsOne);
1711 } else {
1712 // something from class 8
1713 resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1714 }
1715 }
1716 } else { // totlines==0
1717 if (hasclass[2]==1) {
1718 // something from class 2
1719 resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
1720 } else {
1721 // something from class 1
1722 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
1723 }
1724 }
1725 return true;
1726 }
1727
1728 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,
1729 int functionSpaceType_target) const
1730 {
1731 switch(functionSpaceType_source) {
1732 case Nodes:
1733 switch(functionSpaceType_target) {
1734 case Nodes:
1735 case ReducedNodes:
1736 case ReducedDegreesOfFreedom:
1737 case DegreesOfFreedom:
1738 case Elements:
1739 case ReducedElements:
1740 case FaceElements:
1741 case ReducedFaceElements:
1742 case Points:
1743 case ContactElementsZero:
1744 case ReducedContactElementsZero:
1745 case ContactElementsOne:
1746 case ReducedContactElementsOne:
1747 return true;
1748 default:
1749 stringstream temp;
1750 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1751 throw FinleyAdapterException(temp.str());
1752 }
1753 break;
1754 case ReducedNodes:
1755 switch(functionSpaceType_target) {
1756 case ReducedNodes:
1757 case ReducedDegreesOfFreedom:
1758 case Elements:
1759 case ReducedElements:
1760 case FaceElements:
1761 case ReducedFaceElements:
1762 case Points:
1763 case ContactElementsZero:
1764 case ReducedContactElementsZero:
1765 case ContactElementsOne:
1766 case ReducedContactElementsOne:
1767 return true;
1768 case Nodes:
1769 case DegreesOfFreedom:
1770 return false;
1771 default:
1772 stringstream temp;
1773 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1774 throw FinleyAdapterException(temp.str());
1775 }
1776 break;
1777 case Elements:
1778 if (functionSpaceType_target==Elements) {
1779 return true;
1780 } else if (functionSpaceType_target==ReducedElements) {
1781 return true;
1782 } else {
1783 return false;
1784 }
1785 case ReducedElements:
1786 if (functionSpaceType_target==ReducedElements) {
1787 return true;
1788 } else {
1789 return false;
1790 }
1791 case FaceElements:
1792 if (functionSpaceType_target==FaceElements) {
1793 return true;
1794 } else if (functionSpaceType_target==ReducedFaceElements) {
1795 return true;
1796 } else {
1797 return false;
1798 }
1799 case ReducedFaceElements:
1800 if (functionSpaceType_target==ReducedFaceElements) {
1801 return true;
1802 } else {
1803 return false;
1804 }
1805 case Points:
1806 if (functionSpaceType_target==Points) {
1807 return true;
1808 } else {
1809 return false;
1810 }
1811 case ContactElementsZero:
1812 case ContactElementsOne:
1813 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1814 return true;
1815 } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1816 return true;
1817 } else {
1818 return false;
1819 }
1820 case ReducedContactElementsZero:
1821 case ReducedContactElementsOne:
1822 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1823 return true;
1824 } else {
1825 return false;
1826 }
1827 case DegreesOfFreedom:
1828 switch(functionSpaceType_target) {
1829 case ReducedDegreesOfFreedom:
1830 case DegreesOfFreedom:
1831 case Nodes:
1832 case ReducedNodes:
1833 case Elements:
1834 case ReducedElements:
1835 case Points:
1836 case FaceElements:
1837 case ReducedFaceElements:
1838 case ContactElementsZero:
1839 case ReducedContactElementsZero:
1840 case ContactElementsOne:
1841 case ReducedContactElementsOne:
1842 return true;
1843 default:
1844 stringstream temp;
1845 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1846 throw FinleyAdapterException(temp.str());
1847 }
1848 break;
1849 case ReducedDegreesOfFreedom:
1850 switch(functionSpaceType_target) {
1851 case ReducedDegreesOfFreedom:
1852 case ReducedNodes:
1853 case Elements:
1854 case ReducedElements:
1855 case FaceElements:
1856 case ReducedFaceElements:
1857 case Points:
1858 case ContactElementsZero:
1859 case ReducedContactElementsZero:
1860 case ContactElementsOne:
1861 case ReducedContactElementsOne:
1862 return true;
1863 case Nodes:
1864 case DegreesOfFreedom:
1865 return false;
1866 default:
1867 stringstream temp;
1868 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1869 throw FinleyAdapterException(temp.str());
1870 }
1871 break;
1872 default:
1873 stringstream temp;
1874 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1875 throw FinleyAdapterException(temp.str());
1876 break;
1877 }
1878 return false;
1879 }
1880
1881 signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source, int functionSpaceType_target) const
1882 {
1883 if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1884 return 1;
1885
1886 if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1887 return -1;
1888
1889 return 0;
1890 }
1891
1892 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,
1893 const escript::AbstractDomain& targetDomain,
1894 int functionSpaceType_target) const
1895 {
1896 return false;
1897 }
1898
1899 bool MeshAdapter::operator==(const escript::AbstractDomain& other) const
1900 {
1901 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1902 if (temp) {
1903 return (m_finleyMesh==temp->m_finleyMesh);
1904 } else {
1905 return false;
1906 }
1907 }
1908
1909 bool MeshAdapter::operator!=(const escript::AbstractDomain& other) const
1910 {
1911 return !(operator==(other));
1912 }
1913
1914 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1915 {
1916 Finley_Mesh* mesh=m_finleyMesh.get();
1917 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
1918 package, symmetry, mesh->MPIInfo);
1919 }
1920
1921 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1922 {
1923 Finley_Mesh* mesh=m_finleyMesh.get();
1924 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
1925 package, symmetry, mesh->MPIInfo);
1926 }
1927
1928 escript::Data MeshAdapter::getX() const
1929 {
1930 return continuousFunction(*this).getX();
1931 }
1932
1933 escript::Data MeshAdapter::getNormal() const
1934 {
1935 return functionOnBoundary(*this).getNormal();
1936 }
1937
1938 escript::Data MeshAdapter::getSize() const
1939 {
1940 return escript::function(*this).getSize();
1941 }
1942
1943 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1944 {
1945 int *out = NULL;
1946 Finley_Mesh* mesh=m_finleyMesh.get();
1947 switch (functionSpaceType) {
1948 case Nodes:
1949 out=mesh->Nodes->Id;
1950 break;
1951 case ReducedNodes:
1952 out=mesh->Nodes->reducedNodesId;
1953 break;
1954 case Elements:
1955 case ReducedElements:
1956 out=mesh->Elements->Id;
1957 break;
1958 case FaceElements:
1959 case ReducedFaceElements:
1960 out=mesh->FaceElements->Id;
1961 break;
1962 case Points:
1963 out=mesh->Points->Id;
1964 break;
1965 case ContactElementsZero:
1966 case ReducedContactElementsZero:
1967 case ContactElementsOne:
1968 case ReducedContactElementsOne:
1969 out=mesh->ContactElements->Id;
1970 break;
1971 case DegreesOfFreedom:
1972 out=mesh->Nodes->degreesOfFreedomId;
1973 break;
1974 case ReducedDegreesOfFreedom:
1975 out=mesh->Nodes->reducedDegreesOfFreedomId;
1976 break;
1977 default:
1978 stringstream temp;
1979 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1980 throw FinleyAdapterException(temp.str());
1981 break;
1982 }
1983 return out;
1984 }
1985 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1986 {
1987 int out=0;
1988 Finley_Mesh* mesh=m_finleyMesh.get();
1989 switch (functionSpaceType) {
1990 case Nodes:
1991 out=mesh->Nodes->Tag[sampleNo];
1992 break;
1993 case ReducedNodes:
1994 throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1995 break;
1996 case Elements:
1997 case ReducedElements:
1998 out=mesh->Elements->Tag[sampleNo];
1999 break;
2000 case FaceElements:
2001 case ReducedFaceElements:
2002 out=mesh->FaceElements->Tag[sampleNo];
2003 break;
2004 case Points:
2005 out=mesh->Points->Tag[sampleNo];
2006 break;
2007 case ContactElementsZero:
2008 case ReducedContactElementsZero:
2009 case ContactElementsOne:
2010 case ReducedContactElementsOne:
2011 out=mesh->ContactElements->Tag[sampleNo];
2012 break;
2013 case DegreesOfFreedom:
2014 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2015 break;
2016 case ReducedDegreesOfFreedom:
2017 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2018 break;
2019 default:
2020 stringstream temp;
2021 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2022 throw FinleyAdapterException(temp.str());
2023 break;
2024 }
2025 return out;
2026 }
2027
2028
2029 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2030 {
2031 Finley_Mesh* mesh=m_finleyMesh.get();
2032 escriptDataC tmp=mask.getDataC();
2033 switch(functionSpaceType) {
2034 case Nodes:
2035 mesh->Nodes->setTags(newTag, mask);
2036 break;
2037 case ReducedNodes:
2038 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2039 break;
2040 case DegreesOfFreedom:
2041 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2042 break;
2043 case ReducedDegreesOfFreedom:
2044 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2045 break;
2046 case Elements:
2047 case ReducedElements:
2048 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2049 break;
2050 case FaceElements:
2051 case ReducedFaceElements:
2052 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2053 break;
2054 case Points:
2055 Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2056 break;
2057 case ContactElementsZero:
2058 case ReducedContactElementsZero:
2059 case ContactElementsOne:
2060 case ReducedContactElementsOne:
2061 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2062 break;
2063 default:
2064 stringstream temp;
2065 temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2066 throw FinleyAdapterException(temp.str());
2067 }
2068 checkFinleyError();
2069 }
2070
2071 void MeshAdapter::setTagMap(const string& name, int tag)
2072 {
2073 Finley_Mesh* mesh=m_finleyMesh.get();
2074 Finley_Mesh_addTagMap(mesh, name.c_str(), tag);
2075 checkFinleyError();
2076 }
2077
2078 int MeshAdapter::getTag(const string& name) const
2079 {
2080 Finley_Mesh* mesh=m_finleyMesh.get();
2081 int tag = Finley_Mesh_getTag(mesh, name.c_str());
2082 checkFinleyError();
2083 return tag;
2084 }
2085
2086 bool MeshAdapter::isValidTagName(const string& name) const
2087 {
2088 Finley_Mesh* mesh=m_finleyMesh.get();
2089 return Finley_Mesh_isValidTagName(mesh, name.c_str());
2090 }
2091
2092 string MeshAdapter::showTagNames() const
2093 {
2094 stringstream temp;
2095 Finley_Mesh* mesh=m_finleyMesh.get();
2096 Finley_TagMap* tag_map=mesh->TagMap;
2097 while (tag_map) {
2098 temp << tag_map->name;
2099 tag_map=tag_map->next;
2100 if (tag_map) temp << ", ";
2101 }
2102 return temp.str();
2103 }
2104
2105 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2106 {
2107 Finley_Mesh* mesh=m_finleyMesh.get();
2108 switch(functionSpaceCode) {
2109 case Nodes:
2110 return mesh->Nodes->tagsInUse.size();
2111 case ReducedNodes:
2112 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2113 case DegreesOfFreedom:
2114 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2115 case ReducedDegreesOfFreedom:
2116 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2117 case Elements:
2118 case ReducedElements:
2119 return mesh->Elements->tagsInUse.size();
2120 case FaceElements:
2121 case ReducedFaceElements:
2122 return mesh->FaceElements->tagsInUse.size();
2123 case Points:
2124 return mesh->Points->tagsInUse.size();
2125 case ContactElementsZero:
2126 case ReducedContactElementsZero:
2127 case ContactElementsOne:
2128 case ReducedContactElementsOne:
2129 return mesh->ContactElements->tagsInUse.size();
2130 default:
2131 stringstream ss;
2132 ss << "Finley does not know anything about function space type "
2133 << functionSpaceCode;
2134 throw FinleyAdapterException(ss.str());
2135 }
2136 return 0;
2137 }
2138
2139 const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2140 {
2141 Finley_Mesh* mesh=m_finleyMesh.get();
2142 switch(functionSpaceCode) {
2143 case Nodes:
2144 return &mesh->Nodes->tagsInUse[0];
2145 case ReducedNodes:
2146 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2147 case DegreesOfFreedom:
2148 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2149 case ReducedDegreesOfFreedom:
2150 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2151 case Elements:
2152 case ReducedElements:
2153 return &mesh->Elements->tagsInUse[0];
2154 case FaceElements:
2155 case ReducedFaceElements:
2156 return &mesh->FaceElements->tagsInUse[0];
2157 case Points:
2158 return &mesh->Points->tagsInUse[0];
2159 case ContactElementsZero:
2160 case ReducedContactElementsZero:
2161 case ContactElementsOne:
2162 case ReducedContactElementsOne:
2163 return &mesh->ContactElements->tagsInUse[0];
2164 default:
2165 stringstream temp;
2166 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2167 throw FinleyAdapterException(temp.str());
2168 }
2169 return NULL;
2170 }
2171
2172
2173 bool MeshAdapter::canTag(int functionSpaceCode) const
2174 {
2175 switch(functionSpaceCode) {
2176 case Nodes:
2177 case Elements:
2178 case ReducedElements:
2179 case FaceElements:
2180 case ReducedFaceElements:
2181 case Points:
2182 case ContactElementsZero:
2183 case ReducedContactElementsZero:
2184 case ContactElementsOne:
2185 case ReducedContactElementsOne:
2186 return true;
2187 default:
2188 return false;
2189 }
2190 }
2191
2192 escript::AbstractDomain::StatusType MeshAdapter::getStatus() const
2193 {
2194 Finley_Mesh* mesh=m_finleyMesh.get();
2195 return Finley_Mesh_getStatus(mesh);
2196 }
2197
2198 int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2199 {
2200 Finley_Mesh* mesh=m_finleyMesh.get();
2201 int order =-1;
2202 switch(functionSpaceCode) {
2203 case Nodes:
2204 case DegreesOfFreedom:
2205 order=mesh->approximationOrder;
2206 break;
2207 case ReducedNodes:
2208 case ReducedDegreesOfFreedom:
2209 order=mesh->reducedApproximationOrder;
2210 break;
2211 case Elements:
2212 case FaceElements:
2213 case Points:
2214 case ContactElementsZero:
2215 case ContactElementsOne:
2216 order=mesh->integrationOrder;
2217 break;
2218 case ReducedElements:
2219 case ReducedFaceElements:
2220 case ReducedContactElementsZero:
2221 case ReducedContactElementsOne:
2222 order=mesh->reducedIntegrationOrder;
2223 break;
2224 default:
2225 stringstream temp;
2226 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2227 throw FinleyAdapterException(temp.str());
2228 }
2229 return order;
2230 }
2231
2232 bool MeshAdapter::supportsContactElements() const
2233 {
2234 return true;
2235 }
2236
2237 ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id,
2238 index_t order, index_t reducedOrder)
2239 {
2240 m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2241 }
2242
2243 ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2244 {
2245 Finley_ReferenceElementSet_dealloc(m_refSet);
2246 }
2247
2248 void MeshAdapter::addDiracPoints(const vector<double>& points,
2249 const vector<int>& tags) const
2250 {
2251 // points will be flattened
2252 const int dim = getDim();
2253 int numPoints=points.size()/dim;
2254 int numTags=tags.size();
2255 Finley_Mesh* mesh=m_finleyMesh.get();
2256
2257 if ( points.size() % dim != 0 ) {
2258 throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2259 }
2260
2261 if ((numTags > 0) && (numPoints != numTags))
2262 throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2263
2264 Finley_Mesh_addPoints(mesh, numPoints, &points[0], &tags[0]);
2265 checkFinleyError();
2266 }
2267
2268 // void MeshAdapter::addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2269 // {
2270 // const int dim = getDim();
2271 // int numPoints=boost::python::extract<int>(points.attr("__len__")());
2272 // int numTags=boost::python::extract<int>(tags.attr("__len__")());
2273 // Finley_Mesh* mesh=m_finleyMesh.get();
2274 //
2275 // if ( (numTags > 0) && ( numPoints != numTags ) )
2276 // throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2277 //
2278 // double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2279 // int* tags_ptr= TMPMEMALLOC(numPoints, int);
2280 //
2281 // for (int i=0;i<numPoints;++i) {
2282 // int tag_id=-1;
2283 // int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2284 // if ( numComps != dim ) {
2285 // stringstream temp;
2286 // temp << "Error - illegal number of components " << numComps << " for point " << i;
2287 // throw FinleyAdapterException(temp.str());
2288 // }
2289 // points_ptr[ i * dim ] = boost::python::extract<double>(points[i][0]);
2290 // if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2291 // if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2292 //
2293 // if ( numTags > 0) {
2294 // boost::python::extract<string> ex_str(tags[i]);
2295 // if ( ex_str.check() ) {
2296 // tag_id=getTag( ex_str());
2297 // } else {
2298 // boost::python::extract<int> ex_int(tags[i]);
2299 // if ( ex_int.check() ) {
2300 // tag_id=ex_int();
2301 // } else {
2302 // stringstream temp;
2303 // temp << "Error - unable to extract tag for point " << i;
2304 // throw FinleyAdapterException(temp.str());
2305 // }
2306 // }
2307 // }
2308 // tags_ptr[i]=tag_id;
2309 // }
2310 //
2311 // Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2312 // checkPasoError();
2313 //
2314 // TMPMEMFREE(points_ptr);
2315 // TMPMEMFREE(tags_ptr);
2316 // }
2317
2318 /*
2319 void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2320 {
2321 boost::python::list points = boost::python::list();
2322 boost::python::list tags = boost::python::list();
2323 points.append(point);
2324 tags.append(tag);
2325 addDiracPoints(points, tags);
2326 }
2327 */
2328
2329 /*
2330 void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const string& tag) const
2331 {
2332 boost::python::list points = boost::python::list();
2333 boost::python::list tags = boost::python::list();
2334 points.append(point);
2335 tags.append(tag);
2336 addDiracPoints(points, tags);
2337 }
2338 */
2339 } // end of namespace
2340

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26