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

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

Parent Directory Parent Directory | Revision Log Revision Log


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