/[escript]/branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Contents of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2078 - (show annotations)
Thu Nov 20 16:10:10 2008 UTC (10 years, 5 months ago) by phornby
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 85322 byte(s)
Two changes.

1. Move blocktimer from escript to esysUtils.
2. Make it possible to link to paso as a DLL or .so.

Should have no effect on 'nix's

In respect of 1., blocktimer had begun to spring up everywhere, so
for the moment I thought it best to move it to the only other library that
pops up all over the place.

In respect of 2., paso needed to be a DLL in order to use the windows intelc /fast
option, which does aggressive multi-file optimisations. Even in its current form, it either
vectorises or parallelises  hundreds more loops in the esys system than appear in the pragmas.

In achieving 2. I have not been too delicate in adding

PASO_DLL_API

declarations to the .h files in paso/src. Only toward the end of the process of
the conversion, when the number of linker errors dropped below 20, say, did I choosy about what
functions in a header I declared PASO_DLL_API. As a result, there are likely to be many routines
declared as external functions symbols that are in fact internal to the paso DLL. 
Why is this an issue?  It prevents the intelc compiler from getting aggressive on the paso module.
With pain there is sometimes gain. At least all the DLL rules in windows give good
(non-microsoft) compiler writers a chance to really shine.

So, if you should see a PASO_DLL_API on a function in a paso header file,
and think to yourself, "that function is only called in paso, why export it?", then feel free to
delete the PASO_DLL_API export declaration.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26