/[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 4408 - (show annotations)
Tue May 14 06:58:43 2013 UTC (6 years, 1 month ago) by caltinay
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 94960 byte(s)
Mostly no-op - trying to get some consistency in style while going through a
few files.
Also some minor changes related to C->C++ (e.g. no need to copy vector to pass
to functions etc.)

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26