/[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 1754 - (show annotations)
Sun Sep 7 22:07:26 2008 UTC (10 years, 6 months ago) by ksteube
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 84431 byte(s)
Added new test suite run_inputOutput.py to systematically test I/O.
Can determine if two domains are same with Fourier analysis.
Added new method getNumDataPointsGlobal to return number of DPs of a distributed mesh.
Reading of tags in ReadMeshMPI failed occasionally, should be more robust now.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26