/[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 3317 - (hide annotations)
Thu Oct 28 00:50:41 2010 UTC (8 years, 10 months ago) by caltinay
File size: 86021 byte(s)
Removed bogus mpi.h includes, replaced others by our esysUtils wrapper
and rearranged so that the wrapper always comes before netcdf which fixes
linking problems when disabling mpi on systems with netCDF 4.x.

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