/[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 2078 - (hide annotations)
Thu Nov 20 16:10:10 2008 UTC (10 years, 4 months ago) by phornby
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 85322 byte(s)
Two changes.

1. Move blocktimer from escript to esysUtils.
2. Make it possible to link to paso as a DLL or .so.

Should have no effect on 'nix's

In respect of 1., blocktimer had begun to spring up everywhere, so
for the moment I thought it best to move it to the only other library that
pops up all over the place.

In respect of 2., paso needed to be a DLL in order to use the windows intelc /fast
option, which does aggressive multi-file optimisations. Even in its current form, it either
vectorises or parallelises  hundreds more loops in the esys system than appear in the pragmas.

In achieving 2. I have not been too delicate in adding

PASO_DLL_API

declarations to the .h files in paso/src. Only toward the end of the process of
the conversion, when the number of linker errors dropped below 20, say, did I choosy about what
functions in a header I declared PASO_DLL_API. As a result, there are likely to be many routines
declared as external functions symbols that are in fact internal to the paso DLL. 
Why is this an issue?  It prevents the intelc compiler from getting aggressive on the paso module.
With pain there is sometimes gain. At least all the DLL rules in windows give good
(non-microsoft) compiler writers a chance to really shine.

So, if you should see a PASO_DLL_API on a function in a paso header file,
and think to yourself, "that function is only called in paso, why export it?", then feel free to
delete the PASO_DLL_API export declaration.

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