/[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 2151 - (show annotations)
Wed Dec 10 06:34:25 2008 UTC (10 years, 8 months ago) by caltinay
File size: 84042 byte(s)
Moved common functionality from saveDX/VTK into private method. May be used
for future save* methods.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26