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

Contents of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26