/[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 3344 - (hide annotations)
Thu Nov 11 23:26:52 2010 UTC (8 years, 6 months ago) by caltinay
File size: 85874 byte(s)
Phew!
-escript, finley, and dudley now uses weipa's saveVTK implementation
-moved tests from finley to weipa accordingly; dudley still to do
-rebaselined all test files
-fixed a few issues in weipa.saveVTK, e.g. saving metadata without schema
-added a deprecation warning to esys.escript.util.saveVTK
-todo: change doco, tests and other places to use weipa.saveVTK


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