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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4687 - (hide annotations)
Wed Feb 19 00:03:29 2014 UTC (5 years, 4 months ago) by jfenwick
File size: 91259 byte(s)
Remove randomFill python method from ripley domains.
All random data objects (for all domain types) should be generated 
using esys.escript.RandomData()

The only filtered random we have is gaussian on ripley but
it is triggered by passing the tuple as the last arg of RandomData().

While the interface is a bit more complicated (in that you always need
 to pass in shape and functionspace) it does mean we have a 
common interface for all domains. 

Removed randomFill from DataExpanded.
The reasoning behind this is to force domains to call the util function
themselves and enforce whatever consistancy requirements they have.

Added version of blocktools to deal with 2D case in Ripley.
Use blocktools for the 2D transfers [This was cleaner than modifying the
previous implementation to deal with variable shaped points].

Note that under MPI, ripley can not generate random data (even unfiltered)
if any of its per rank dimensions is <4 elements on any side.

Unit tests for these calls are in but some extra checks still needed.



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