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