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

Annotation of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1802 - (hide annotations)
Tue Sep 23 01:03:29 2008 UTC (10 years, 8 months ago) by jfenwick
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 84417 byte(s)
Added canTag methods to FunctionSpace and AbstractDomain (and its 
offspring).
This checks to see if the domain supports tags for the given type of 
function space.

Constructors for DataTagged now throw exceptions if you attempt to make 
a DataTagged with a FunctionSpace which does not support tags.
To allow the default constructor to work, NullDomain has a single 
functioncode which "supports" tagging.

Fixed a bug in DataTagged::toString and DataTypes::pointToString.

Added FunctionSpace::getListOfTagsSTL.

algorithm(DataTagged, BinaryFunction) in DataAlgorithm now only 
processes tags known to be in use.
This fixes mantis issue #0000186.

Added comment to Data.h intro warning about holding references if the 
underlying DataAbstract changes.

_python_ unit tests have been updated to test TaggedData with invalid 
FunctionSpaces and to give the correct answers to Lsup etc.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26