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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5136 - (hide annotations)
Tue Sep 9 07:13:55 2014 UTC (4 years, 6 months ago) by caltinay
File size: 91420 byte(s)
ripley now supports paso solvers again and returns an appropriate matrix type
id. Changed the getSystemMatrixTypeId() method to take a full SolverBuddy
instance and made some other simplifications.

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