/[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 1755 - (show annotations)
Mon Sep 8 01:34:40 2008 UTC (10 years, 7 months ago) by ksteube
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 85272 byte(s)
Added new field to NetCDF dump of mesh

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26