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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2748 - (hide annotations)
Tue Nov 17 07:32:59 2009 UTC (9 years, 6 months ago) by gross
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 89901 byte(s)
Macro elements are implemented now. VTK writer for macro elements still needs testing.
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 gross 2748 if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
244 ksteube 1347 throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);
245 gross 2748 if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
246 ksteube 1347 throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);
247 gross 2748 if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
248 ksteube 1347 throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);
249 gross 2748 if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
250 ksteube 1347 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 gross 2748 numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
697 phornby 1628 }
698     break;
699     case(ReducedElements):
700     if (mesh->Elements!=NULL) {
701     numSamples=mesh->Elements->numElements;
702 gross 2748 numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
703 phornby 1628 }
704     break;
705     case(FaceElements):
706     if (mesh->FaceElements!=NULL) {
707 gross 2748 numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
708 phornby 1628 numSamples=mesh->FaceElements->numElements;
709     }
710     break;
711     case(ReducedFaceElements):
712     if (mesh->FaceElements!=NULL) {
713 gross 2748 numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
714 phornby 1628 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 gross 2748 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
726 phornby 1628 numSamples=mesh->ContactElements->numElements;
727     }
728     break;
729     case(ReducedContactElementsZero):
730     if (mesh->ContactElements!=NULL) {
731 gross 2748 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
732 phornby 1628 numSamples=mesh->ContactElements->numElements;
733     }
734     break;
735     case(ContactElementsOne):
736     if (mesh->ContactElements!=NULL) {
737 gross 2748 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
738 phornby 1628 numSamples=mesh->ContactElements->numElements;
739     }
740     break;
741     case(ReducedContactElementsOne):
742     if (mesh->ContactElements!=NULL) {
743 gross 2748 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
744 phornby 1628 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 jfenwick 2642 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1558     {
1559 jfenwick 2644 #ifdef PASO_MPI
1560 jfenwick 2642 index_t myFirstNode=0, myLastNode=0, k=0;
1561     index_t* globalNodeIndex=0;
1562     Finley_Mesh* mesh_p=m_finleyMesh.get();
1563     if (fs_code == FINLEY_REDUCED_NODES)
1564     {
1565     myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1566     myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1567     globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1568     }
1569     else
1570     {
1571     myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1572     myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1573     globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1574     }
1575     k=globalNodeIndex[id];
1576     return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1577 jfenwick 2644 #endif
1578     return true;
1579 jfenwick 2642 }
1580    
1581    
1582    
1583 caltinay 2151 //
1584     // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1585     //
1586 jgs 82 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1587 phornby 1628 const int row_blocksize,
1588     const escript::FunctionSpace& row_functionspace,
1589     const int column_blocksize,
1590     const escript::FunctionSpace& column_functionspace,
1591     const int type) const
1592 jgs 82 {
1593 phornby 1628 int reduceRowOrder=0;
1594     int reduceColOrder=0;
1595     // is the domain right?
1596 jfenwick 1872 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1597 phornby 1628 if (row_domain!=*this)
1598     throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1599 jfenwick 1872 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1600 phornby 1628 if (col_domain!=*this)
1601     throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1602     // is the function space type right
1603     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1604     reduceRowOrder=0;
1605     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1606     reduceRowOrder=1;
1607     } else {
1608     throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1609     }
1610     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1611     reduceColOrder=0;
1612     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1613     reduceColOrder=1;
1614     } else {
1615     throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1616     }
1617     // generate matrix:
1618    
1619     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1620     checkFinleyError();
1621     Paso_SystemMatrix* fsystemMatrix;
1622     int trilinos = 0;
1623     if (trilinos) {
1624 ksteube 1312 #ifdef TRILINOS
1625     /* Allocation Epetra_VrbMatrix here */
1626     #endif
1627 phornby 1628 }
1628     else {
1629 gross 2551 fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1630 phornby 1628 }
1631     checkPasoError();
1632     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1633     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1634 jgs 82 }
1635 caltinay 2151
1636     //
1637 gross 1366 // creates a TransportProblemAdapter
1638 caltinay 2151 //
1639 gross 1366 TransportProblemAdapter MeshAdapter::newTransportProblem(
1640 phornby 1628 const double theta,
1641     const int blocksize,
1642     const escript::FunctionSpace& functionspace,
1643     const int type) const
1644 gross 1366 {
1645 phornby 1628 int reduceOrder=0;
1646     // is the domain right?
1647 jfenwick 1872 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1648 phornby 1628 if (domain!=*this)
1649     throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1650     // is the function space type right
1651     if (functionspace.getTypeCode()==DegreesOfFreedom) {
1652     reduceOrder=0;
1653     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1654     reduceOrder=1;
1655     } else {
1656     throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1657     }
1658     // generate matrix:
1659    
1660     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1661     checkFinleyError();
1662     Paso_FCTransportProblem* transportProblem;
1663     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);
1664     checkPasoError();
1665     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1666     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);
1667 gross 1366 }
1668 jgs 149
1669 jgs 82 //
1670     // vtkObject MeshAdapter::createVtkObject() const
1671     // TODO:
1672     //
1673 jgs 149
1674 jgs 82 //
1675     // returns true if data at the atom_type is considered as being cell centered:
1676     bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1677     {
1678 phornby 1628 switch(functionSpaceCode) {
1679     case(Nodes):
1680     case(DegreesOfFreedom):
1681     case(ReducedDegreesOfFreedom):
1682     return false;
1683     break;
1684     case(Elements):
1685     case(FaceElements):
1686     case(Points):
1687     case(ContactElementsZero):
1688     case(ContactElementsOne):
1689     case(ReducedElements):
1690     case(ReducedFaceElements):
1691     case(ReducedContactElementsZero):
1692     case(ReducedContactElementsOne):
1693     return true;
1694     break;
1695     default:
1696     stringstream temp;
1697     temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1698     throw FinleyAdapterException(temp.str());
1699     break;
1700     }
1701     checkFinleyError();
1702     return false;
1703 jgs 82 }
1704 jgs 149
1705 jfenwick 2635 bool
1706     MeshAdapter::commonFunctionSpace(const std::vector<int>& fs, int& resultcode) const
1707     {
1708     /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1709     class 1: DOF <-> Nodes
1710     class 2: ReducedDOF <-> ReducedNodes
1711     class 3: Points
1712     class 4: Elements
1713     class 5: ReducedElements
1714     class 6: FaceElements
1715     class 7: ReducedFaceElements
1716     class 8: ContactElementZero <-> ContactElementOne
1717     class 9: ReducedContactElementZero <-> ReducedContactElementOne
1718    
1719     There is also a set of lines. Interpolation is possible down a line but not between lines.
1720     class 1 and 2 belong to all lines so aren't considered.
1721     line 0: class 3
1722     line 1: class 4,5
1723     line 2: class 6,7
1724     line 3: class 8,9
1725    
1726     For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1727     eg hasnodes is true if we have at least one instance of Nodes.
1728     */
1729     if (fs.size()==0)
1730     {
1731     return false;
1732     }
1733     std::vector<int> hasclass(10);
1734     std::vector<int> hasline(4);
1735     bool hasnodes=false;
1736     bool hasrednodes=false;
1737     bool hascez=false;
1738     bool hasrcez=false;
1739     for (int i=0;i<fs.size();++i)
1740     {
1741     switch(fs[i])
1742     {
1743     case(Nodes): hasnodes=true; // no break is deliberate
1744     case(DegreesOfFreedom):
1745     hasclass[1]=1;
1746     break;
1747     case(ReducedNodes): hasrednodes=true; // no break is deliberate
1748     case(ReducedDegreesOfFreedom):
1749     hasclass[2]=1;
1750     break;
1751     case(Points):
1752     hasline[0]=1;
1753     hasclass[3]=1;
1754     break;
1755     case(Elements):
1756     hasclass[4]=1;
1757     hasline[1]=1;
1758     break;
1759     case(ReducedElements):
1760     hasclass[5]=1;
1761     hasline[1]=1;
1762     break;
1763     case(FaceElements):
1764     hasclass[6]=1;
1765     hasline[2]=1;
1766     break;
1767     case(ReducedFaceElements):
1768     hasclass[7]=1;
1769     hasline[2]=1;
1770     break;
1771     case(ContactElementsZero): hascez=true; // no break is deliberate
1772     case(ContactElementsOne):
1773     hasclass[8]=1;
1774     hasline[3]=1;
1775     break;
1776     case(ReducedContactElementsZero): hasrcez=true; // no break is deliberate
1777     case(ReducedContactElementsOne):
1778     hasclass[9]=1;
1779     hasline[3]=1;
1780     break;
1781     default:
1782     return false;
1783     }
1784     }
1785     int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1786     // fail if we have more than one leaf group
1787    
1788     if (totlines>1)
1789     {
1790     return false; // there are at least two branches we can't interpolate between
1791     }
1792     else if (totlines==1)
1793     {
1794     if (hasline[0]==1) // we have points
1795     {
1796     resultcode=Points;
1797     }
1798     else if (hasline[1]==1)
1799     {
1800     if (hasclass[5]==1)
1801     {
1802     resultcode=ReducedElements;
1803     }
1804     else
1805     {
1806     resultcode=Elements;
1807     }
1808     }
1809     else if (hasline[2]==1)
1810     {
1811 jfenwick 2644 if (hasclass[7]==1)
1812 jfenwick 2635 {
1813     resultcode=ReducedFaceElements;
1814     }
1815     else
1816     {
1817     resultcode=FaceElements;
1818     }
1819     }
1820     else // so we must be in line3
1821     {
1822     if (hasclass[9]==1)
1823     {
1824     // need something from class 9
1825     resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1826     }
1827     else
1828     {
1829     // something from class 8
1830     resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1831     }
1832     }
1833     }
1834     else // totlines==0
1835     {
1836     if (hasclass[2]==1)
1837     {
1838     // something from class 2
1839     resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1840     }
1841     else
1842     { // something from class 1
1843     resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1844     }
1845     }
1846     return true;
1847     }
1848    
1849 jgs 82 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1850     {
1851 phornby 1628 switch(functionSpaceType_source) {
1852     case(Nodes):
1853 jfenwick 2635 switch(functionSpaceType_target) {
1854     case(Nodes):
1855     case(ReducedNodes):
1856     case(ReducedDegreesOfFreedom):
1857     case(DegreesOfFreedom):
1858     case(Elements):
1859     case(ReducedElements):
1860     case(FaceElements):
1861     case(ReducedFaceElements):
1862     case(Points):
1863     case(ContactElementsZero):
1864     case(ReducedContactElementsZero):
1865     case(ContactElementsOne):
1866     case(ReducedContactElementsOne):
1867     return true;
1868     default:
1869     stringstream temp;
1870     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1871     throw FinleyAdapterException(temp.str());
1872 phornby 1628 }
1873     break;
1874     case(ReducedNodes):
1875 jfenwick 2635 switch(functionSpaceType_target) {
1876     case(ReducedNodes):
1877     case(ReducedDegreesOfFreedom):
1878     case(Elements):
1879     case(ReducedElements):
1880     case(FaceElements):
1881     case(ReducedFaceElements):
1882     case(Points):
1883     case(ContactElementsZero):
1884     case(ReducedContactElementsZero):
1885     case(ContactElementsOne):
1886     case(ReducedContactElementsOne):
1887     return true;
1888     case(Nodes):
1889     case(DegreesOfFreedom):
1890     return false;
1891     default:
1892     stringstream temp;
1893     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1894     throw FinleyAdapterException(temp.str());
1895 phornby 1628 }
1896     break;
1897     case(Elements):
1898 jfenwick 2635 if (functionSpaceType_target==Elements) {
1899     return true;
1900     } else if (functionSpaceType_target==ReducedElements) {
1901     return true;
1902     } else {
1903     return false;
1904     }
1905 phornby 1628 case(ReducedElements):
1906 jfenwick 2635 if (functionSpaceType_target==ReducedElements) {
1907     return true;
1908     } else {
1909     return false;
1910     }
1911 phornby 1628 case(FaceElements):
1912 jfenwick 2635 if (functionSpaceType_target==FaceElements) {
1913     return true;
1914     } else if (functionSpaceType_target==ReducedFaceElements) {
1915     return true;
1916     } else {
1917     return false;
1918     }
1919 phornby 1628 case(ReducedFaceElements):
1920 jfenwick 2635 if (functionSpaceType_target==ReducedFaceElements) {
1921     return true;
1922     } else {
1923     return false;
1924     }
1925 phornby 1628 case(Points):
1926 jfenwick 2635 if (functionSpaceType_target==Points) {
1927     return true;
1928     } else {
1929     return false;
1930     }
1931 phornby 1628 case(ContactElementsZero):
1932     case(ContactElementsOne):
1933 jfenwick 2635 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1934     return true;
1935     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1936     return true;
1937     } else {
1938     return false;
1939     }
1940 phornby 1628 case(ReducedContactElementsZero):
1941     case(ReducedContactElementsOne):
1942 jfenwick 2635 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1943     return true;
1944     } else {
1945     return false;
1946     }
1947 phornby 1628 case(DegreesOfFreedom):
1948 jfenwick 2635 switch(functionSpaceType_target) {
1949     case(ReducedDegreesOfFreedom):
1950     case(DegreesOfFreedom):
1951     case(Nodes):
1952     case(ReducedNodes):
1953     case(Elements):
1954     case(ReducedElements):
1955     case(Points):
1956     case(FaceElements):
1957     case(ReducedFaceElements):
1958     case(ContactElementsZero):
1959     case(ReducedContactElementsZero):
1960     case(ContactElementsOne):
1961     case(ReducedContactElementsOne):
1962     return true;
1963     default:
1964     stringstream temp;
1965     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1966     throw FinleyAdapterException(temp.str());
1967     }
1968     break;
1969 phornby 1628 case(ReducedDegreesOfFreedom):
1970     switch(functionSpaceType_target) {
1971 jfenwick 2635 case(ReducedDegreesOfFreedom):
1972     case(ReducedNodes):
1973     case(Elements):
1974     case(ReducedElements):
1975     case(FaceElements):
1976     case(ReducedFaceElements):
1977     case(Points):
1978     case(ContactElementsZero):
1979     case(ReducedContactElementsZero):
1980     case(ContactElementsOne):
1981     case(ReducedContactElementsOne):
1982     return true;
1983     case(Nodes):
1984     case(DegreesOfFreedom):
1985     return false;
1986     default:
1987     stringstream temp;
1988     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1989     throw FinleyAdapterException(temp.str());
1990     }
1991     break;
1992 phornby 1628 default:
1993     stringstream temp;
1994     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1995     throw FinleyAdapterException(temp.str());
1996     break;
1997     }
1998     checkFinleyError();
1999     return false;
2000 jgs 82 }
2001 jgs 149
2002 jgs 82 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
2003     {
2004 phornby 1628 return false;
2005 jgs 82 }
2006 jgs 104
2007 jgs 121 bool MeshAdapter::operator==(const AbstractDomain& other) const
2008 jgs 82 {
2009 phornby 1628 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2010     if (temp!=0) {
2011     return (m_finleyMesh==temp->m_finleyMesh);
2012     } else {
2013     return false;
2014     }
2015 jgs 82 }
2016    
2017 jgs 121 bool MeshAdapter::operator!=(const AbstractDomain& other) const
2018 jgs 104 {
2019 phornby 1628 return !(operator==(other));
2020 jgs 104 }
2021    
2022 gross 1859 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2023 jgs 102 {
2024 gross 2315 Finley_Mesh* mesh=m_finleyMesh.get();
2025     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2026 jgs 150 checkPasoError();
2027 jgs 102 return out;
2028     }
2029 gross 1859 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2030     {
2031 gross 2315 Finley_Mesh* mesh=m_finleyMesh.get();
2032     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2033 gross 1859 checkPasoError();
2034     return out;
2035     }
2036 jgs 149
2037 jgs 480 escript::Data MeshAdapter::getX() const
2038 jgs 102 {
2039 phornby 1628 return continuousFunction(asAbstractContinuousDomain()).getX();
2040 jgs 102 }
2041 jgs 149
2042 jgs 480 escript::Data MeshAdapter::getNormal() const
2043 jgs 102 {
2044 phornby 1628 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
2045 jgs 102 }
2046 jgs 149
2047 jgs 480 escript::Data MeshAdapter::getSize() const
2048 jgs 102 {
2049 jfenwick 1872 return escript::function(asAbstractContinuousDomain()).getSize();
2050 jgs 102 }
2051    
2052 jfenwick 2487 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2053 gross 1062 {
2054 phornby 1628 int *out = NULL;
2055     Finley_Mesh* mesh=m_finleyMesh.get();
2056     switch (functionSpaceType) {
2057     case(Nodes):
2058     out=mesh->Nodes->Id;
2059     break;
2060     case(ReducedNodes):
2061     out=mesh->Nodes->reducedNodesId;
2062     break;
2063     case(Elements):
2064     out=mesh->Elements->Id;
2065     break;
2066     case(ReducedElements):
2067     out=mesh->Elements->Id;
2068     break;
2069     case(FaceElements):
2070     out=mesh->FaceElements->Id;
2071     break;
2072     case(ReducedFaceElements):
2073     out=mesh->FaceElements->Id;
2074     break;
2075     case(Points):
2076     out=mesh->Points->Id;
2077     break;
2078     case(ContactElementsZero):
2079     out=mesh->ContactElements->Id;
2080     break;
2081     case(ReducedContactElementsZero):
2082     out=mesh->ContactElements->Id;
2083     break;
2084     case(ContactElementsOne):
2085     out=mesh->ContactElements->Id;
2086     break;
2087     case(ReducedContactElementsOne):
2088     out=mesh->ContactElements->Id;
2089     break;
2090     case(DegreesOfFreedom):
2091     out=mesh->Nodes->degreesOfFreedomId;
2092     break;
2093     case(ReducedDegreesOfFreedom):
2094     out=mesh->Nodes->reducedDegreesOfFreedomId;
2095     break;
2096     default:
2097 gross 1062 stringstream temp;
2098     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2099     throw FinleyAdapterException(temp.str());
2100     break;
2101 phornby 1628 }
2102     return out;
2103 gross 1062 }
2104 jgs 110 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2105     {
2106 phornby 1628 int out=0;
2107     Finley_Mesh* mesh=m_finleyMesh.get();
2108     switch (functionSpaceType) {
2109     case(Nodes):
2110     out=mesh->Nodes->Tag[sampleNo];
2111     break;
2112     case(ReducedNodes):
2113     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2114     break;
2115     case(Elements):
2116     out=mesh->Elements->Tag[sampleNo];
2117     break;
2118     case(ReducedElements):
2119     out=mesh->Elements->Tag[sampleNo];
2120     break;
2121     case(FaceElements):
2122     out=mesh->FaceElements->Tag[sampleNo];
2123     break;
2124     case(ReducedFaceElements):
2125     out=mesh->FaceElements->Tag[sampleNo];
2126     break;
2127     case(Points):
2128     out=mesh->Points->Tag[sampleNo];
2129     break;
2130     case(ContactElementsZero):
2131     out=mesh->ContactElements->Tag[sampleNo];
2132     break;
2133     case(ReducedContactElementsZero):
2134     out=mesh->ContactElements->Tag[sampleNo];
2135     break;
2136     case(ContactElementsOne):
2137     out=mesh->ContactElements->Tag[sampleNo];
2138     break;
2139     case(ReducedContactElementsOne):
2140     out=mesh->ContactElements->Tag[sampleNo];
2141     break;
2142     case(DegreesOfFreedom):
2143     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2144     break;
2145     case(ReducedDegreesOfFreedom):
2146     throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2147     break;
2148     default:
2149     stringstream temp;
2150     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2151     throw FinleyAdapterException(temp.str());
2152     break;
2153     }
2154     return out;
2155 jgs 110 }
2156    
2157 gross 1062
2158 gross 767