/[escript]/branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Contents of /branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2635 - (show annotations)
Thu Aug 27 04:54:41 2009 UTC (9 years, 7 months ago) by jfenwick
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 88534 byte(s)
A bunch of changes related to saveDataCSV.
[Not completed or unit tested yet]

Added saveDataCSV to util.py
AbstractDomain (and MeshAdapter) have a commonFunctionSpace method to 
take a group of FunctionSpaces and return something they can all be interpolated to.

Added pointToStream() in DataTypes to help print points.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26