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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3914 - (hide annotations)
Wed Jun 20 04:51:47 2012 UTC (6 years, 9 months ago) by jfenwick
File size: 90718 byte(s)
Only allow setX from ContinuousFunctions.
This addresses mantis 644.
1 jgs 472
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 3911 * Copyright (c) 2003-2012 by University of Queensland
5 ksteube 1811 * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 ksteube 1312
14 jfenwick 3675 #include <pasowrap/PasoException.h>
15     #include <pasowrap/TransportProblemAdapter.h>
16 jgs 203 #include "MeshAdapter.h"
17 robwdcock 682 #include "escript/Data.h"
18     #include "escript/DataFactory.h"
19 ksteube 1343 #ifdef USE_NETCDF
20     #include <netcdfcpp.h>
21     #endif
22 ksteube 1312 extern "C" {
23 phornby 2078 #include "esysUtils/blocktimer.h"
24 ksteube 1312 }
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 caltinay 2151 //
1524     // Helper for the save* methods. Extracts optional data variable names and the
1525     // corresponding pointers from python dictionary. Caller must free arrays.
1526     //
1527     void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1528 jgs 153 {
1529 caltinay 2151 numData = boost::python::extract<int>(arg.attr("__len__")());
1530 phornby 1628 /* win32 refactor */
1531 caltinay 2151 names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1532     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1533     dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1534 jgs 153
1535 phornby 1628 boost::python::list keys=arg.keys();
1536 caltinay 2151 for (int i=0; i<numData; ++i) {
1537 caltinay 2778 string n=boost::python::extract<string>(keys[i]);
1538 phornby 1628 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1539 jfenwick 1872 if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1540 caltinay 2151 throw FinleyAdapterException("Error: Data must be defined on same Domain");
1541     data[i] = d.getDataC();
1542     dataPtr[i] = &(data[i]);
1543 caltinay 2150 names[i] = TMPMEMALLOC(n.length()+1, char);
1544     strcpy(names[i], n.c_str());
1545 phornby 1628 }
1546 caltinay 2151 }
1547    
1548     //
1549     // saves mesh and optionally data arrays in openDX format
1550     //
1551 caltinay 2778 void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1552 caltinay 2151 {
1553     int num_data;
1554     char **names;
1555     escriptDataC *data;
1556     escriptDataC **ptr_data;
1557    
1558     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1559     Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1560 phornby 1628 checkFinleyError();
1561    
1562     /* win32 refactor */
1563     TMPMEMFREE(data);
1564     TMPMEMFREE(ptr_data);
1565 caltinay 2150 for(int i=0; i<num_data; i++)
1566 phornby 1628 {
1567     TMPMEMFREE(names[i]);
1568     }
1569     TMPMEMFREE(names);
1570 woo409 757
1571 phornby 1628 return;
1572 jgs 82 }
1573 jgs 149
1574 caltinay 2151 //
1575     // saves mesh and optionally data arrays in VTK format
1576     //
1577 caltinay 2778 void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg, const string& metadata, const string& metadata_schema) const
1578 jgs 153 {
1579 caltinay 3355 boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1580     pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1581     metadata, metadata_schema, arg);
1582 caltinay 2151 }
1583 woo409 757
1584 jfenwick 2642 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1585     {
1586 jfenwick 3259 #ifdef ESYS_MPI
1587 jfenwick 2642 index_t myFirstNode=0, myLastNode=0, k=0;
1588     index_t* globalNodeIndex=0;
1589     Finley_Mesh* mesh_p=m_finleyMesh.get();
1590     if (fs_code == FINLEY_REDUCED_NODES)
1591     {
1592     myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1593     myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1594     globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1595     }
1596     else
1597     {
1598     myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1599     myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1600     globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1601     }
1602     k=globalNodeIndex[id];
1603     return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1604 jfenwick 2644 #endif
1605     return true;
1606 jfenwick 2642 }
1607    
1608    
1609    
1610 caltinay 2151 //
1611     // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1612     //
1613 jfenwick 3259 ASM_ptr MeshAdapter::newSystemMatrix(
1614 phornby 1628 const int row_blocksize,
1615     const escript::FunctionSpace& row_functionspace,
1616     const int column_blocksize,
1617     const escript::FunctionSpace& column_functionspace,
1618     const int type) const
1619 jgs 82 {
1620 phornby 1628 int reduceRowOrder=0;
1621     int reduceColOrder=0;
1622     // is the domain right?
1623 jfenwick 1872 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1624 phornby 1628 if (row_domain!=*this)
1625     throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1626 jfenwick 1872 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1627 phornby 1628 if (col_domain!=*this)
1628     throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1629     // is the function space type right
1630     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1631     reduceRowOrder=0;
1632     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1633     reduceRowOrder=1;
1634     } else {
1635     throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1636     }
1637     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1638     reduceColOrder=0;
1639     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1640     reduceColOrder=1;
1641     } else {
1642     throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1643     }
1644     // generate matrix:
1645    
1646     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1647     checkFinleyError();
1648     Paso_SystemMatrix* fsystemMatrix;
1649     int trilinos = 0;
1650     if (trilinos) {
1651 ksteube 1312 #ifdef TRILINOS
1652     /* Allocation Epetra_VrbMatrix here */
1653     #endif
1654 phornby 1628 }
1655     else {
1656 gross 2551 fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1657 phornby 1628 }
1658     checkPasoError();
1659     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1660 jfenwick 3259 SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1661     return ASM_ptr(sma);
1662     // return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1663 jgs 82 }
1664 caltinay 2151
1665     //
1666 gross 1366 // creates a TransportProblemAdapter
1667 caltinay 2151 //
1668 jfenwick 3259 ATP_ptr MeshAdapter::newTransportProblem(
1669 phornby 1628 const int blocksize,
1670     const escript::FunctionSpace& functionspace,
1671     const int type) const
1672 gross 1366 {
1673 phornby 1628 int reduceOrder=0;
1674     // is the domain right?
1675 jfenwick 1872 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1676 phornby 1628 if (domain!=*this)
1677     throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1678     // is the function space type right
1679     if (functionspace.getTypeCode()==DegreesOfFreedom) {
1680     reduceOrder=0;
1681     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1682     reduceOrder=1;
1683     } else {
1684     throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1685     }
1686     // generate matrix:
1687    
1688     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1689     checkFinleyError();
1690 gross 2987 Paso_TransportProblem* transportProblem;
1691 gross 3793 transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1692 phornby 1628 checkPasoError();
1693     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1694 gross 3793 TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1695 jfenwick 3259 return ATP_ptr(tpa);
1696 gross 3793 // return TransportProblemAdapter(transportProblem, blocksize,functionspace);
1697 gross 1366 }
1698 jgs 149
1699 jgs 82 //
1700     // vtkObject MeshAdapter::createVtkObject() const
1701     // TODO:
1702     //
1703 jgs 149
1704 jgs 82 //
1705     // returns true if data at the atom_type is considered as being cell centered:
1706     bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1707     {
1708 phornby 1628 switch(functionSpaceCode) {
1709     case(Nodes):
1710     case(DegreesOfFreedom):
1711     case(ReducedDegreesOfFreedom):
1712     return false;
1713     break;
1714     case(Elements):
1715     case(FaceElements):
1716     case(Points):
1717     case(ContactElementsZero):
1718     case(ContactElementsOne):
1719     case(ReducedElements):
1720     case(ReducedFaceElements):
1721     case(ReducedContactElementsZero):
1722     case(ReducedContactElementsOne):
1723     return true;
1724     break;
1725     default:
1726     stringstream temp;
1727     temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1728     throw FinleyAdapterException(temp.str());
1729     break;
1730     }
1731     checkFinleyError();
1732     return false;
1733 jgs 82 }
1734 jgs 149
1735 jfenwick 2635 bool
1736 caltinay 2778 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1737 jfenwick 2635 {
1738     /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1739     class 1: DOF <-> Nodes
1740     class 2: ReducedDOF <-> ReducedNodes
1741     class 3: Points
1742     class 4: Elements
1743     class 5: ReducedElements
1744     class 6: FaceElements
1745     class 7: ReducedFaceElements
1746     class 8: ContactElementZero <-> ContactElementOne
1747     class 9: ReducedContactElementZero <-> ReducedContactElementOne
1748    
1749     There is also a set of lines. Interpolation is possible down a line but not between lines.
1750     class 1 and 2 belong to all lines so aren't considered.
1751     line 0: class 3
1752     line 1: class 4,5
1753     line 2: class 6,7
1754     line 3: class 8,9
1755    
1756     For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1757     eg hasnodes is true if we have at least one instance of Nodes.
1758     */
1759 caltinay 2940 if (fs.empty())
1760 jfenwick 2635 {
1761 caltinay 2778 return false;
1762 jfenwick 2635 }
1763 caltinay 2778 vector<int> hasclass(10);
1764     vector<int> hasline(4);
1765 jfenwick 2635 bool hasnodes=false;
1766     bool hasrednodes=false;
1767     bool hascez=false;
1768     bool hasrcez=false;
1769     for (int i=0;i<fs.size();++i)
1770     {
1771     switch(fs[i])
1772     {
1773     case(Nodes): hasnodes=true; // no break is deliberate
1774     case(DegreesOfFreedom):
1775     hasclass[1]=1;
1776     break;
1777     case(ReducedNodes): hasrednodes=true; // no break is deliberate
1778     case(ReducedDegreesOfFreedom):
1779     hasclass[2]=1;
1780     break;
1781     case(Points):
1782     hasline[0]=1;
1783     hasclass[3]=1;
1784     break;
1785     case(Elements):
1786     hasclass[4]=1;
1787     hasline[1]=1;
1788     break;
1789     case(ReducedElements):
1790     hasclass[5]=1;
1791     hasline[1]=1;
1792     break;
1793     case(FaceElements):
1794     hasclass[6]=1;
1795     hasline[2]=1;
1796     break;
1797     case(ReducedFaceElements):
1798     hasclass[7]=1;
1799     hasline[2]=1;
1800     break;
1801     case(ContactElementsZero): hascez=true; // no break is deliberate
1802     case(ContactElementsOne):
1803     hasclass[8]=1;
1804     hasline[3]=1;
1805     break;
1806     case(ReducedContactElementsZero): hasrcez=true; // no break is deliberate
1807     case(ReducedContactElementsOne):
1808     hasclass[9]=1;
1809     hasline[3]=1;
1810     break;
1811     default:
1812     return false;
1813     }
1814     }
1815     int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1816     // fail if we have more than one leaf group
1817    
1818     if (totlines>1)
1819     {
1820     return false; // there are at least two branches we can't interpolate between
1821     }
1822     else if (totlines==1)
1823     {
1824     if (hasline[0]==1) // we have points
1825     {
1826     resultcode=Points;
1827     }
1828     else if (hasline[1]==1)
1829     {
1830     if (hasclass[5]==1)
1831     {
1832     resultcode=ReducedElements;
1833     }
1834     else
1835     {
1836     resultcode=Elements;
1837     }
1838     }
1839     else if (hasline[2]==1)
1840     {
1841 jfenwick 2644 if (hasclass[7]==1)
1842 jfenwick 2635 {
1843     resultcode=ReducedFaceElements;
1844     }
1845     else
1846     {
1847     resultcode=FaceElements;
1848     }
1849     }
1850     else // so we must be in line3
1851     {
1852     if (hasclass[9]==1)
1853     {
1854     // need something from class 9
1855     resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1856     }
1857     else
1858     {
1859     // something from class 8
1860     resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1861     }
1862     }
1863     }
1864     else // totlines==0
1865     {
1866     if (hasclass[2]==1)
1867     {
1868     // something from class 2
1869     resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1870     }
1871     else
1872     { // something from class 1
1873     resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1874     }
1875     }
1876     return true;
1877     }
1878    
1879 jgs 82 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1880     {
1881 phornby 1628 switch(functionSpaceType_source) {
1882     case(Nodes):
1883 jfenwick 2635 switch(functionSpaceType_target) {
1884     case(Nodes):
1885     case(ReducedNodes):
1886     case(ReducedDegreesOfFreedom):
1887     case(DegreesOfFreedom):
1888     case(Elements):
1889     case(ReducedElements):
1890     case(FaceElements):
1891     case(ReducedFaceElements):
1892     case(Points):
1893     case(ContactElementsZero):
1894     case(ReducedContactElementsZero):
1895     case(ContactElementsOne):
1896     case(ReducedContactElementsOne):
1897     return true;
1898     default:
1899     stringstream temp;
1900     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1901     throw FinleyAdapterException(temp.str());
1902 phornby 1628 }
1903     break;
1904     case(ReducedNodes):
1905 jfenwick 2635 switch(functionSpaceType_target) {
1906     case(ReducedNodes):
1907     case(ReducedDegreesOfFreedom):
1908     case(Elements):
1909     case(ReducedElements):
1910     case(FaceElements):
1911     case(ReducedFaceElements):
1912     case(Points):
1913     case(ContactElementsZero):
1914     case(ReducedContactElementsZero):
1915     case(ContactElementsOne):
1916     case(ReducedContactElementsOne):
1917     return true;
1918     case(Nodes):
1919     case(DegreesOfFreedom):
1920     return false;
1921     default:
1922     stringstream temp;
1923     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1924     throw FinleyAdapterException(temp.str());
1925 phornby 1628 }
1926     break;
1927     case(Elements):
1928 jfenwick 2635 if (functionSpaceType_target==Elements) {
1929     return true;
1930     } else if (functionSpaceType_target==ReducedElements) {
1931     return true;
1932     } else {
1933     return false;
1934     }
1935 phornby 1628 case(ReducedElements):
1936 jfenwick 2635 if (functionSpaceType_target==ReducedElements) {
1937     return true;
1938     } else {
1939     return false;
1940     }
1941 phornby 1628 case(FaceElements):
1942 jfenwick 2635 if (functionSpaceType_target==FaceElements) {
1943     return true;
1944     } else if (functionSpaceType_target==ReducedFaceElements) {
1945     return true;
1946     } else {
1947     return false;
1948     }
1949 phornby 1628 case(ReducedFaceElements):
1950 jfenwick 2635 if (functionSpaceType_target==ReducedFaceElements) {
1951     return true;
1952     } else {
1953     return false;
1954     }
1955 phornby 1628 case(Points):
1956 jfenwick 2635 if (functionSpaceType_target==Points) {
1957     return true;
1958     } else {
1959     return false;
1960     }
1961 phornby 1628 case(ContactElementsZero):
1962     case(ContactElementsOne):
1963 jfenwick 2635 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1964     return true;
1965     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1966     return true;
1967     } else {
1968     return false;
1969     }
1970 phornby 1628 case(ReducedContactElementsZero):
1971     case(ReducedContactElementsOne):
1972 jfenwick 2635 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1973     return true;
1974     } else {
1975     return false;
1976     }
1977 phornby 1628 case(DegreesOfFreedom):
1978 jfenwick 2635 switch(functionSpaceType_target) {
1979     case(ReducedDegreesOfFreedom):
1980     case(DegreesOfFreedom):
1981     case(Nodes):
1982     case(ReducedNodes):
1983     case(Elements):
1984     case(ReducedElements):
1985     case(Points):
1986     case(FaceElements):
1987     case(ReducedFaceElements):
1988     case(ContactElementsZero):
1989     case(ReducedContactElementsZero):
1990     case(ContactElementsOne):
1991     case(ReducedContactElementsOne):
1992     return true;
1993     default:
1994     stringstream temp;
1995     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1996     throw FinleyAdapterException(temp.str());
1997     }
1998     break;
1999 phornby 1628 case(ReducedDegreesOfFreedom):
2000     switch(functionSpaceType_target) {
2001 jfenwick 2635 case(ReducedDegreesOfFreedom):
2002     case(ReducedNodes):
2003     case(Elements):
2004     case(ReducedElements):
2005     case(FaceElements):
2006     case(ReducedFaceElements):
2007     case(Points):
2008     case(ContactElementsZero):
2009     case(ReducedContactElementsZero):
2010     case(ContactElementsOne):
2011     case(ReducedContactElementsOne):
2012     return true;
2013     case(Nodes):
2014     case(DegreesOfFreedom):
2015     return false;
2016     default:
2017     stringstream temp;
2018     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
2019     throw FinleyAdapterException(temp.str());
2020     }
2021     break;
2022 phornby 1628 default:
2023     stringstream temp;
2024     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
2025     throw FinleyAdapterException(temp.str());
2026     break;
2027     }
2028     checkFinleyError();
2029     return false;
2030 jgs 82 }
2031 jgs 149
2032 jgs 82 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
2033     {
2034 phornby 1628 return false;
2035 jgs 82 }
2036 jgs 104
2037 jgs 121 bool MeshAdapter::operator==(const AbstractDomain& other) const
2038 jgs 82 {
2039 phornby 1628 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2040     if (temp!=0) {
2041     return (m_finleyMesh==temp->m_finleyMesh);
2042     } else {
2043     return false;
2044     }
2045 jgs 82 }
2046    
2047 jgs 121 bool MeshAdapter::operator!=(const AbstractDomain& other) const
2048 jgs 104 {
2049 phornby 1628 return !(operator==(other));
2050 jgs 104 }
2051    
2052 gross 1859 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2053 jgs 102 {
2054 gross 2315 Finley_Mesh* mesh=m_finleyMesh.get();
2055 caltinay 3681 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2056     package, symmetry, mesh->MPIInfo);
2057 jgs 102 }
2058 caltinay 3681
2059 gross 1859 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2060     {
2061 gross 2315 Finley_Mesh* mesh=m_finleyMesh.get();
2062 caltinay 3681 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2063     package, symmetry, mesh->MPIInfo);
2064 gross 1859 }
2065 jgs 149
2066 jgs 480 escript::Data MeshAdapter::getX() const
2067 jgs 102 {
2068 jfenwick 3259 return continuousFunction(*this).getX();
2069 jgs 102 }
2070 jgs 149
2071 jgs 480 escript::Data MeshAdapter::getNormal() const
2072 jgs 102 {
2073 jfenwick 3259 return functionOnBoundary(*this).getNormal();
2074 jgs 102 }
2075 jgs 149
2076 jgs 480 escript::Data MeshAdapter::getSize() const
2077 jgs 102 {
2078 jfenwick 3259 return escript::function(*this).getSize();
2079 jgs 102 }
2080    
2081 jfenwick 2487 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2082 gross 1062 {
2083 phornby 1628 int *out = NULL;
2084     Finley_Mesh* mesh=m_finleyMesh.get();
2085     switch (functionSpaceType) {
2086     case(Nodes):
2087     out=mesh->Nodes->Id;
2088     break;
2089     case(ReducedNodes):
2090     out=mesh->Nodes->reducedNodesId;
2091     break;
2092     case(Elements):
2093     out=mesh->Elements->Id;
2094     break;
2095     case(ReducedElements):
2096     out=mesh->Elements->Id;
2097     break;
2098     case(FaceElements):
2099     out=mesh->FaceElements->Id;
2100     break;
2101     case(ReducedFaceElements):
2102     out=mesh->FaceElements->Id;
2103     break;
2104     case(Points):
2105     out=mesh->Points->Id;
2106     break;
2107     case(ContactElementsZero):
2108     out=mesh->ContactElements->Id;
2109     break;
2110     case(ReducedContactElementsZero):
2111     out=mesh->ContactElements->Id;
2112     break;
2113     case(ContactElementsOne):
2114     out=mesh->ContactElements->Id;
2115     break;
2116     case(ReducedContactElementsOne):
2117     out=mesh->ContactElements->Id;
2118     break;
2119     case(DegreesOfFreedom):
2120     out=mesh->Nodes->degreesOfFreedomId;
2121     break;
2122     case(ReducedDegreesOfFreedom):
2123     out=mesh->Nodes->reducedDegreesOfFreedomId;
2124     break;
2125     default:
2126 gross 1062 stringstream temp;
2127     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2128     throw FinleyAdapterException(temp.str());
2129     break;
2130 phornby 1628 }
2131     return out;
2132 gross 1062 }
2133 jgs 110 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2134     {
2135 phornby 1628 int out=0;
2136     Finley_Mesh* mesh=m_finleyMesh.get();
2137     switch (functionSpaceType) {
2138     case(Nodes):
2139     out=mesh->Nodes->Tag[sampleNo];
2140     break;
2141     case(ReducedNodes):
2142     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2143     break;
2144     case(Elements):
2145     out=mesh->Elements->Tag[sampleNo];
2146     break;
2147     case(ReducedElements):
2148     out=mesh->Elements->Tag[sampleNo];
2149     break;
2150     case(FaceElements):
2151     out=mesh->FaceElements->Tag[sampleNo];
2152     break;
2153     case(ReducedFaceElements):
2154     out=