/[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 2641 - (show annotations)
Mon Aug 31 07:41:49 2009 UTC (9 years, 6 months ago) by jfenwick
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 88704 byte(s)
Fixed some of my stupids related to MPI compile errors.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26