/[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 1800 - (show annotations)
Thu Sep 18 05:28:20 2008 UTC (11 years ago) by ksteube
File size: 83999 byte(s)
Serialized parallel I/O when writing mesh or data to NetCDF file on multiple MPI processors.
Added domain method getMPIComm() to complement getMPISize() and getMPIRank().

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26