/[escript]/branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Annotation of /branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2635 - (hide annotations)
Thu Aug 27 04:54:41 2009 UTC (9 years, 7 months ago) by jfenwick
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 88534 byte(s)
A bunch of changes related to saveDataCSV.
[Not completed or unit tested yet]

Added saveDataCSV to util.py
AbstractDomain (and MeshAdapter) have a commonFunctionSpace method to 
take a group of FunctionSpaces and return something they can all be interpolated to.

Added pointToStream() in DataTypes to help print points.

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