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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1748 - (show annotations)
Wed Sep 3 06:10:39 2008 UTC (10 years, 8 months ago) by ksteube
File size: 84232 byte(s)
MPI parallelism for Data().dump and load.  Use multiple NetCDF
files, one file per MPI process

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 per sample and the number of samples
652 // needed to represent data on a parts of the mesh.
653 //
654 pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
655 {
656 int numDataPointsPerSample=0;
657 int numSamples=0;
658 Finley_Mesh* mesh=m_finleyMesh.get();
659 switch (functionSpaceCode) {
660 case(Nodes):
661 numDataPointsPerSample=1;
662 numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);
663 break;
664 case(ReducedNodes):
665 numDataPointsPerSample=1;
666 numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
667 break;
668 case(Elements):
669 if (mesh->Elements!=NULL) {
670 numSamples=mesh->Elements->numElements;
671 numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
672 }
673 break;
674 case(ReducedElements):
675 if (mesh->Elements!=NULL) {
676 numSamples=mesh->Elements->numElements;
677 numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;
678 }
679 break;
680 case(FaceElements):
681 if (mesh->FaceElements!=NULL) {
682 numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
683 numSamples=mesh->FaceElements->numElements;
684 }
685 break;
686 case(ReducedFaceElements):
687 if (mesh->FaceElements!=NULL) {
688 numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;
689 numSamples=mesh->FaceElements->numElements;
690 }
691 break;
692 case(Points):
693 if (mesh->Points!=NULL) {
694 numDataPointsPerSample=1;
695 numSamples=mesh->Points->numElements;
696 }
697 break;
698 case(ContactElementsZero):
699 if (mesh->ContactElements!=NULL) {
700 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
701 numSamples=mesh->ContactElements->numElements;
702 }
703 break;
704 case(ReducedContactElementsZero):
705 if (mesh->ContactElements!=NULL) {
706 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
707 numSamples=mesh->ContactElements->numElements;
708 }
709 break;
710 case(ContactElementsOne):
711 if (mesh->ContactElements!=NULL) {
712 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
713 numSamples=mesh->ContactElements->numElements;
714 }
715 break;
716 case(ReducedContactElementsOne):
717 if (mesh->ContactElements!=NULL) {
718 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
719 numSamples=mesh->ContactElements->numElements;
720 }
721 break;
722 case(DegreesOfFreedom):
723 if (mesh->Nodes!=NULL) {
724 numDataPointsPerSample=1;
725 numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
726 }
727 break;
728 case(ReducedDegreesOfFreedom):
729 if (mesh->Nodes!=NULL) {
730 numDataPointsPerSample=1;
731 numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
732 }
733 break;
734 default:
735 stringstream temp;
736 temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
737 throw FinleyAdapterException(temp.str());
738 break;
739 }
740 return pair<int,int>(numDataPointsPerSample,numSamples);
741 }
742
743 //
744 // adds linear PDE of second order into a given stiffness matrix and right hand side:
745 //
746 void MeshAdapter::addPDEToSystem(
747 SystemMatrixAdapter& mat, escript::Data& rhs,
748 const escript::Data& A, const escript::Data& B, const escript::Data& C,const escript::Data& D,const escript::Data& X,const escript::Data& Y,
749 const escript::Data& d, const escript::Data& y,
750 const escript::Data& d_contact,const escript::Data& y_contact) const
751 {
752 escriptDataC _rhs=rhs.getDataC();
753 escriptDataC _A =A.getDataC();
754 escriptDataC _B=B.getDataC();
755 escriptDataC _C=C.getDataC();
756 escriptDataC _D=D.getDataC();
757 escriptDataC _X=X.getDataC();
758 escriptDataC _Y=Y.getDataC();
759 escriptDataC _d=d.getDataC();
760 escriptDataC _y=y.getDataC();
761 escriptDataC _d_contact=d_contact.getDataC();
762 escriptDataC _y_contact=y_contact.getDataC();
763
764 Finley_Mesh* mesh=m_finleyMesh.get();
765
766 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
767 checkFinleyError();
768
769 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
770 checkFinleyError();
771
772 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
773 checkFinleyError();
774 }
775
776 void MeshAdapter::addPDEToLumpedSystem(
777 escript::Data& mat,
778 const escript::Data& D,
779 const escript::Data& d) const
780 {
781 escriptDataC _mat=mat.getDataC();
782 escriptDataC _D=D.getDataC();
783 escriptDataC _d=d.getDataC();
784
785 Finley_Mesh* mesh=m_finleyMesh.get();
786
787 Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
788 Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
789
790 checkFinleyError();
791 }
792
793
794 //
795 // adds linear PDE of second order into the right hand side only
796 //
797 void MeshAdapter::addPDEToRHS( escript::Data& rhs, const escript::Data& X,const escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
798 {
799 Finley_Mesh* mesh=m_finleyMesh.get();
800
801 escriptDataC _rhs=rhs.getDataC();
802 escriptDataC _X=X.getDataC();
803 escriptDataC _Y=Y.getDataC();
804 escriptDataC _y=y.getDataC();
805 escriptDataC _y_contact=y_contact.getDataC();
806
807 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
808 checkFinleyError();
809
810 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
811 checkFinleyError();
812
813 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
814 checkFinleyError();
815 }
816 //
817 // adds PDE of second order into a transport problem
818 //
819 void MeshAdapter::addPDEToTransportProblem(
820 TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,
821 const escript::Data& A, const escript::Data& B, const escript::Data& C,
822 const escript::Data& D,const escript::Data& X,const escript::Data& Y,
823 const escript::Data& d, const escript::Data& y,
824 const escript::Data& d_contact,const escript::Data& y_contact) const
825 {
826 DataArrayView::ShapeType shape;
827 source.expand();
828 escriptDataC _source=source.getDataC();
829 escriptDataC _M=M.getDataC();
830 escriptDataC _A=A.getDataC();
831 escriptDataC _B=B.getDataC();
832 escriptDataC _C=C.getDataC();
833 escriptDataC _D=D.getDataC();
834 escriptDataC _X=X.getDataC();
835 escriptDataC _Y=Y.getDataC();
836 escriptDataC _d=d.getDataC();
837 escriptDataC _y=y.getDataC();
838 escriptDataC _d_contact=d_contact.getDataC();
839 escriptDataC _y_contact=y_contact.getDataC();
840
841 Finley_Mesh* mesh=m_finleyMesh.get();
842 Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();
843
844 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
845 checkFinleyError();
846
847 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
848 checkFinleyError();
849
850 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
851 checkFinleyError();
852
853 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
854 checkFinleyError();
855 }
856
857 //
858 // interpolates data between different function spaces:
859 //
860 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
861 {
862 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
863 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
864 if (inDomain!=*this)
865 throw FinleyAdapterException("Error - Illegal domain of interpolant.");
866 if (targetDomain!=*this)
867 throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
868
869 Finley_Mesh* mesh=m_finleyMesh.get();
870 escriptDataC _target=target.getDataC();
871 escriptDataC _in=in.getDataC();
872 switch(in.getFunctionSpace().getTypeCode()) {
873 case(Nodes):
874 switch(target.getFunctionSpace().getTypeCode()) {
875 case(Nodes):
876 case(ReducedNodes):
877 case(DegreesOfFreedom):
878 case(ReducedDegreesOfFreedom):
879 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
880 break;
881 case(Elements):
882 case(ReducedElements):
883 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
884 break;
885 case(FaceElements):
886 case(ReducedFaceElements):
887 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
888 break;
889 case(Points):
890 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
891 break;
892 case(ContactElementsZero):
893 case(ReducedContactElementsZero):
894 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
895 break;
896 case(ContactElementsOne):
897 case(ReducedContactElementsOne):
898 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
899 break;
900 default:
901 stringstream temp;
902 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
903 throw FinleyAdapterException(temp.str());
904 break;
905 }
906 break;
907 case(ReducedNodes):
908 switch(target.getFunctionSpace().getTypeCode()) {
909 case(Nodes):
910 case(ReducedNodes):
911 case(DegreesOfFreedom):
912 case(ReducedDegreesOfFreedom):
913 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
914 break;
915 case(Elements):
916 case(ReducedElements):
917 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
918 break;
919 case(FaceElements):
920 case(ReducedFaceElements):
921 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
922 break;
923 case(Points):
924 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
925 break;
926 case(ContactElementsZero):
927 case(ReducedContactElementsZero):
928 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
929 break;
930 case(ContactElementsOne):
931 case(ReducedContactElementsOne):
932 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
933 break;
934 default:
935 stringstream temp;
936 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
937 throw FinleyAdapterException(temp.str());
938 break;
939 }
940 break;
941 case(Elements):
942 if (target.getFunctionSpace().getTypeCode()==Elements) {
943 Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
944 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
945 Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
946 } else {
947 throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
948 }
949 break;
950 case(ReducedElements):
951 if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
952 Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
953 } else {
954 throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
955 }
956 break;
957 case(FaceElements):
958 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
959 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
960 } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
961 Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
962 } else {
963 throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
964 }
965 break;
966 case(ReducedFaceElements):
967 if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
968 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
969 } else {
970 throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
971 }
972 break;
973 case(Points):
974 if (target.getFunctionSpace().getTypeCode()==Points) {
975 Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
976 } else {
977 throw FinleyAdapterException("Error - No interpolation with data on points possible.");
978 }
979 break;
980 case(ContactElementsZero):
981 case(ContactElementsOne):
982 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
983 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
984 } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
985 Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
986 } else {
987 throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
988 }
989 break;
990 case(ReducedContactElementsZero):
991 case(ReducedContactElementsOne):
992 if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
993 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
994 } else {
995 throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
996 }
997 break;
998 case(DegreesOfFreedom):
999 switch(target.getFunctionSpace().getTypeCode()) {
1000 case(ReducedDegreesOfFreedom):
1001 case(DegreesOfFreedom):
1002 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1003 break;
1004
1005 case(Nodes):
1006 case(ReducedNodes):
1007 if (getMPISize()>1) {
1008 escript::Data temp=escript::Data(in);
1009 temp.expand();
1010 escriptDataC _in2 = temp.getDataC();
1011 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1012 } else {
1013 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1014 }
1015 break;
1016 case(Elements):
1017 case(ReducedElements):
1018 if (getMPISize()>1) {
1019 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
1020 escriptDataC _in2 = temp.getDataC();
1021 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1022 } else {
1023 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1024 }
1025 break;
1026 case(FaceElements):
1027 case(ReducedFaceElements):
1028 if (getMPISize()>1) {
1029 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
1030 escriptDataC _in2 = temp.getDataC();
1031 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1032
1033 } else {
1034 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1035 }
1036 break;
1037 case(Points):
1038 if (getMPISize()>1) {
1039 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
1040 escriptDataC _in2 = temp.getDataC();
1041 } else {
1042 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1043 }
1044 break;
1045 case(ContactElementsZero):
1046 case(ContactElementsOne):
1047 case(ReducedContactElementsZero):
1048 case(ReducedContactElementsOne):
1049 if (getMPISize()>1) {
1050 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
1051 escriptDataC _in2 = temp.getDataC();
1052 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1053 } else {
1054 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1055 }
1056 break;
1057 default:
1058 stringstream temp;
1059 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1060 throw FinleyAdapterException(temp.str());
1061 break;
1062 }
1063 break;
1064 case(ReducedDegreesOfFreedom):
1065 switch(target.getFunctionSpace().getTypeCode()) {
1066 case(Nodes):
1067 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1068 break;
1069 case(ReducedNodes):
1070 if (getMPISize()>1) {
1071 escript::Data temp=escript::Data(in);
1072 temp.expand();
1073 escriptDataC _in2 = temp.getDataC();
1074 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1075 } else {
1076 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1077 }
1078 break;
1079 case(DegreesOfFreedom):
1080 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1081 break;
1082 case(ReducedDegreesOfFreedom):
1083 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1084 break;
1085 case(Elements):
1086 case(ReducedElements):
1087 if (getMPISize()>1) {
1088 escript::Data temp=escript::Data( in, reducedContinuousFunction(asAbstractContinuousDomain()) );
1089 escriptDataC _in2 = temp.getDataC();
1090 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1091 } else {
1092 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1093 }
1094 break;
1095 case(FaceElements):
1096 case(ReducedFaceElements):
1097 if (getMPISize()>1) {
1098 escript::Data temp=escript::Data( in, reducedContinuousFunction(asAbstractContinuousDomain()) );
1099 escriptDataC _in2 = temp.getDataC();
1100 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1101 } else {
1102 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1103 }
1104 break;
1105 case(Points):
1106 if (getMPISize()>1) {
1107 escript::Data temp=escript::Data( in, reducedContinuousFunction(asAbstractContinuousDomain()) );
1108 escriptDataC _in2 = temp.getDataC();
1109 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1110 } else {
1111 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1112 }
1113 break;
1114 case(ContactElementsZero):
1115 case(ContactElementsOne):
1116 case(ReducedContactElementsZero):
1117 case(ReducedContactElementsOne):
1118 if (getMPISize()>1) {
1119 escript::Data temp=escript::Data( in, reducedContinuousFunction(asAbstractContinuousDomain()) );
1120 escriptDataC _in2 = temp.getDataC();
1121 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1122 } else {
1123 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1124 }
1125 break;
1126 default:
1127 stringstream temp;
1128 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1129 throw FinleyAdapterException(temp.str());
1130 break;
1131 }
1132 break;
1133 default:
1134 stringstream temp;
1135 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1136 throw FinleyAdapterException(temp.str());
1137 break;
1138 }
1139 checkFinleyError();
1140 }
1141
1142 //
1143 // copies the locations of sample points into x:
1144 //
1145 void MeshAdapter::setToX(escript::Data& arg) const
1146 {
1147 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
1148 if (argDomain!=*this)
1149 throw FinleyAdapterException("Error - Illegal domain of data point locations");
1150 Finley_Mesh* mesh=m_finleyMesh.get();
1151 // in case of values node coordinates we can do the job directly:
1152 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1153 escriptDataC _arg=arg.getDataC();
1154 Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1155 } else {
1156 escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
1157 escriptDataC _tmp_data=tmp_data.getDataC();
1158 Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1159 // this is then interpolated onto arg:
1160 interpolateOnDomain(arg,tmp_data);
1161 }
1162 checkFinleyError();
1163 }
1164
1165 //
1166 // return the normal vectors at the location of data points as a Data object:
1167 //
1168 void MeshAdapter::setToNormal(escript::Data& normal) const
1169 {
1170 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
1171 if (normalDomain!=*this)
1172 throw FinleyAdapterException("Error - Illegal domain of normal locations");
1173 Finley_Mesh* mesh=m_finleyMesh.get();
1174 escriptDataC _normal=normal.getDataC();
1175 switch(normal.getFunctionSpace().getTypeCode()) {
1176 case(Nodes):
1177 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1178 break;
1179 case(ReducedNodes):
1180 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1181 break;
1182 case(Elements):
1183 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1184 break;
1185 case(ReducedElements):
1186 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1187 break;
1188 case (FaceElements):
1189 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1190 break;
1191 case (ReducedFaceElements):
1192 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1193 break;
1194 case(Points):
1195 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1196 break;
1197 case (ContactElementsOne):
1198 case (ContactElementsZero):
1199 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1200 break;
1201 case (ReducedContactElementsOne):
1202 case (ReducedContactElementsZero):
1203 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1204 break;
1205 case(DegreesOfFreedom):
1206 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1207 break;
1208 case(ReducedDegreesOfFreedom):
1209 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1210 break;
1211 default:
1212 stringstream temp;
1213 temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1214 throw FinleyAdapterException(temp.str());
1215 break;
1216 }
1217 checkFinleyError();
1218 }
1219
1220 //
1221 // interpolates data to other domain:
1222 //
1223 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1224 {
1225 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
1226 if (targetDomain!=*this)
1227 throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1228
1229 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
1230 }
1231
1232 //
1233 // calculates the integral of a function defined of arg:
1234 //
1235 void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
1236 {
1237 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
1238 if (argDomain!=*this)
1239 throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1240
1241 double blocktimer_start = blocktimer_time();
1242 Finley_Mesh* mesh=m_finleyMesh.get();
1243 escriptDataC _temp;
1244 escript::Data temp;
1245 escriptDataC _arg=arg.getDataC();
1246 switch(arg.getFunctionSpace().getTypeCode()) {
1247 case(Nodes):
1248 temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
1249 _temp=temp.getDataC();
1250 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1251 break;
1252 case(ReducedNodes):
1253 temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
1254 _temp=temp.getDataC();
1255 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1256 break;
1257 case(Elements):
1258 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1259 break;
1260 case(ReducedElements):
1261 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1262 break;
1263 case(FaceElements):
1264 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1265 break;
1266 case(ReducedFaceElements):
1267 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1268 break;
1269 case(Points):
1270 throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1271 break;
1272 case(ContactElementsZero):
1273 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1274 break;
1275 case(ReducedContactElementsZero):
1276 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1277 break;
1278 case(ContactElementsOne):
1279 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1280 break;
1281 case(ReducedContactElementsOne):
1282 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1283 break;
1284 case(DegreesOfFreedom):
1285 temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
1286 _temp=temp.getDataC();
1287 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1288 break;
1289 case(ReducedDegreesOfFreedom):
1290 temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
1291 _temp=temp.getDataC();
1292 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1293 break;
1294 default:
1295 stringstream temp;
1296 temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1297 throw FinleyAdapterException(temp.str());
1298 break;
1299 }
1300 checkFinleyError();
1301 blocktimer_increment("integrate()", blocktimer_start);
1302 }
1303
1304 //
1305 // calculates the gradient of arg:
1306 //
1307 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1308 {
1309 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
1310 if (argDomain!=*this)
1311 throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1312 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());
1313 if (gradDomain!=*this)
1314 throw FinleyAdapterException("Error - Illegal domain of gradient");
1315
1316 Finley_Mesh* mesh=m_finleyMesh.get();
1317 escriptDataC _grad=grad.getDataC();
1318 escriptDataC nodeDataC;
1319 escript::Data temp;
1320 if (getMPISize()>1) {
1321 if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1322 temp=escript::Data( arg, continuousFunction(asAbstractContinuousDomain()) );
1323 nodeDataC = temp.getDataC();
1324 } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1325 temp=escript::Data( arg, reducedContinuousFunction(asAbstractContinuousDomain()) );
1326 nodeDataC = temp.getDataC();
1327 } else {
1328 nodeDataC = arg.getDataC();
1329 }
1330 } else {
1331 nodeDataC = arg.getDataC();
1332 }
1333 switch(grad.getFunctionSpace().getTypeCode()) {
1334 case(Nodes):
1335 throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1336 break;
1337 case(ReducedNodes):
1338 throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1339 break;
1340 case(Elements):
1341 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1342 break;
1343 case(ReducedElements):
1344 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1345 break;
1346 case(FaceElements):
1347 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1348 break;
1349 case(ReducedFaceElements):
1350 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1351 break;
1352 case(Points):
1353 throw FinleyAdapterException("Error - Gradient at points is not supported.");
1354 break;
1355 case(ContactElementsZero):
1356 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1357 break;
1358 case(ReducedContactElementsZero):
1359 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1360 break;
1361 case(ContactElementsOne):
1362 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1363 break;
1364 case(ReducedContactElementsOne):
1365 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1366 break;
1367 case(DegreesOfFreedom):
1368 throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1369 break;
1370 case(ReducedDegreesOfFreedom):
1371 throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1372 break;
1373 default:
1374 stringstream temp;
1375 temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1376 throw FinleyAdapterException(temp.str());
1377 break;
1378 }
1379 checkFinleyError();
1380 }
1381
1382 //
1383 // returns the size of elements:
1384 //
1385 void MeshAdapter::setToSize(escript::Data& size) const
1386 {
1387 Finley_Mesh* mesh=m_finleyMesh.get();
1388 escriptDataC tmp=size.getDataC();
1389 switch(size.getFunctionSpace().getTypeCode()) {
1390 case(Nodes):
1391 throw FinleyAdapterException("Error - Size of nodes is not supported.");
1392 break;
1393 case(ReducedNodes):
1394 throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1395 break;
1396 case(Elements):
1397 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1398 break;
1399 case(ReducedElements):
1400 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1401 break;
1402 case(FaceElements):
1403 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1404 break;
1405 case(ReducedFaceElements):
1406 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1407 break;
1408 case(Points):
1409 throw FinleyAdapterException("Error - Size of point elements is not supported.");
1410 break;
1411 case(ContactElementsZero):
1412 case(ContactElementsOne):
1413 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1414 break;
1415 case(ReducedContactElementsZero):
1416 case(ReducedContactElementsOne):
1417 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1418 break;
1419 case(DegreesOfFreedom):
1420 throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1421 break;
1422 case(ReducedDegreesOfFreedom):
1423 throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1424 break;
1425 default:
1426 stringstream temp;
1427 temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1428 throw FinleyAdapterException(temp.str());
1429 break;
1430 }
1431 checkFinleyError();
1432 }
1433
1434 // sets the location of nodes:
1435 void MeshAdapter::setNewX(const escript::Data& new_x)
1436 {
1437 Finley_Mesh* mesh=m_finleyMesh.get();
1438 escriptDataC tmp;
1439 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
1440 if (newDomain!=*this)
1441 throw FinleyAdapterException("Error - Illegal domain of new point locations");
1442 tmp = new_x.getDataC();
1443 Finley_Mesh_setCoordinates(mesh,&tmp);
1444 checkFinleyError();
1445 }
1446
1447 // saves a data array in openDX format:
1448 void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
1449 {
1450 unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;
1451 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
1452 /* win32 refactor */
1453 char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1454 for(int i=0;i<num_data;i++)
1455 {
1456 names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
1457 }
1458
1459 escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
1460 escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
1461
1462 boost::python::list keys=arg.keys();
1463 for (int i=0;i<num_data;++i) {
1464 std::string n=boost::python::extract<std::string>(keys[i]);
1465 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1466 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
1467 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
1468 data[i]=d.getDataC();
1469 ptr_data[i]= &(data[i]);
1470 if (n.length()>MAX_namelength-1) {
1471 strncpy(names[i],n.c_str(),MAX_namelength-1);
1472 names[i][MAX_namelength-1]='\0';
1473 } else {
1474 strcpy(names[i],n.c_str());
1475 }
1476 }
1477 Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);
1478 checkFinleyError();
1479
1480 /* win32 refactor */
1481 TMPMEMFREE(data);
1482 TMPMEMFREE(ptr_data);
1483 for(int i=0;i<num_data;i++)
1484 {
1485 TMPMEMFREE(names[i]);
1486 }
1487 TMPMEMFREE(names);
1488
1489 return;
1490 }
1491
1492 // saves a data array in openVTK format:
1493 void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
1494 {
1495 unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;
1496 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
1497 /* win32 refactor */
1498 char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1499 for(int i=0;i<num_data;i++)
1500 {
1501 names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
1502 }
1503
1504 escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
1505 escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
1506
1507 boost::python::list keys=arg.keys();
1508 for (int i=0;i<num_data;++i) {
1509 std::string n=boost::python::extract<std::string>(keys[i]);
1510 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1511 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
1512 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
1513 data[i]=d.getDataC();
1514 ptr_data[i]=&(data[i]);
1515 if (n.length()>MAX_namelength-1) {
1516 strncpy(names[i],n.c_str(),MAX_namelength-1);
1517 names[i][MAX_namelength-1]='\0';
1518 } else {
1519 strcpy(names[i],n.c_str());
1520 }
1521 }
1522 Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);
1523
1524 checkFinleyError();
1525 /* win32 refactor */
1526 TMPMEMFREE(data);
1527 TMPMEMFREE(ptr_data);
1528 for(int i=0;i<num_data;i++)
1529 {
1530 TMPMEMFREE(names[i]);
1531 }
1532 TMPMEMFREE(names);
1533
1534 return;
1535 }
1536
1537
1538 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
1539 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1540 const int row_blocksize,
1541 const escript::FunctionSpace& row_functionspace,
1542 const int column_blocksize,
1543 const escript::FunctionSpace& column_functionspace,
1544 const int type) const
1545 {
1546 int reduceRowOrder=0;
1547 int reduceColOrder=0;
1548 // is the domain right?
1549 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());
1550 if (row_domain!=*this)
1551 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1552 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());
1553 if (col_domain!=*this)
1554 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1555 // is the function space type right
1556 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1557 reduceRowOrder=0;
1558 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1559 reduceRowOrder=1;
1560 } else {
1561 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1562 }
1563 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1564 reduceColOrder=0;
1565 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1566 reduceColOrder=1;
1567 } else {
1568 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1569 }
1570 // generate matrix:
1571
1572 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1573 checkFinleyError();
1574 Paso_SystemMatrix* fsystemMatrix;
1575 int trilinos = 0;
1576 if (trilinos) {
1577 #ifdef TRILINOS
1578 /* Allocation Epetra_VrbMatrix here */
1579 #endif
1580 }
1581 else {
1582 fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
1583 }
1584 checkPasoError();
1585 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1586 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1587 }
1588 // creates a TransportProblemAdapter
1589 TransportProblemAdapter MeshAdapter::newTransportProblem(
1590 const double theta,
1591 const int blocksize,
1592 const escript::FunctionSpace& functionspace,
1593 const int type) const
1594 {
1595 int reduceOrder=0;
1596 // is the domain right?
1597 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(functionspace.getDomain());
1598 if (domain!=*this)
1599 throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1600 // is the function space type right
1601 if (functionspace.getTypeCode()==DegreesOfFreedom) {
1602 reduceOrder=0;
1603 } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1604 reduceOrder=1;
1605 } else {
1606 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1607 }
1608 // generate matrix:
1609
1610 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1611 checkFinleyError();
1612 Paso_FCTransportProblem* transportProblem;
1613 transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);
1614 checkPasoError();
1615 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1616 return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);
1617 }
1618
1619 //
1620 // vtkObject MeshAdapter::createVtkObject() const
1621 // TODO:
1622 //
1623
1624 //
1625 // returns true if data at the atom_type is considered as being cell centered:
1626 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1627 {
1628 switch(functionSpaceCode) {
1629 case(Nodes):
1630 case(DegreesOfFreedom):
1631 case(ReducedDegreesOfFreedom):
1632 return false;
1633 break;
1634 case(Elements):
1635 case(FaceElements):
1636 case(Points):
1637 case(ContactElementsZero):
1638 case(ContactElementsOne):
1639 case(ReducedElements):
1640 case(ReducedFaceElements):
1641 case(ReducedContactElementsZero):
1642 case(ReducedContactElementsOne):
1643 return true;
1644 break;
1645 default:
1646 stringstream temp;
1647 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1648 throw FinleyAdapterException(temp.str());
1649 break;
1650 }
1651 checkFinleyError();
1652 return false;
1653 }
1654
1655 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1656 {
1657 switch(functionSpaceType_source) {
1658 case(Nodes):
1659 switch(functionSpaceType_target) {
1660 case(Nodes):
1661 case(ReducedNodes):
1662 case(ReducedDegreesOfFreedom):
1663 case(DegreesOfFreedom):
1664 case(Elements):
1665 case(ReducedElements):
1666 case(FaceElements):
1667 case(ReducedFaceElements):
1668 case(Points):
1669 case(ContactElementsZero):
1670 case(ReducedContactElementsZero):
1671 case(ContactElementsOne):
1672 case(ReducedContactElementsOne):
1673 return true;
1674 default:
1675 stringstream temp;
1676 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1677 throw FinleyAdapterException(temp.str());
1678 }
1679 break;
1680 case(ReducedNodes):
1681 switch(functionSpaceType_target) {
1682 case(ReducedNodes):
1683 case(ReducedDegreesOfFreedom):
1684 case(Elements):
1685 case(ReducedElements):
1686 case(FaceElements):
1687 case(ReducedFaceElements):
1688 case(Points):
1689 case(ContactElementsZero):
1690 case(ReducedContactElementsZero):
1691 case(ContactElementsOne):
1692 case(ReducedContactElementsOne):
1693 return true;
1694 case(Nodes):
1695 case(DegreesOfFreedom):
1696 return false;
1697 default:
1698 stringstream temp;
1699 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1700 throw FinleyAdapterException(temp.str());
1701 }
1702 break;
1703 case(Elements):
1704 if (functionSpaceType_target==Elements) {
1705 return true;
1706 } else if (functionSpaceType_target==ReducedElements) {
1707 return true;
1708 } else {
1709 return false;
1710 }
1711 case(ReducedElements):
1712 if (functionSpaceType_target==ReducedElements) {
1713 return true;
1714 } else {
1715 return false;
1716 }
1717 case(FaceElements):
1718 if (functionSpaceType_target==FaceElements) {
1719 return true;
1720 } else if (functionSpaceType_target==ReducedFaceElements) {
1721 return true;
1722 } else {
1723 return false;
1724 }
1725 case(ReducedFaceElements):
1726 if (functionSpaceType_target==ReducedFaceElements) {
1727 return true;
1728 } else {
1729 return false;
1730 }
1731 case(Points):
1732 if (functionSpaceType_target==Points) {
1733 return true;
1734 } else {
1735 return false;
1736 }
1737 case(ContactElementsZero):
1738 case(ContactElementsOne):
1739 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1740 return true;
1741 } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1742 return true;
1743 } else {
1744 return false;
1745 }
1746 case(ReducedContactElementsZero):
1747 case(ReducedContactElementsOne):
1748 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1749 return true;
1750 } else {
1751 return false;
1752 }
1753 case(DegreesOfFreedom):
1754 switch(functionSpaceType_target) {
1755 case(ReducedDegreesOfFreedom):
1756 case(DegreesOfFreedom):
1757 case(Nodes):
1758 case(ReducedNodes):
1759 case(Elements):
1760 case(ReducedElements):
1761 case(Points):
1762 case(FaceElements):
1763 case(ReducedFaceElements):
1764 case(ContactElementsZero):
1765 case(ReducedContactElementsZero):
1766 case(ContactElementsOne):
1767 case(ReducedContactElementsOne):
1768 return true;
1769 default:
1770 stringstream temp;
1771 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1772 throw FinleyAdapterException(temp.str());
1773 }
1774 break;
1775 case(ReducedDegreesOfFreedom):
1776 switch(functionSpaceType_target) {
1777 case(ReducedDegreesOfFreedom):
1778 case(ReducedNodes):
1779 case(Elements):
1780 case(ReducedElements):
1781 case(FaceElements):
1782 case(ReducedFaceElements):
1783 case(Points):
1784 case(ContactElementsZero):
1785 case(ReducedContactElementsZero):
1786 case(ContactElementsOne):
1787 case(ReducedContactElementsOne):
1788 return true;
1789 case(Nodes):
1790 case(DegreesOfFreedom):
1791 return false;
1792 default:
1793 stringstream temp;
1794 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1795 throw FinleyAdapterException(temp.str());
1796 }
1797 break;
1798 default:
1799 stringstream temp;
1800 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1801 throw FinleyAdapterException(temp.str());
1802 break;
1803 }
1804 checkFinleyError();
1805 return false;
1806 }
1807
1808 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1809 {
1810 return false;
1811 }
1812
1813 bool MeshAdapter::operator==(const AbstractDomain& other) const
1814 {
1815 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1816 if (temp!=0) {
1817 return (m_finleyMesh==temp->m_finleyMesh);
1818 } else {
1819 return false;
1820 }
1821 }
1822
1823 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1824 {
1825 return !(operator==(other));
1826 }
1827
1828 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1829 {
1830 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1831 checkPasoError();
1832 return out;
1833 }
1834
1835 escript::Data MeshAdapter::getX() const
1836 {
1837 return continuousFunction(asAbstractContinuousDomain()).getX();
1838 }
1839
1840 escript::Data MeshAdapter::getNormal() const
1841 {
1842 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1843 }
1844
1845 escript::Data MeshAdapter::getSize() const
1846 {
1847 return function(asAbstractContinuousDomain()).getSize();
1848 }
1849
1850 int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1851 {
1852 int *out = NULL;
1853 Finley_Mesh* mesh=m_finleyMesh.get();
1854 switch (functionSpaceType) {
1855 case(Nodes):
1856 out=mesh->Nodes->Id;
1857 break;
1858 case(ReducedNodes):
1859 out=mesh->Nodes->reducedNodesId;
1860 break;
1861 case(Elements):
1862 out=mesh->Elements->Id;
1863 break;
1864 case(ReducedElements):
1865 out=mesh->Elements->Id;
1866 break;
1867 case(FaceElements):
1868 out=mesh->FaceElements->Id;
1869 break;
1870 case(ReducedFaceElements):
1871 out=mesh->FaceElements->Id;
1872 break;
1873 case(Points):
1874 out=mesh->Points->Id;
1875 break;
1876 case(ContactElementsZero):
1877 out=mesh->ContactElements->Id;
1878 break;
1879 case(ReducedContactElementsZero):
1880 out=mesh->ContactElements->Id;
1881 break;
1882 case(ContactElementsOne):
1883 out=mesh->ContactElements->Id;
1884 break;
1885 case(ReducedContactElementsOne):
1886 out=mesh->ContactElements->Id;
1887 break;
1888 case(DegreesOfFreedom):
1889 out=mesh->Nodes->degreesOfFreedomId;
1890 break;
1891 case(ReducedDegreesOfFreedom):
1892 out=mesh->Nodes->reducedDegreesOfFreedomId;
1893 break;
1894 default:
1895 stringstream temp;
1896 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1897 throw FinleyAdapterException(temp.str());
1898 break;
1899 }
1900 return out;
1901 }
1902 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1903 {
1904 int out=0;
1905 Finley_Mesh* mesh=m_finleyMesh.get();
1906 switch (functionSpaceType) {
1907 case(Nodes):
1908 out=mesh->Nodes->Tag[sampleNo];
1909 break;
1910 case(ReducedNodes):
1911 throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1912 break;
1913 case(Elements):
1914 out=mesh->Elements->Tag[sampleNo];
1915 break;
1916 case(ReducedElements):
1917 out=mesh->Elements->Tag[sampleNo];
1918 break;
1919 case(FaceElements):
1920 out=mesh->FaceElements->Tag[sampleNo];
1921 break;
1922 case(ReducedFaceElements):
1923 out=mesh->FaceElements->Tag[sampleNo];
1924 break;
1925 case(Points):
1926 out=mesh->Points->Tag[sampleNo];
1927 break;
1928 case(ContactElementsZero):
1929 out=mesh->ContactElements->Tag[sampleNo];
1930 break;
1931 case(ReducedContactElementsZero):
1932 out=mesh->ContactElements->Tag[sampleNo];
1933 break;
1934 case(ContactElementsOne):
1935 out=mesh->ContactElements->Tag[sampleNo];
1936 break;
1937 case(ReducedContactElementsOne):
1938 out=mesh->ContactElements->Tag[sampleNo];
1939 break;
1940 case(DegreesOfFreedom):
1941 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1942 break;
1943 case(ReducedDegreesOfFreedom):
1944 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1945 break;
1946 default:
1947 stringstream temp;
1948 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1949 throw FinleyAdapterException(temp.str());
1950 break;
1951 }
1952 return out;
1953 }
1954
1955
1956 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1957 {
1958 Finley_Mesh* mesh=m_finleyMesh.get();
1959 escriptDataC tmp=mask.getDataC();
1960 switch(functionSpaceType) {
1961 case(Nodes):
1962 Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1963 break;
1964 case(ReducedNodes):
1965 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
1966 break;
1967 case(DegreesOfFreedom):
1968 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1969 break;
1970 case(ReducedDegreesOfFreedom):
1971 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1972 break;
1973 case(Elements):
1974 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1975 break;
1976 case(ReducedElements):
1977 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1978 break;
1979 case(FaceElements):
1980 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1981 break;
1982 case(ReducedFaceElements):
1983 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1984 break;
1985 case(Points):
1986 Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1987 break;
1988 case(ContactElementsZero):
1989 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1990 break;
1991 case(ReducedContactElementsZero):
1992 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1993 break;
1994 case(ContactElementsOne):
1995 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1996 break;
1997 case(ReducedContactElementsOne):
1998 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1999 break;
2000 default:
2001 stringstream temp;
2002 temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2003 throw FinleyAdapterException(temp.str());
2004 }
2005 checkFinleyError();
2006 return;
2007 }
2008
2009 void MeshAdapter::setTagMap(const std::string& name, int tag)
2010 {
2011 Finley_Mesh* mesh=m_finleyMesh.get();
2012 Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2013 checkPasoError();
2014 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2015 }
2016
2017 int MeshAdapter::getTag(const std::string& name) const
2018 {
2019 Finley_Mesh* mesh=m_finleyMesh.get();
2020 int tag=0;
2021 tag=Finley_Mesh_getTag(mesh, name.c_str());
2022 checkPasoError();
2023 // throwStandardException("MeshAdapter::getTag is not implemented.");
2024 return tag;
2025 }
2026
2027 bool MeshAdapter::isValidTagName(const std::string& name) const
2028 {
2029 Finley_Mesh* mesh=m_finleyMesh.get();
2030 return Finley_Mesh_isValidTagName(mesh,name.c_str());
2031 }
2032
2033 std::string MeshAdapter::showTagNames() const
2034 {
2035 stringstream temp;
2036 Finley_Mesh* mesh=m_finleyMesh.get();
2037 Finley_TagMap* tag_map=mesh->TagMap;
2038 while (tag_map) {
2039 temp << tag_map->name;
2040 tag_map=tag_map->next;
2041 if (tag_map) temp << ", ";
2042 }
2043 return temp.str();
2044 }
2045
2046 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2047 {
2048 Finley_Mesh* mesh=m_finleyMesh.get();
2049 dim_t numTags=0;
2050 switch(functionSpaceCode) {
2051 case(Nodes):
2052 numTags=mesh->Nodes->numTagsInUse;
2053 break;
2054 case(ReducedNodes):
2055 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2056 break;
2057 case(DegreesOfFreedom):
2058 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2059 break;
2060 case(ReducedDegreesOfFreedom):
2061 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2062 break;
2063 case(Elements):
2064 case(ReducedElements):
2065 numTags=mesh->Elements->numTagsInUse;
2066 break;
2067 case(FaceElements):
2068 case(ReducedFaceElements):
2069 numTags=mesh->FaceElements->numTagsInUse;
2070 break;
2071 case(Points):
2072 numTags=mesh->Points->numTagsInUse;
2073 break;
2074 case(ContactElementsZero):
2075 case(ReducedContactElementsZero):
2076 case(ContactElementsOne):
2077 case(ReducedContactElementsOne):
2078 numTags=mesh->ContactElements->numTagsInUse;
2079 break;
2080 default:
2081 stringstream temp;
2082 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2083 throw FinleyAdapterException(temp.str());
2084 }
2085 return numTags;
2086 }
2087 int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2088 {
2089 Finley_Mesh* mesh=m_finleyMesh.get();
2090 index_t* tags=NULL;
2091 switch(functionSpaceCode) {
2092 case(Nodes):
2093 tags=mesh->Nodes->tagsInUse;
2094 break;
2095 case(ReducedNodes):
2096 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2097 break;
2098 case(DegreesOfFreedom):
2099 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2100 break;
2101 case(ReducedDegreesOfFreedom):
2102 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2103 break;
2104 case(Elements):
2105 case(ReducedElements):
2106 tags=mesh->Elements->tagsInUse;
2107 break;
2108 case(FaceElements):
2109 case(ReducedFaceElements):
2110 tags=mesh->FaceElements->tagsInUse;
2111 break;
2112 case(Points):
2113 tags=mesh->Points->tagsInUse;
2114 break;
2115 case(ContactElementsZero):
2116 case(ReducedContactElementsZero):
2117 case(ContactElementsOne):
2118 case(ReducedContactElementsOne):
2119 tags=mesh->ContactElements->tagsInUse;
2120 break;
2121 default:
2122 stringstream temp;
2123 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2124 throw FinleyAdapterException(temp.str());
2125 }
2126 return tags;
2127 }
2128
2129
2130 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26