/[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 2748 - (show annotations)
Tue Nov 17 07:32:59 2009 UTC (9 years, 5 months ago) by gross
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 89901 byte(s)
Macro elements are implemented now. VTK writer for macro elements still needs testing.
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->referenceElementSet->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->referenceElementSet->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->referenceElementSet->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->referenceElementSet->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->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
697 }
698 break;
699 case(ReducedElements):
700 if (mesh->Elements!=NULL) {
701 numSamples=mesh->Elements->numElements;
702 numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
703 }
704 break;
705 case(FaceElements):
706 if (mesh->FaceElements!=NULL) {
707 numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
708 numSamples=mesh->FaceElements->numElements;
709 }
710 break;
711 case(ReducedFaceElements):
712 if (mesh->FaceElements!=NULL) {
713 numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->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->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
726 numSamples=mesh->ContactElements->numElements;
727 }
728 break;
729 case(ReducedContactElementsZero):
730 if (mesh->ContactElements!=NULL) {
731 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
732 numSamples=mesh->ContactElements->numElements;
733 }
734 break;
735 case(ContactElementsOne):
736 if (mesh->ContactElements!=NULL) {
737 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
738 numSamples=mesh->ContactElements->numElements;
739 }
740 break;
741 case(ReducedContactElementsOne):
742 if (mesh->ContactElements!=NULL) {
743 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->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 #ifdef PASO_MPI
1560 index_t myFirstNode=0, myLastNode=0, k=0;
1561 index_t* globalNodeIndex=0;
1562 Finley_Mesh* mesh_p=m_finleyMesh.get();
1563 if (fs_code == FINLEY_REDUCED_NODES)
1564 {
1565 myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1566 myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1567 globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1568 }
1569 else
1570 {
1571 myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1572 myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1573 globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1574 }
1575 k=globalNodeIndex[id];
1576 return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1577 #endif
1578 return true;
1579 }
1580
1581
1582
1583 //
1584 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1585 //
1586 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1587 const int row_blocksize,
1588 const escript::FunctionSpace& row_functionspace,
1589 const int column_blocksize,
1590 const escript::FunctionSpace& column_functionspace,
1591 const int type) const
1592 {
1593 int reduceRowOrder=0;
1594 int reduceColOrder=0;
1595 // is the domain right?
1596 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1597 if (row_domain!=*this)
1598 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1599 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1600 if (col_domain!=*this)
1601 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1602 // is the function space type right
1603 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1604 reduceRowOrder=0;
1605 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1606 reduceRowOrder=1;
1607 } else {
1608 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1609 }
1610 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1611 reduceColOrder=0;
1612 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1613 reduceColOrder=1;
1614 } else {
1615 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1616 }
1617 // generate matrix:
1618
1619 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1620 checkFinleyError();
1621 Paso_SystemMatrix* fsystemMatrix;
1622 int trilinos = 0;
1623 if (trilinos) {
1624 #ifdef TRILINOS
1625 /* Allocation Epetra_VrbMatrix here */
1626 #endif
1627 }
1628 else {
1629 fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1630 }
1631 checkPasoError();
1632 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1633 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1634 }
1635
1636 //
1637 // creates a TransportProblemAdapter
1638 //
1639 TransportProblemAdapter MeshAdapter::newTransportProblem(
1640 const double theta,
1641 const int blocksize,
1642 const escript::FunctionSpace& functionspace,
1643 const int type) const
1644 {
1645 int reduceOrder=0;
1646 // is the domain right?
1647 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1648 if (domain!=*this)
1649 throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1650 // is the function space type right
1651 if (functionspace.getTypeCode()==DegreesOfFreedom) {
1652 reduceOrder=0;
1653 } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1654 reduceOrder=1;
1655 } else {
1656 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1657 }
1658 // generate matrix:
1659
1660 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1661 checkFinleyError();
1662 Paso_FCTransportProblem* transportProblem;
1663 transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);
1664 checkPasoError();
1665 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1666 return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);
1667 }
1668
1669 //
1670 // vtkObject MeshAdapter::createVtkObject() const
1671 // TODO:
1672 //
1673
1674 //
1675 // returns true if data at the atom_type is considered as being cell centered:
1676 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1677 {
1678 switch(functionSpaceCode) {
1679 case(Nodes):
1680 case(DegreesOfFreedom):
1681 case(ReducedDegreesOfFreedom):
1682 return false;
1683 break;
1684 case(Elements):
1685 case(FaceElements):
1686 case(Points):
1687 case(ContactElementsZero):
1688 case(ContactElementsOne):
1689 case(ReducedElements):
1690 case(ReducedFaceElements):
1691 case(ReducedContactElementsZero):
1692 case(ReducedContactElementsOne):
1693 return true;
1694 break;
1695 default:
1696 stringstream temp;
1697 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1698 throw FinleyAdapterException(temp.str());
1699 break;
1700 }
1701 checkFinleyError();
1702 return false;
1703 }
1704
1705 bool
1706 MeshAdapter::commonFunctionSpace(const std::vector<int>& fs, int& resultcode) const
1707 {
1708 /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1709 class 1: DOF <-> Nodes
1710 class 2: ReducedDOF <-> ReducedNodes
1711 class 3: Points
1712 class 4: Elements
1713 class 5: ReducedElements
1714 class 6: FaceElements
1715 class 7: ReducedFaceElements
1716 class 8: ContactElementZero <-> ContactElementOne
1717 class 9: ReducedContactElementZero <-> ReducedContactElementOne
1718
1719 There is also a set of lines. Interpolation is possible down a line but not between lines.
1720 class 1 and 2 belong to all lines so aren't considered.
1721 line 0: class 3
1722 line 1: class 4,5
1723 line 2: class 6,7
1724 line 3: class 8,9
1725
1726 For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1727 eg hasnodes is true if we have at least one instance of Nodes.
1728 */
1729 if (fs.size()==0)
1730 {
1731 return false;
1732 }
1733 std::vector<int> hasclass(10);
1734 std::vector<int> hasline(4);
1735 bool hasnodes=false;
1736 bool hasrednodes=false;
1737 bool hascez=false;
1738 bool hasrcez=false;
1739 for (int i=0;i<fs.size();++i)
1740 {
1741 switch(fs[i])
1742 {
1743 case(Nodes): hasnodes=true; // no break is deliberate
1744 case(DegreesOfFreedom):
1745 hasclass[1]=1;
1746 break;
1747 case(ReducedNodes): hasrednodes=true; // no break is deliberate
1748 case(ReducedDegreesOfFreedom):
1749 hasclass[2]=1;
1750 break;
1751 case(Points):
1752 hasline[0]=1;
1753 hasclass[3]=1;
1754 break;
1755 case(Elements):
1756 hasclass[4]=1;
1757 hasline[1]=1;
1758 break;
1759 case(ReducedElements):
1760 hasclass[5]=1;
1761 hasline[1]=1;
1762 break;
1763 case(FaceElements):
1764 hasclass[6]=1;
1765 hasline[2]=1;
1766 break;
1767 case(ReducedFaceElements):
1768 hasclass[7]=1;
1769 hasline[2]=1;
1770 break;
1771 case(ContactElementsZero): hascez=true; // no break is deliberate
1772 case(ContactElementsOne):
1773 hasclass[8]=1;
1774 hasline[3]=1;
1775 break;
1776 case(ReducedContactElementsZero): hasrcez=true; // no break is deliberate
1777 case(ReducedContactElementsOne):
1778 hasclass[9]=1;
1779 hasline[3]=1;
1780 break;
1781 default:
1782 return false;
1783 }
1784 }
1785 int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1786 // fail if we have more than one leaf group
1787
1788 if (totlines>1)
1789 {
1790 return false; // there are at least two branches we can't interpolate between
1791 }
1792 else if (totlines==1)
1793 {
1794 if (hasline[0]==1) // we have points
1795 {
1796 resultcode=Points;
1797 }
1798 else if (hasline[1]==1)
1799 {
1800 if (hasclass[5]==1)
1801 {
1802 resultcode=ReducedElements;
1803 }
1804 else
1805 {
1806 resultcode=Elements;
1807 }
1808 }
1809 else if (hasline[2]==1)
1810 {
1811 if (hasclass[7]==1)
1812 {
1813 resultcode=ReducedFaceElements;
1814 }
1815 else
1816 {
1817 resultcode=FaceElements;
1818 }
1819 }
1820 else // so we must be in line3
1821 {
1822 if (hasclass[9]==1)
1823 {
1824 // need something from class 9
1825 resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1826 }
1827 else
1828 {
1829 // something from class 8
1830 resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1831 }
1832 }
1833 }
1834 else // totlines==0
1835 {
1836 if (hasclass[2]==1)
1837 {
1838 // something from class 2
1839 resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1840 }
1841 else
1842 { // something from class 1
1843 resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1844 }
1845 }
1846 return true;
1847 }
1848
1849 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1850 {
1851 switch(functionSpaceType_source) {
1852 case(Nodes):
1853 switch(functionSpaceType_target) {
1854 case(Nodes):
1855 case(ReducedNodes):
1856 case(ReducedDegreesOfFreedom):
1857 case(DegreesOfFreedom):
1858 case(Elements):
1859 case(ReducedElements):
1860 case(FaceElements):
1861 case(ReducedFaceElements):
1862 case(Points):
1863 case(ContactElementsZero):
1864 case(ReducedContactElementsZero):
1865 case(ContactElementsOne):
1866 case(ReducedContactElementsOne):
1867 return true;
1868 default:
1869 stringstream temp;
1870 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1871 throw FinleyAdapterException(temp.str());
1872 }
1873 break;
1874 case(ReducedNodes):
1875 switch(functionSpaceType_target) {
1876 case(ReducedNodes):
1877 case(ReducedDegreesOfFreedom):
1878 case(Elements):
1879 case(ReducedElements):
1880 case(FaceElements):
1881 case(ReducedFaceElements):
1882 case(Points):
1883 case(ContactElementsZero):
1884 case(ReducedContactElementsZero):
1885 case(ContactElementsOne):
1886 case(ReducedContactElementsOne):
1887 return true;
1888 case(Nodes):
1889 case(DegreesOfFreedom):
1890 return false;
1891 default:
1892 stringstream temp;
1893 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1894 throw FinleyAdapterException(temp.str());
1895 }
1896 break;
1897 case(Elements):
1898 if (functionSpaceType_target==Elements) {
1899 return true;
1900 } else if (functionSpaceType_target==ReducedElements) {
1901 return true;
1902 } else {
1903 return false;
1904 }
1905 case(ReducedElements):
1906 if (functionSpaceType_target==ReducedElements) {
1907 return true;
1908 } else {
1909 return false;
1910 }
1911 case(FaceElements):
1912 if (functionSpaceType_target==FaceElements) {
1913 return true;
1914 } else if (functionSpaceType_target==ReducedFaceElements) {
1915 return true;
1916 } else {
1917 return false;
1918 }
1919 case(ReducedFaceElements):
1920 if (functionSpaceType_target==ReducedFaceElements) {
1921 return true;
1922 } else {
1923 return false;
1924 }
1925 case(Points):
1926 if (functionSpaceType_target==Points) {
1927 return true;
1928 } else {
1929 return false;
1930 }
1931 case(ContactElementsZero):
1932 case(ContactElementsOne):
1933 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1934 return true;
1935 } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1936 return true;
1937 } else {
1938 return false;
1939 }
1940 case(ReducedContactElementsZero):
1941 case(ReducedContactElementsOne):
1942 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1943 return true;
1944 } else {
1945 return false;
1946 }
1947 case(DegreesOfFreedom):
1948 switch(functionSpaceType_target) {
1949 case(ReducedDegreesOfFreedom):
1950 case(DegreesOfFreedom):
1951 case(Nodes):
1952 case(ReducedNodes):
1953 case(Elements):
1954 case(ReducedElements):
1955 case(Points):
1956 case(FaceElements):
1957 case(ReducedFaceElements):
1958 case(ContactElementsZero):
1959 case(ReducedContactElementsZero):
1960 case(ContactElementsOne):
1961 case(ReducedContactElementsOne):
1962 return true;
1963 default:
1964 stringstream temp;
1965 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1966 throw FinleyAdapterException(temp.str());
1967 }
1968 break;
1969 case(ReducedDegreesOfFreedom):
1970 switch(functionSpaceType_target) {
1971 case(ReducedDegreesOfFreedom):
1972 case(ReducedNodes):
1973 case(Elements):
1974 case(ReducedElements):
1975 case(FaceElements):
1976 case(ReducedFaceElements):
1977 case(Points):
1978 case(ContactElementsZero):
1979 case(ReducedContactElementsZero):
1980 case(ContactElementsOne):
1981 case(ReducedContactElementsOne):
1982 return true;
1983 case(Nodes):
1984 case(DegreesOfFreedom):
1985 return false;
1986 default:
1987 stringstream temp;
1988 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1989 throw FinleyAdapterException(temp.str());
1990 }
1991 break;
1992 default:
1993 stringstream temp;
1994 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1995 throw FinleyAdapterException(temp.str());
1996 break;
1997 }
1998 checkFinleyError();
1999 return false;
2000 }
2001
2002 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
2003 {
2004 return false;
2005 }
2006
2007 bool MeshAdapter::operator==(const AbstractDomain& other) const
2008 {
2009 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2010 if (temp!=0) {
2011 return (m_finleyMesh==temp->m_finleyMesh);
2012 } else {
2013 return false;
2014 }
2015 }
2016
2017 bool MeshAdapter::operator!=(const AbstractDomain& other) const
2018 {
2019 return !(operator==(other));
2020 }
2021
2022 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2023 {
2024 Finley_Mesh* mesh=m_finleyMesh.get();
2025 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2026 checkPasoError();
2027 return out;
2028 }
2029 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2030 {
2031 Finley_Mesh* mesh=m_finleyMesh.get();
2032 int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2033 checkPasoError();
2034 return out;
2035 }
2036
2037 escript::Data MeshAdapter::getX() const
2038 {
2039 return continuousFunction(asAbstractContinuousDomain()).getX();
2040 }
2041
2042 escript::Data MeshAdapter::getNormal() const
2043 {
2044 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
2045 }
2046
2047 escript::Data MeshAdapter::getSize() const
2048 {
2049 return escript::function(asAbstractContinuousDomain()).getSize();
2050 }
2051
2052 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2053 {
2054 int *out = NULL;
2055 Finley_Mesh* mesh=m_finleyMesh.get();
2056 switch (functionSpaceType) {
2057 case(Nodes):
2058 out=mesh->Nodes->Id;
2059 break;
2060 case(ReducedNodes):
2061 out=mesh->Nodes->reducedNodesId;
2062 break;
2063 case(Elements):
2064 out=mesh->Elements->Id;
2065 break;
2066 case(ReducedElements):
2067 out=mesh->Elements->Id;
2068 break;
2069 case(FaceElements):
2070 out=mesh->FaceElements->Id;
2071 break;
2072 case(ReducedFaceElements):
2073 out=mesh->FaceElements->Id;
2074 break;
2075 case(Points):
2076 out=mesh->Points->Id;
2077 break;
2078 case(ContactElementsZero):
2079 out=mesh->ContactElements->Id;
2080 break;
2081 case(ReducedContactElementsZero):
2082 out=mesh->ContactElements->Id;
2083 break;
2084 case(ContactElementsOne):
2085 out=mesh->ContactElements->Id;
2086 break;
2087 case(ReducedContactElementsOne):
2088 out=mesh->ContactElements->Id;
2089 break;
2090 case(DegreesOfFreedom):
2091 out=mesh->Nodes->degreesOfFreedomId;
2092 break;
2093 case(ReducedDegreesOfFreedom):
2094 out=mesh->Nodes->reducedDegreesOfFreedomId;
2095 break;
2096 default:
2097 stringstream temp;
2098 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2099 throw FinleyAdapterException(temp.str());
2100 break;
2101 }
2102 return out;
2103 }
2104 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2105 {
2106 int out=0;
2107 Finley_Mesh* mesh=m_finleyMesh.get();
2108 switch (functionSpaceType) {
2109 case(Nodes):
2110 out=mesh->Nodes->Tag[sampleNo];
2111 break;
2112 case(ReducedNodes):
2113 throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2114 break;
2115 case(Elements):
2116 out=mesh->Elements->Tag[sampleNo];
2117 break;
2118 case(ReducedElements):
2119 out=mesh->Elements->Tag[sampleNo];
2120 break;
2121 case(FaceElements):
2122 out=mesh->FaceElements->Tag[sampleNo];
2123 break;
2124 case(ReducedFaceElements):
2125 out=mesh->FaceElements->Tag[sampleNo];
2126 break;
2127 case(Points):
2128 out=mesh->Points->Tag[sampleNo];
2129 break;
2130 case(ContactElementsZero):
2131 out=mesh->ContactElements->Tag[sampleNo];
2132 break;
2133 case(ReducedContactElementsZero):
2134 out=mesh->ContactElements->Tag[sampleNo];
2135 break;
2136 case(ContactElementsOne):
2137 out=mesh->ContactElements->Tag[sampleNo];
2138 break;
2139 case(ReducedContactElementsOne):
2140 out=mesh->ContactElements->Tag[sampleNo];
2141 break;
2142 case(DegreesOfFreedom):
2143 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2144 break;
2145 case(ReducedDegreesOfFreedom):
2146 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2147 break;
2148 default:
2149 stringstream temp;
2150 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2151 throw FinleyAdapterException(temp.str());
2152 break;
2153 }
2154 return out;
2155 }
2156
2157
2158 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2159 {
2160 Finley_Mesh* mesh=m_finleyMesh.get();
2161 escriptDataC tmp=mask.getDataC();
2162 switch(functionSpaceType) {
2163 case(Nodes):
2164 Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2165 break;
2166 case(ReducedNodes):
2167 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2168 break;
2169 case(DegreesOfFreedom):
2170 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2171 break;
2172 case(ReducedDegreesOfFreedom):
2173 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2174 break;
2175 case(Elements):
2176 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2177 break;
2178 case(ReducedElements):
2179 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2180 break;
2181 case(FaceElements):
2182 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2183 break;
2184 case(ReducedFaceElements):
2185 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2186 break;
2187 case(Points):
2188 Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2189 break;
2190 case(ContactElementsZero):
2191 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2192 break;
2193 case(ReducedContactElementsZero):
2194 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2195 break;
2196 case(ContactElementsOne):
2197 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2198 break;
2199 case(ReducedContactElementsOne):
2200 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2201 break;
2202 default:
2203 stringstream temp;
2204 temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2205 throw FinleyAdapterException(temp.str());
2206 }
2207 checkFinleyError();
2208 return;
2209 }
2210
2211 void MeshAdapter::setTagMap(const std::string& name, int tag)
2212 {
2213 Finley_Mesh* mesh=m_finleyMesh.get();
2214 Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2215 checkPasoError();
2216 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2217 }
2218
2219 int MeshAdapter::getTag(const std::string& name) const
2220 {
2221 Finley_Mesh* mesh=m_finleyMesh.get();
2222 int tag=0;
2223 tag=Finley_Mesh_getTag(mesh, name.c_str());
2224 checkPasoError();
2225 // throwStandardException("MeshAdapter::getTag is not implemented.");
2226 return tag;
2227 }
2228
2229 bool MeshAdapter::isValidTagName(const std::string& name) const
2230 {
2231 Finley_Mesh* mesh=m_finleyMesh.get();
2232 return Finley_Mesh_isValidTagName(mesh,name.c_str());
2233 }
2234
2235 std::string MeshAdapter::showTagNames() const
2236 {
2237 stringstream temp;
2238 Finley_Mesh* mesh=m_finleyMesh.get();
2239 Finley_TagMap* tag_map=mesh->TagMap;
2240 while (tag_map) {
2241 temp << tag_map->name;
2242 tag_map=tag_map->next;
2243 if (tag_map) temp << ", ";
2244 }
2245 return temp.str();
2246 }
2247
2248 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2249 {
2250 Finley_Mesh* mesh=m_finleyMesh.get();
2251 dim_t numTags=0;
2252 switch(functionSpaceCode) {
2253 case(Nodes):
2254 numTags=mesh->Nodes->numTagsInUse;
2255 break;
2256 case(ReducedNodes):
2257 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2258 break;
2259 case(DegreesOfFreedom):
2260 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2261 break;
2262 case(ReducedDegreesOfFreedom):
2263 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2264 break;
2265 case(Elements):
2266 case(ReducedElements):
2267 numTags=mesh->Elements->numTagsInUse;
2268 break;
2269 case(FaceElements):
2270 case(ReducedFaceElements):
2271 numTags=mesh->FaceElements->numTagsInUse;
2272 break;
2273 case(Points):
2274 numTags=mesh->Points->numTagsInUse;
2275 break;
2276 case(ContactElementsZero):
2277 case(ReducedContactElementsZero):
2278 case(ContactElementsOne):
2279 case(ReducedContactElementsOne):
2280 numTags=mesh->ContactElements->numTagsInUse;
2281 break;
2282 default:
2283 stringstream temp;
2284 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2285 throw FinleyAdapterException(temp.str());
2286 }
2287 return numTags;
2288 }
2289
2290 const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2291 {
2292 Finley_Mesh* mesh=m_finleyMesh.get();
2293 index_t* tags=NULL;
2294 switch(functionSpaceCode) {
2295 case(Nodes):
2296 tags=mesh->Nodes->tagsInUse;
2297 break;
2298 case(ReducedNodes):
2299 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2300 break;
2301 case(DegreesOfFreedom):
2302 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2303 break;
2304 case(ReducedDegreesOfFreedom):
2305 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2306 break;
2307 case(Elements):
2308 case(ReducedElements):
2309 tags=mesh->Elements->tagsInUse;
2310 break;
2311 case(FaceElements):
2312 case(ReducedFaceElements):
2313 tags=mesh->FaceElements->tagsInUse;
2314 break;
2315 case(Points):
2316 tags=mesh->Points->tagsInUse;
2317 break;
2318 case(ContactElementsZero):
2319 case(ReducedContactElementsZero):
2320 case(ContactElementsOne):
2321 case(ReducedContactElementsOne):
2322 tags=mesh->ContactElements->tagsInUse;
2323 break;
2324 default:
2325 stringstream temp;
2326 temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2327 throw FinleyAdapterException(temp.str());
2328 }
2329 return tags;
2330 }
2331
2332
2333 bool MeshAdapter::canTag(int functionSpaceCode) const
2334 {
2335 switch(functionSpaceCode) {
2336 case(Nodes):
2337 case(Elements):
2338 case(ReducedElements):
2339 case(FaceElements):
2340 case(ReducedFaceElements):
2341 case(Points):
2342 case(ContactElementsZero):
2343 case(ReducedContactElementsZero):
2344 case(ContactElementsOne):
2345 case(ReducedContactElementsOne):
2346 return true;
2347 case(ReducedNodes):
2348 case(DegreesOfFreedom):
2349 case(ReducedDegreesOfFreedom):
2350 return false;
2351 default:
2352 return false;
2353 }
2354 }
2355
2356 AbstractDomain::StatusType MeshAdapter::getStatus() const
2357 {
2358 Finley_Mesh* mesh=m_finleyMesh.get();
2359 return Finley_Mesh_getStatus(mesh);
2360 }
2361
2362
2363
2364 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26