/[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 3993 - (hide annotations)
Wed Sep 26 07:00:55 2012 UTC (6 years, 9 months ago) by caltinay
File size: 91124 byte(s)
Addressing mantis 651. Test for finley will fail under MPI since there is
sth wrong with the way ReducedContinuousFunction samples are interrogated.

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