/[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 5122 - (hide annotations)
Thu Aug 21 04:32:39 2014 UTC (4 years, 7 months ago) by caltinay
File size: 91592 byte(s)
Fast-forward to latest trunk to be able to use both Paso or cusp.

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