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

Annotation of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2641 - (hide annotations)
Mon Aug 31 07:41:49 2009 UTC (9 years, 9 months ago) by jfenwick
File size: 88704 byte(s)
Fixed some of my stupids related to MPI compile errors.

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