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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26