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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3269 - (hide annotations)
Wed Oct 13 03:21:50 2010 UTC (8 years, 7 months ago) by jfenwick
File size: 86093 byte(s)
Fixed some intel compiler warnings.
Put linearPDEs back the way it was and present a common interface for dudley and finley (as per Lutz)

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