/[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 1872 - (show annotations)
Mon Oct 13 00:18:55 2008 UTC (10 years, 6 months ago) by jfenwick
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 85032 byte(s)
Closing the moreshared branch

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26