/[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 1801 - (show annotations)
Fri Sep 19 01:37:09 2008 UTC (10 years, 9 months ago) by ksteube
File size: 83909 byte(s)
Fixed serialization of I/O for MPI...code didn't compile without MPI

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26