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

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

Parent Directory Parent Directory | Revision Log Revision Log


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