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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4496 - (hide annotations)
Mon Jul 15 06:53:44 2013 UTC (5 years, 8 months ago) by caltinay
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 90629 byte(s)
finley (WIP):
-moved all of finley into its namespace
-introduced some shared pointers
-Mesh is now a class
-other bits and pieces...

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     Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with netCDF. Please contact your installation manager.");
504     #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 caltinay 4441 Assemble_PDE(mesh->Nodes, mesh->FaceElements, S, rhs,
758     escript::Data(), escript::Data(), escript::Data(), d,
759     escript::Data(), y);
760 caltinay 4408 checkFinleyError();
761 gross 3522
762 caltinay 4441 Assemble_PDE(mesh->Nodes, mesh->ContactElements, S, rhs,
763     escript::Data(), escript::Data(), escript::Data(), d_contact,
764     escript::Data(), y_contact);
765 caltinay 4408 checkFinleyError();
766    
767 caltinay 4441 Assemble_PDE(mesh->Nodes, mesh->Points, S, rhs, escript::Data(),
768     escript::Data(), escript::Data(), d_dirac, escript::Data(), y_dirac);
769 caltinay 4408 checkFinleyError();
770 jgs 82 }
771 jgs 149
772 caltinay 4408 void MeshAdapter::addPDEToLumpedSystem(escript::Data& mat,
773     const escript::Data& D, const escript::Data& d,
774     const escript::Data& d_dirac, const bool useHRZ) const
775 gross 1204 {
776 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
777 caltinay 4441 Assemble_LumpedSystem(mesh->Nodes, mesh->Elements, mat, D, useHRZ);
778 caltinay 4408 checkFinleyError();
779 gross 1204
780 caltinay 4441 Assemble_LumpedSystem(mesh->Nodes, mesh->FaceElements, mat, d, useHRZ);
781 caltinay 4408 checkFinleyError();
782 gross 3522
783 caltinay 4441 Assemble_LumpedSystem(mesh->Nodes, mesh->Points, mat, d_dirac, useHRZ);
784 caltinay 4408 checkFinleyError();
785 gross 1204 }
786    
787 jgs 82 //
788 jgs 102 // adds linear PDE of second order into the right hand side only
789     //
790 caltinay 4408 void MeshAdapter::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
791     const escript::Data& Y, const escript::Data& y,
792     const escript::Data& y_contact, const escript::Data& y_dirac) const
793 jgs 102 {
794 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
795 caltinay 4441 Assemble_PDE(mesh->Nodes, mesh->Elements, 0, rhs,
796     escript::Data(), escript::Data(), escript::Data(), escript::Data(),
797     X, Y);
798 caltinay 4408 checkFinleyError();
799 jgs 147
800 caltinay 4441 Assemble_PDE(mesh->Nodes, mesh->FaceElements, 0, rhs,
801     escript::Data(), escript::Data(), escript::Data(), escript::Data(),
802     escript::Data(), y);
803 caltinay 4408 checkFinleyError();
804 jgs 147
805 caltinay 4441 Assemble_PDE(mesh->Nodes, mesh->ContactElements, 0, rhs,
806     escript::Data(), escript::Data(), escript::Data(),
807     escript::Data(), escript::Data(), y_contact);
808 caltinay 4408 checkFinleyError();
809 gross 3522
810 caltinay 4441 Assemble_PDE(mesh->Nodes, mesh->Points, 0, rhs,
811     escript::Data(), escript::Data(), escript::Data(), escript::Data(),
812     escript::Data(), y_dirac);
813 caltinay 4408 checkFinleyError();
814     }
815 gross 3522
816 gross 1367 //
817     // adds PDE of second order into a transport problem
818     //
819     void MeshAdapter::addPDEToTransportProblem(
820 caltinay 4408 escript::AbstractTransportProblem& tp, escript::Data& source,
821     const escript::Data& M, const escript::Data& A, const escript::Data& B,
822     const escript::Data& C, const escript::Data& D, const escript::Data& X,
823     const escript::Data& Y, const escript::Data& d, const escript::Data& y,
824     const escript::Data& d_contact, const escript::Data& y_contact,
825     const escript::Data& d_dirac, const escript::Data& y_dirac) const
826 gross 1367 {
827 caltinay 4408 TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
828     if (!tpa)
829     throw FinleyAdapterException("finley only supports Paso transport problems.");
830 jfenwick 3269
831 caltinay 4408 source.expand();
832 jfenwick 3269
833 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
834 caltinay 4408 Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
835 jgs 149
836 caltinay 4441 Assemble_PDE(mesh->Nodes, mesh->Elements, _tp->mass_matrix, source,
837     escript::Data(), escript::Data(), escript::Data(),
838     M, escript::Data(), escript::Data());
839 caltinay 4408 checkFinleyError();
840 gross 3522
841 caltinay 4441 Assemble_PDE(mesh->Nodes, mesh->Elements, _tp->transport_matrix,
842     source, A, B, C, D, X, Y);
843 caltinay 4408 checkFinleyError();
844 phornby 1628
845 caltinay 4441 Assemble_PDE(mesh->Nodes, mesh->FaceElements, _tp->transport_matrix,
846     source, escript::Data(), escript::Data(),
847     escript::Data(), d, escript::Data(), y);
848 caltinay 4408 checkFinleyError();
849 gross 1367
850 caltinay 4441 Assemble_PDE(mesh->Nodes, mesh->ContactElements,
851     _tp->transport_matrix, source, escript::Data(),
852     escript::Data(), escript::Data(), d_contact,
853     escript::Data(), y_contact);
854 caltinay 4408 checkFinleyError();
855 gross 1367
856 caltinay 4441 Assemble_PDE(mesh->Nodes, mesh->Points, _tp->transport_matrix,
857     source, escript::Data(), escript::Data(),
858     escript::Data(), d_dirac, escript::Data(), y_dirac);
859 caltinay 4408 checkFinleyError();
860 gross 1367 }
861    
862 jgs 102 //
863 caltinay 4408 // interpolates data between different function spaces
864 jgs 82 //
865 caltinay 4408 void MeshAdapter::interpolateOnDomain(escript::Data& target, const escript::Data& in) const
866 jgs 82 {
867 caltinay 4408 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
868     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
869     if (inDomain!=*this)
870     throw FinleyAdapterException("Error - Illegal domain of interpolant.");
871     if (targetDomain!=*this)
872     throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
873 jgs 82
874 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
875 caltinay 4408 switch(in.getFunctionSpace().getTypeCode()) {
876     case Nodes:
877     switch(target.getFunctionSpace().getTypeCode()) {
878     case Nodes:
879     case ReducedNodes:
880     case DegreesOfFreedom:
881     case ReducedDegreesOfFreedom:
882 caltinay 4441 Assemble_CopyNodalData(mesh->Nodes, target, in);
883 caltinay 4408 break;
884     case Elements:
885     case ReducedElements:
886 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->Elements, in,target);
887 caltinay 4408 break;
888     case FaceElements:
889     case ReducedFaceElements:
890 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
891 caltinay 4408 break;
892     case Points:
893 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
894 caltinay 4408 break;
895     case ContactElementsZero:
896     case ReducedContactElementsZero:
897     case ContactElementsOne:
898     case ReducedContactElementsOne:
899 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in, target);
900 caltinay 4408 break;
901     default:
902     stringstream temp;
903     temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
904     throw FinleyAdapterException(temp.str());
905     break;
906     }
907     break;
908     case ReducedNodes:
909     switch(target.getFunctionSpace().getTypeCode()) {
910     case Nodes:
911     case ReducedNodes:
912     case DegreesOfFreedom:
913     case ReducedDegreesOfFreedom:
914 caltinay 4441 Assemble_CopyNodalData(mesh->Nodes, target, in);
915 caltinay 4408 break;
916     case Elements:
917     case ReducedElements:
918 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->Elements, in, target);
919 caltinay 4408 break;
920     case FaceElements:
921     case ReducedFaceElements:
922 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
923 caltinay 4408 break;
924     case Points:
925 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
926 caltinay 4408 break;
927     case ContactElementsZero:
928     case ReducedContactElementsZero:
929 caltinay 4410 case ContactElementsOne:
930     case ReducedContactElementsOne:
931 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in, target);
932 caltinay 4408 break;
933     default:
934     stringstream temp;
935     temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
936     throw FinleyAdapterException(temp.str());
937     break;
938     }
939     break;
940     case Elements:
941     if (target.getFunctionSpace().getTypeCode()==Elements) {
942 caltinay 4441 Assemble_CopyElementData(mesh->Elements, target, in);
943 caltinay 4408 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
944 caltinay 4441 Assemble_AverageElementData(mesh->Elements, target, in);
945 caltinay 4408 } else {
946     throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
947     }
948     break;
949     case ReducedElements:
950     if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
951 caltinay 4441 Assemble_CopyElementData(mesh->Elements, target, in);
952 caltinay 4408 } else {
953     throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
954     }
955     break;
956     case FaceElements:
957     if (target.getFunctionSpace().getTypeCode()==FaceElements) {
958 caltinay 4441 Assemble_CopyElementData(mesh->FaceElements, target, in);
959 caltinay 4408 } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
960 caltinay 4441 Assemble_AverageElementData(mesh->FaceElements, target, in);
961 caltinay 4408 } else {
962     throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
963     }
964     break;
965     case ReducedFaceElements:
966     if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
967 caltinay 4441 Assemble_CopyElementData(mesh->FaceElements, target, in);
968 caltinay 4408 } else {
969     throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
970     }
971     break;
972     case Points:
973     if (target.getFunctionSpace().getTypeCode()==Points) {
974 caltinay 4441 Assemble_CopyElementData(mesh->Points, target, in);
975 caltinay 4408 } else {
976     throw FinleyAdapterException("Error - No interpolation with data on points possible.");
977     }
978     break;
979     case ContactElementsZero:
980     case ContactElementsOne:
981     if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
982 caltinay 4441 Assemble_CopyElementData(mesh->ContactElements, target, in);
983 caltinay 4408 } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
984 caltinay 4441 Assemble_AverageElementData(mesh->ContactElements, target, in);
985 caltinay 4408 } else {
986     throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
987     }
988     break;
989     case ReducedContactElementsZero:
990     case ReducedContactElementsOne:
991     if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
992 caltinay 4441 Assemble_CopyElementData(mesh->ContactElements, target, in);
993 caltinay 4408 } else {
994     throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
995     }
996     break;
997     case DegreesOfFreedom:
998     switch(target.getFunctionSpace().getTypeCode()) {
999     case ReducedDegreesOfFreedom:
1000     case DegreesOfFreedom:
1001 caltinay 4441 Assemble_CopyNodalData(mesh->Nodes, target, in);
1002 caltinay 4408 break;
1003    
1004     case Nodes:
1005     case ReducedNodes:
1006     if (getMPISize()>1) {
1007 caltinay 4441 escript::Data in2=escript::Data(in);
1008     in2.expand();
1009     Assemble_CopyNodalData(mesh->Nodes, target, in2);
1010 caltinay 4408 } else {
1011 caltinay 4441 Assemble_CopyNodalData(mesh->Nodes, target, in);
1012 caltinay 4408 }
1013     break;
1014     case Elements:
1015     case ReducedElements:
1016     if (getMPISize()>1) {
1017 caltinay 4441 escript::Data in2=escript::Data(in, continuousFunction(*this));
1018     Assemble_interpolate(mesh->Nodes, mesh->Elements, in2, target);
1019 caltinay 4408 } else {
1020 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->Elements, in, target);
1021 caltinay 4408 }
1022     break;
1023     case FaceElements:
1024     case ReducedFaceElements:
1025     if (getMPISize()>1) {
1026 caltinay 4441 escript::Data in2=escript::Data(in, continuousFunction(*this));
1027     Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in2, target);
1028 caltinay 4408 } else {
1029 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
1030 caltinay 4408 }
1031     break;
1032     case Points:
1033     if (getMPISize()>1) {
1034 caltinay 4441 //escript::Data in2=escript::Data(in, continuousFunction(*this) );
1035 caltinay 4408 } else {
1036 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
1037 caltinay 4408 }
1038     break;
1039     case ContactElementsZero:
1040     case ContactElementsOne:
1041     case ReducedContactElementsZero:
1042     case ReducedContactElementsOne:
1043     if (getMPISize()>1) {
1044 caltinay 4441 escript::Data in2=escript::Data(in, continuousFunction(*this));
1045     Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in2, target);
1046 caltinay 4408 } else {
1047 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in, target);
1048 caltinay 4408 }
1049     break;
1050     default:
1051     stringstream temp;
1052     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1053     throw FinleyAdapterException(temp.str());
1054     }
1055     break;
1056     case ReducedDegreesOfFreedom:
1057     switch(target.getFunctionSpace().getTypeCode()) {
1058     case Nodes:
1059     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1060     case ReducedNodes:
1061     if (getMPISize()>1) {
1062 caltinay 4441 escript::Data in2=escript::Data(in);
1063     in2.expand();
1064     Assemble_CopyNodalData(mesh->Nodes, target, in2);
1065 caltinay 4408 } else {
1066 caltinay 4441 Assemble_CopyNodalData(mesh->Nodes, target, in);
1067 caltinay 4408 }
1068     break;
1069     case DegreesOfFreedom:
1070     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1071     break;
1072     case ReducedDegreesOfFreedom:
1073 caltinay 4441 Assemble_CopyNodalData(mesh->Nodes, target, in);
1074 caltinay 4408 break;
1075     case Elements:
1076     case ReducedElements:
1077     if (getMPISize()>1) {
1078 caltinay 4441 escript::Data in2=escript::Data(in, reducedContinuousFunction(*this) );
1079     Assemble_interpolate(mesh->Nodes, mesh->Elements, in2, target);
1080 caltinay 4408 } else {
1081 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->Elements, in, target);
1082 caltinay 4408 }
1083     break;
1084     case FaceElements:
1085     case ReducedFaceElements:
1086     if (getMPISize()>1) {
1087 caltinay 4441 escript::Data in2=escript::Data(in, reducedContinuousFunction(*this) );
1088     Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in2, target);
1089 caltinay 4408 } else {
1090 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
1091 caltinay 4408 }
1092     break;
1093     case Points:
1094     if (getMPISize()>1) {
1095 caltinay 4441 escript::Data in2=escript::Data(in, reducedContinuousFunction(*this));
1096     Assemble_interpolate(mesh->Nodes, mesh->Points, in2, target);
1097 caltinay 4408 } else {
1098 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
1099 caltinay 4408 }
1100     break;
1101     case ContactElementsZero:
1102     case ContactElementsOne:
1103     case ReducedContactElementsZero:
1104     case ReducedContactElementsOne:
1105     if (getMPISize()>1) {
1106 caltinay 4441 escript::Data in2=escript::Data(in, reducedContinuousFunction(*this));
1107     Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in2, target);
1108 caltinay 4408 } else {
1109 caltinay 4441 Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in, target);
1110 caltinay 4408 }
1111     break;
1112     default:
1113     stringstream temp;
1114     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1115     throw FinleyAdapterException(temp.str());
1116     break;
1117     }
1118     break;
1119     default:
1120     stringstream temp;
1121     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1122     throw FinleyAdapterException(temp.str());
1123     break;
1124     }
1125     checkFinleyError();
1126 jgs 82 }
1127    
1128     //
1129 caltinay 4408 // copies the locations of sample points into x
1130 jgs 82 //
1131 jgs 480 void MeshAdapter::setToX(escript::Data& arg) const
1132 jgs 82 {
1133 caltinay 4408 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1134     if (argDomain!=*this)
1135     throw FinleyAdapterException("Error - Illegal domain of data point locations");
1136 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
1137 caltinay 4408 // in case of appropriate function space we can do the job directly:
1138     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1139 caltinay 4441 Assemble_NodeCoordinates(mesh->Nodes, arg);
1140 caltinay 4408 } else {
1141     escript::Data tmp_data=Vector(0., continuousFunction(*this), true);
1142 caltinay 4441 Assemble_NodeCoordinates(mesh->Nodes, tmp_data);
1143 caltinay 4408 // this is then interpolated onto arg:
1144     interpolateOnDomain(arg, tmp_data);
1145     }
1146     checkFinleyError();
1147 jgs 82 }
1148 jgs 149
1149 jgs 82 //
1150 caltinay 4408 // return the normal vectors at the location of data points as a Data object
1151 jgs 82 //
1152 jgs 480 void MeshAdapter::setToNormal(escript::Data& normal) const
1153 jgs 82 {
1154 caltinay 4408 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1155     if (normalDomain!=*this)
1156     throw FinleyAdapterException("Error - Illegal domain of normal locations");
1157 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
1158 caltinay 4408 switch(normal.getFunctionSpace().getTypeCode()) {
1159     case Nodes:
1160     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1161     break;
1162     case ReducedNodes:
1163     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1164     break;
1165     case Elements:
1166     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1167     break;
1168     case ReducedElements:
1169     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1170     break;
1171     case FaceElements:
1172     case ReducedFaceElements:
1173 caltinay 4496 Assemble_getNormal(mesh->Nodes, mesh->FaceElements, normal);
1174 caltinay 4408 break;
1175     case Points:
1176     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1177     break;
1178     case ContactElementsOne:
1179     case ContactElementsZero:
1180     case ReducedContactElementsOne:
1181     case ReducedContactElementsZero:
1182 caltinay 4496 Assemble_getNormal(mesh->Nodes, mesh->ContactElements, normal);
1183 caltinay 4408 break;
1184     case DegreesOfFreedom:
1185     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1186     break;
1187     case ReducedDegreesOfFreedom:
1188     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1189     break;
1190     default:
1191     stringstream temp;
1192     temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1193     throw FinleyAdapterException(temp.str());
1194     break;
1195     }
1196     checkFinleyError();
1197 jgs 82 }
1198 jgs 149
1199 jgs 82 //
1200 caltinay 4408 // interpolates data to other domain
1201 jgs 82 //
1202 caltinay 4441 void MeshAdapter::interpolateACross(escript::Data& target, const escript::Data& source) const
1203 jgs 82 {
1204 caltinay 4441 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains.");
1205 jgs 82 }
1206 jgs 149
1207 jgs 82 //
1208 caltinay 4408 // calculates the integral of a function defined on arg
1209 jgs 82 //
1210 caltinay 4408 void MeshAdapter::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
1211 jgs 82 {
1212 caltinay 4408 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1213     if (argDomain!=*this)
1214     throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1215 jgs 82
1216 caltinay 4408 double blocktimer_start = blocktimer_time();
1217 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
1218 caltinay 4408 switch(arg.getFunctionSpace().getTypeCode()) {
1219     case Nodes:
1220 caltinay 4441 {
1221     escript::Data temp(arg, escript::function(*this));
1222     Assemble_integrate(mesh->Nodes, mesh->Elements, temp, &integrals[0]);
1223     }
1224 caltinay 4408 break;
1225     case ReducedNodes:
1226 caltinay 4441 {
1227     escript::Data temp(arg, escript::function(*this));
1228     Assemble_integrate(mesh->Nodes, mesh->Elements, temp, &integrals[0]);
1229     }
1230 caltinay 4408 break;
1231     case Elements:
1232     case ReducedElements:
1233 caltinay 4441 Assemble_integrate(mesh->Nodes, mesh->Elements, arg, &integrals[0]);
1234 caltinay 4408 break;
1235     case FaceElements:
1236     case ReducedFaceElements:
1237 caltinay 4441 Assemble_integrate(mesh->Nodes, mesh->FaceElements, arg, &integrals[0]);
1238 caltinay 4408 break;
1239     case Points:
1240     throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1241     break;
1242     case ContactElementsZero:
1243     case ReducedContactElementsZero:
1244     case ContactElementsOne:
1245     case ReducedContactElementsOne:
1246 caltinay 4441 Assemble_integrate(mesh->Nodes, mesh->ContactElements, arg, &integrals[0]);
1247 caltinay 4408 break;
1248     case DegreesOfFreedom:
1249     case ReducedDegreesOfFreedom:
1250 caltinay 4441 {
1251     escript::Data temp(arg, escript::function(*this));
1252     Assemble_integrate(mesh->Nodes, mesh->Elements, temp, &integrals[0]);
1253     }
1254 caltinay 4408 break;
1255     default:
1256     stringstream temp;
1257     temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1258     throw FinleyAdapterException(temp.str());
1259     break;
1260     }
1261     checkFinleyError();
1262     blocktimer_increment("integrate()", blocktimer_start);
1263 jgs 82 }
1264 jgs 149
1265 jgs 82 //
1266 caltinay 4408 // calculates the gradient of arg
1267 jgs 82 //
1268 caltinay 4408 void MeshAdapter::setToGradient(escript::Data& grad, const escript::Data& arg) const
1269 jgs 82 {
1270 caltinay 4408 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1271     if (argDomain!=*this)
1272     throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1273     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1274     if (gradDomain!=*this)
1275     throw FinleyAdapterException("Error - Illegal domain of gradient");
1276 jgs 82
1277 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
1278 caltinay 4441 escript::Data nodeData;
1279 caltinay 4408 if (getMPISize()>1) {
1280     if (arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom) {
1281 caltinay 4441 nodeData=escript::Data(arg, continuousFunction(*this));
1282 caltinay 4408 } else if(arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom) {
1283 caltinay 4441 nodeData=escript::Data(arg, reducedContinuousFunction(*this));
1284 caltinay 4408 } else {
1285 caltinay 4441 nodeData = arg;
1286 caltinay 4408 }
1287     } else {
1288 caltinay 4441 nodeData = arg;
1289 caltinay 4408 }
1290     switch(grad.getFunctionSpace().getTypeCode()) {
1291     case Nodes:
1292     throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1293     break;
1294     case ReducedNodes:
1295     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1296     break;
1297     case Elements:
1298     case ReducedElements:
1299 caltinay 4441 Assemble_gradient(mesh->Nodes, mesh->Elements, grad, nodeData);
1300 caltinay 4408 break;
1301     case FaceElements:
1302     case ReducedFaceElements:
1303 caltinay 4441 Assemble_gradient(mesh->Nodes, mesh->FaceElements, grad, nodeData);
1304 caltinay 4408 break;
1305     case Points:
1306     throw FinleyAdapterException("Error - Gradient at points is not supported.");
1307     break;
1308     case ContactElementsZero:
1309     case ReducedContactElementsZero:
1310     case ContactElementsOne:
1311     case ReducedContactElementsOne:
1312 caltinay 4441 Assemble_gradient(mesh->Nodes, mesh->ContactElements, grad, nodeData);
1313 caltinay 4408 break;
1314     case DegreesOfFreedom:
1315     throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1316     break;
1317     case ReducedDegreesOfFreedom:
1318     throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1319     break;
1320     default:
1321     stringstream temp;
1322     temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1323     throw FinleyAdapterException(temp.str());
1324     break;
1325     }
1326     checkFinleyError();
1327 jgs 82 }
1328 jgs 149
1329 jgs 82 //
1330 caltinay 4408 // returns the size of elements
1331 jgs 82 //
1332 jgs 480 void MeshAdapter::setToSize(escript::Data& size) const
1333 jgs 82 {
1334 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
1335 caltinay 4408 switch(size.getFunctionSpace().getTypeCode()) {
1336     case Nodes:
1337     throw FinleyAdapterException("Error - Size of nodes is not supported.");
1338     break;
1339     case ReducedNodes:
1340     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1341     break;
1342     case Elements:
1343     case ReducedElements:
1344 caltinay 4441 Assemble_getSize(mesh->Nodes, mesh->Elements, size);
1345 caltinay 4408 break;
1346     case FaceElements:
1347     case ReducedFaceElements:
1348 caltinay 4441 Assemble_getSize(mesh->Nodes, mesh->FaceElements, size);
1349 caltinay 4408 break;
1350     case Points:
1351     throw FinleyAdapterException("Error - Size of point elements is not supported.");
1352     break;
1353     case ContactElementsZero:
1354     case ContactElementsOne:
1355     case ReducedContactElementsZero:
1356     case ReducedContactElementsOne:
1357 caltinay 4441 Assemble_getSize(mesh->Nodes,mesh->ContactElements,size);
1358 caltinay 4408 break;
1359     case DegreesOfFreedom:
1360     throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1361     break;
1362     case ReducedDegreesOfFreedom:
1363     throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1364     break;
1365     default:
1366     stringstream temp;
1367     temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1368     throw FinleyAdapterException(temp.str());
1369     break;
1370     }
1371     checkFinleyError();
1372 jgs 82 }
1373 jgs 149
1374 caltinay 2151 //
1375     // sets the location of nodes
1376     //
1377 jgs 480 void MeshAdapter::setNewX(const escript::Data& new_x)
1378 jgs 82 {
1379 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
1380 caltinay 4408 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1381     if (newDomain!=*this)
1382     throw FinleyAdapterException("Error - Illegal domain of new point locations");
1383     if (new_x.getFunctionSpace() == continuousFunction(*this)) {
1384 caltinay 4496 mesh->setCoordinates(new_x);
1385 caltinay 4408 } else {
1386     throw FinleyAdapterException("As of escript version 3.3 SetX() only accepts ContinuousFunction arguments. Please interpolate.");
1387     }
1388     checkFinleyError();
1389 jgs 82 }
1390 jgs 149
1391 jfenwick 2642 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1392     {
1393 caltinay 3993 if (getMPISize() > 1) {
1394 jfenwick 3259 #ifdef ESYS_MPI
1395 caltinay 4441 int myFirstNode=0, myLastNode=0, k=0;
1396     int* globalNodeIndex=0;
1397 caltinay 4496 Mesh* mesh_p=m_finleyMesh.get();
1398 caltinay 3998 /*
1399     * this method is only used by saveDataCSV which would use the returned
1400     * values for reduced nodes wrongly so this case is disabled for now
1401 caltinay 4408 if (fs_code == FINLEY_REDUCED_NODES) {
1402 caltinay 4496 myFirstNode = NodeFile_getFirstReducedNode(mesh_p->Nodes);
1403     myLastNode = NodeFile_getLastReducedNode(mesh_p->Nodes);
1404     globalNodeIndex = NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1405 caltinay 4408 } else
1406 caltinay 3998 */
1407 caltinay 4408 if (fs_code == FINLEY_NODES) {
1408 caltinay 4428 myFirstNode = mesh_p->Nodes->getFirstNode();
1409     myLastNode = mesh_p->Nodes->getLastNode();
1410     globalNodeIndex = mesh_p->Nodes->borrowGlobalNodesIndex();
1411 caltinay 4408 } else {
1412 caltinay 3993 throw FinleyAdapterException("Unsupported function space type for ownSample()");
1413     }
1414    
1415     k=globalNodeIndex[id];
1416 caltinay 4441 return ((myFirstNode <= k) && (k < myLastNode));
1417 caltinay 3993 #endif
1418 jfenwick 2642 }
1419 jfenwick 2644 return true;
1420 jfenwick 2642 }
1421    
1422    
1423 caltinay 2151 //
1424     // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1425     //
1426 caltinay 4408 escript::ASM_ptr MeshAdapter::newSystemMatrix(const int row_blocksize,
1427     const escript::FunctionSpace& row_functionspace,
1428     const int column_blocksize,
1429     const escript::FunctionSpace& column_functionspace,
1430     const int type) const
1431 jgs 82 {
1432 caltinay 4408 // is the domain right?
1433     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1434     if (row_domain!=*this)
1435     throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1436     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1437     if (col_domain!=*this)
1438     throw FinleyAdapterException("Error - domain of column function space does not match the domain of matrix generator.");
1439    
1440     int reduceRowOrder=0;
1441     int reduceColOrder=0;
1442     // is the function space type right?
1443     if (row_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1444     reduceRowOrder=1;
1445     } else if (row_functionspace.getTypeCode() != DegreesOfFreedom) {
1446     throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1447     }
1448     if (column_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1449     reduceColOrder=1;
1450     } else if (column_functionspace.getTypeCode() != DegreesOfFreedom) {
1451     throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1452     }
1453    
1454     // generate matrix:
1455 caltinay 4496 Paso_SystemMatrixPattern* fsystemMatrixPattern=
1456     getFinley_Mesh()->getPattern(reduceRowOrder, reduceColOrder);
1457 caltinay 4408 checkFinleyError();
1458     Paso_SystemMatrix* fsystemMatrix;
1459     const int trilinos = 0;
1460     if (trilinos) {
1461 ksteube 1312 #ifdef TRILINOS
1462 caltinay 4408 // FIXME: Allocation Epetra_VrbMatrix here...
1463 ksteube 1312 #endif
1464 caltinay 4408 } else {
1465     fsystemMatrix=Paso_SystemMatrix_alloc(type, fsystemMatrixPattern,
1466     row_blocksize, column_blocksize, FALSE);
1467     }
1468     checkPasoError();
1469     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1470     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1471     return escript::ASM_ptr(sma);
1472 jgs 82 }
1473 caltinay 2151
1474     //
1475 gross 1366 // creates a TransportProblemAdapter
1476 caltinay 2151 //
1477 caltinay 4408 escript::ATP_ptr MeshAdapter::newTransportProblem(const int blocksize,
1478     const escript::FunctionSpace& functionspace, const int type) const
1479 gross 1366 {
1480 caltinay 4408 // is the domain right?
1481     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1482     if (domain!=*this)
1483     throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1484    
1485     // is the function space type right?
1486     int reduceOrder=0;
1487     if (functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1488     reduceOrder=1;
1489     } else if (functionspace.getTypeCode() != DegreesOfFreedom) {
1490     throw FinleyAdapterException("Error - illegal function space type for transport problem.");
1491     }
1492    
1493     // generate transport problem:
1494 caltinay 4496 Paso_SystemMatrixPattern* fsystemMatrixPattern=
1495     getFinley_Mesh()->getPattern(reduceOrder, reduceOrder);
1496 caltinay 4408 checkFinleyError();
1497     Paso_TransportProblem* transportProblem;
1498     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern, blocksize);
1499     checkPasoError();
1500     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1501     TransportProblemAdapter* tpa=new TransportProblemAdapter(
1502     transportProblem, blocksize, functionspace);
1503     return escript::ATP_ptr(tpa);
1504 gross 1366 }
1505 jgs 149
1506 jgs 82 //
1507 caltinay 4408 // returns true if data on functionSpaceCode is considered as being cell centered
1508 jgs 82 //
1509     bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1510     {
1511 caltinay 4408 switch(functionSpaceCode) {
1512     case Nodes:
1513     case DegreesOfFreedom:
1514     case ReducedDegreesOfFreedom:
1515     return false;
1516     case Elements:
1517     case FaceElements:
1518     case Points:
1519     case ContactElementsZero:
1520     case ContactElementsOne:
1521     case ReducedElements:
1522     case ReducedFaceElements:
1523     case ReducedContactElementsZero:
1524     case ReducedContactElementsOne:
1525     return true;
1526     default:
1527     stringstream temp;
1528     temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1529     throw FinleyAdapterException(temp.str());
1530     break;
1531     }
1532     return false;
1533 jgs 82 }
1534 jgs 149
1535 jfenwick 2635 bool
1536 caltinay 2778 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1537 jfenwick 2635 {
1538 caltinay 4408 /* The idea is to use equivalence classes. [Types which can be interpolated
1539     back and forth]:
1540     class 1: DOF <-> Nodes
1541     class 2: ReducedDOF <-> ReducedNodes
1542     class 3: Points
1543     class 4: Elements
1544     class 5: ReducedElements
1545     class 6: FaceElements
1546     class 7: ReducedFaceElements
1547     class 8: ContactElementZero <-> ContactElementOne
1548     class 9: ReducedContactElementZero <-> ReducedContactElementOne
1549 jfenwick 2635
1550 caltinay 4408 There is also a set of lines. Interpolation is possible down a line but
1551     not between lines.
1552     class 1 and 2 belong to all lines so aren't considered.
1553     line 0: class 3
1554     line 1: class 4,5
1555     line 2: class 6,7
1556     line 3: class 8,9
1557 jfenwick 2635
1558 caltinay 4408 For classes with multiple members (eg class 2) we have vars to record if
1559     there is at least one instance.
1560     e.g. hasnodes is true if we have at least one instance of Nodes.
1561     */
1562 caltinay 2940 if (fs.empty())
1563 caltinay 2778 return false;
1564 caltinay 4408
1565 caltinay 2778 vector<int> hasclass(10);
1566 caltinay 4408 vector<int> hasline(4);
1567 jfenwick 2635 bool hasnodes=false;
1568     bool hasrednodes=false;
1569     bool hascez=false;
1570     bool hasrcez=false;
1571 caltinay 4408 for (int i=0;i<fs.size();++i) {
1572     switch(fs[i]) {
1573     case Nodes:
1574     hasnodes=true; // fall through
1575     case DegreesOfFreedom:
1576     hasclass[1]=1;
1577     break;
1578     case ReducedNodes:
1579     hasrednodes=true; // fall through
1580     case ReducedDegreesOfFreedom:
1581     hasclass[2]=1;
1582     break;
1583     case Points:
1584     hasline[0]=1;
1585     hasclass[3]=1;
1586     break;
1587     case Elements:
1588     hasclass[4]=1;
1589     hasline[1]=1;
1590     break;
1591     case ReducedElements:
1592     hasclass[5]=1;
1593     hasline[1]=1;
1594     break;
1595     case FaceElements:
1596     hasclass[6]=1;
1597     hasline[2]=1;
1598     break;
1599     case ReducedFaceElements:
1600     hasclass[7]=1;
1601     hasline[2]=1;
1602     break;
1603     case ContactElementsZero:
1604     hascez=true; // fall through
1605     case ContactElementsOne:
1606     hasclass[8]=1;
1607     hasline[3]=1;
1608     break;
1609     case ReducedContactElementsZero:
1610     hasrcez=true; // fall through
1611     case ReducedContactElementsOne:
1612     hasclass[9]=1;
1613     hasline[3]=1;
1614     break;
1615     default:
1616     return false;
1617     }
1618 jfenwick 2635 }
1619     int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1620 caltinay 4408
1621 jfenwick 2635 // fail if we have more than one leaf group
1622     if (totlines>1)
1623 caltinay 4408 return false; // there are at least two branches we can't interpolate between
1624     else if (totlines==1) {
1625     if (hasline[0]==1) // we have points
1626     resultcode=Points;
1627     else if (hasline[1]==1) {
1628     if (hasclass[5]==1)
1629     resultcode=ReducedElements;
1630     else
1631     resultcode=Elements;
1632     } else if (hasline[2]==1) {
1633     if (hasclass[7]==1)
1634     resultcode=ReducedFaceElements;
1635     else
1636     resultcode=FaceElements;
1637     } else { // so we must be in line3
1638     if (hasclass[9]==1) {
1639     // need something from class 9
1640     resultcode=(hasrcez ? ReducedContactElementsZero : ReducedContactElementsOne);
1641     } else {
1642     // something from class 8
1643     resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1644     }
1645     }
1646     } else { // totlines==0
1647     if (hasclass[2]==1) {
1648     // something from class 2
1649     resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
1650     } else {
1651     // something from class 1
1652     resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
1653     }
1654 jfenwick 2635 }
1655     return true;
1656     }
1657    
1658 caltinay 4408 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,
1659     int functionSpaceType_target) const
1660 jgs 82 {
1661 caltinay 4408 switch(functionSpaceType_source) {
1662     case Nodes:
1663     switch(functionSpaceType_target) {
1664     case Nodes:
1665     case ReducedNodes:
1666     case ReducedDegreesOfFreedom:
1667     case DegreesOfFreedom:
1668     case Elements:
1669     case ReducedElements:
1670     case FaceElements:
1671     case ReducedFaceElements:
1672     case Points:
1673     case ContactElementsZero:
1674     case ReducedContactElementsZero:
1675     case ContactElementsOne:
1676     case ReducedContactElementsOne:
1677     return true;
1678     default:
1679     stringstream temp;
1680     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1681     throw FinleyAdapterException(temp.str());
1682     }
1683     break;
1684     case ReducedNodes:
1685     switch(functionSpaceType_target) {
1686     case ReducedNodes:
1687     case ReducedDegreesOfFreedom:
1688     case Elements:
1689     case ReducedElements:
1690     case FaceElements:
1691     case ReducedFaceElements:
1692     case Points:
1693     case ContactElementsZero:
1694     case ReducedContactElementsZero:
1695     case ContactElementsOne:
1696     case ReducedContactElementsOne:
1697     return true;
1698     case Nodes:
1699     case DegreesOfFreedom:
1700     return false;
1701     default:
1702     stringstream temp;
1703     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1704     throw FinleyAdapterException(temp.str());
1705     }
1706     break;
1707     case Elements:
1708     if (functionSpaceType_target==Elements) {
1709     return true;
1710     } else if (functionSpaceType_target==ReducedElements) {
1711     return true;
1712     } else {
1713     return false;
1714     }
1715     case ReducedElements:
1716     if (functionSpaceType_target==ReducedElements) {
1717     return true;
1718     } else {
1719     return false;
1720     }
1721     case FaceElements:
1722     if (functionSpaceType_target==FaceElements) {
1723     return true;
1724     } else if (functionSpaceType_target==ReducedFaceElements) {
1725     return true;
1726     } else {
1727     return false;
1728     }
1729     case ReducedFaceElements:
1730     if (functionSpaceType_target==ReducedFaceElements) {
1731     return true;
1732     } else {
1733     return false;
1734     }
1735     case Points:
1736     if (functionSpaceType_target==Points) {
1737     return true;
1738     } else {
1739     return false;
1740     }
1741     case ContactElementsZero:
1742     case ContactElementsOne:
1743     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1744     return true;
1745     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1746     return true;
1747     } else {
1748     return false;
1749     }
1750     case ReducedContactElementsZero:
1751     case ReducedContactElementsOne:
1752     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1753     return true;
1754     } else {
1755     return false;
1756     }
1757     case DegreesOfFreedom:
1758     switch(functionSpaceType_target) {
1759     case ReducedDegreesOfFreedom:
1760     case DegreesOfFreedom:
1761     case Nodes:
1762     case ReducedNodes:
1763     case Elements:
1764     case ReducedElements:
1765     case Points:
1766     case FaceElements:
1767     case ReducedFaceElements:
1768     case ContactElementsZero:
1769     case ReducedContactElementsZero:
1770     case ContactElementsOne:
1771     case ReducedContactElementsOne:
1772     return true;
1773     default:
1774     stringstream temp;
1775     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1776     throw FinleyAdapterException(temp.str());
1777     }
1778     break;
1779     case ReducedDegreesOfFreedom:
1780     switch(functionSpaceType_target) {
1781     case ReducedDegreesOfFreedom:
1782     case ReducedNodes:
1783     case Elements:
1784     case ReducedElements:
1785     case FaceElements:
1786     case ReducedFaceElements:
1787     case Points:
1788     case ContactElementsZero:
1789     case ReducedContactElementsZero:
1790     case ContactElementsOne:
1791     case ReducedContactElementsOne:
1792     return true;
1793     case Nodes:
1794     case DegreesOfFreedom:
1795     return false;
1796     default:
1797     stringstream temp;
1798     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1799     throw FinleyAdapterException(temp.str());
1800     }
1801     break;
1802     default:
1803     stringstream temp;
1804     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1805     throw FinleyAdapterException(temp.str());
1806     break;
1807     }
1808     return false;
1809 jgs 82 }
1810 jgs 149
1811 caltinay 4408 signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source, int functionSpaceType_target) const
1812 jfenwick 4255 {
1813     if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1814     return 1;
1815 caltinay 4408
1816     if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1817 jfenwick 4255 return -1;
1818 caltinay 4408
1819 jfenwick 4255 return 0;
1820 caltinay 4408 }
1821    
1822     bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,
1823     const escript::AbstractDomain& targetDomain,
1824     int functionSpaceType_target) const
1825 jgs 82 {
1826 caltinay 4408 return false;
1827 jgs 82 }
1828 jgs 104
1829 caltinay 4408 bool MeshAdapter::operator==(const escript::AbstractDomain& other) const
1830 jgs 82 {
1831 caltinay 4408 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1832     if (temp) {
1833     return (m_finleyMesh==temp->m_finleyMesh);
1834     } else {
1835     return false;
1836     }
1837 jgs 82 }
1838    
1839 caltinay 4408 bool MeshAdapter::operator!=(const escript::AbstractDomain& other) const
1840 jgs 104 {
1841 caltinay 4408 return !(operator==(other));
1842 jgs 104 }
1843    
1844 gross 1859 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1845 jgs 102 {
1846 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
1847 caltinay 4408 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
1848     package, symmetry, mesh->MPIInfo);
1849 jgs 102 }
1850 caltinay 3681
1851 gross 1859 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1852     {
1853 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
1854 caltinay 4408 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
1855     package, symmetry, mesh->MPIInfo);
1856 gross 1859 }
1857 jgs 149
1858 jgs 480 escript::Data MeshAdapter::getX() const
1859 jgs 102 {
1860 caltinay 4408 return continuousFunction(*this).getX();
1861 jgs 102 }
1862 jgs 149
1863 jgs 480 escript::Data MeshAdapter::getNormal() const
1864 jgs 102 {
1865 caltinay 4408 return functionOnBoundary(*this).getNormal();
1866 jgs 102 }
1867 jgs 149
1868 jgs 480 escript::Data MeshAdapter::getSize() const
1869 jgs 102 {
1870 caltinay 4408 return escript::function(*this).getSize();
1871 jgs 102 }
1872    
1873 jfenwick 2487 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1874 gross 1062 {
1875 caltinay 4408 int *out = NULL;
1876 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
1877 caltinay 4408 switch (functionSpaceType) {
1878     case Nodes:
1879     out=mesh->Nodes->Id;
1880     break;
1881     case ReducedNodes:
1882     out=mesh->Nodes->reducedNodesId;
1883     break;
1884     case Elements:
1885     case ReducedElements:
1886     out=mesh->Elements->Id;
1887     break;
1888     case FaceElements:
1889     case ReducedFaceElements:
1890     out=mesh->FaceElements->Id;
1891     break;
1892     case Points:
1893     out=mesh->Points->Id;
1894     break;
1895     case ContactElementsZero:
1896     case ReducedContactElementsZero:
1897     case ContactElementsOne:
1898     case ReducedContactElementsOne:
1899     out=mesh->ContactElements->Id;
1900     break;
1901     case DegreesOfFreedom:
1902     out=mesh->Nodes->degreesOfFreedomId;
1903     break;
1904     case ReducedDegreesOfFreedom:
1905     out=mesh->Nodes->reducedDegreesOfFreedomId;
1906     break;
1907     default:
1908     stringstream temp;
1909     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1910     throw FinleyAdapterException(temp.str());
1911     break;
1912     }
1913     return out;
1914 gross 1062 }
1915 jgs 110 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1916     {
1917 caltinay 4408 int out=0;
1918 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
1919 caltinay 4408 switch (functionSpaceType) {
1920     case Nodes:
1921     out=mesh->Nodes->Tag[sampleNo];
1922     break;
1923     case ReducedNodes:
1924     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1925     break;
1926     case Elements:
1927     case ReducedElements:
1928     out=mesh->Elements->Tag[sampleNo];
1929     break;
1930     case FaceElements:
1931     case ReducedFaceElements:
1932     out=mesh->FaceElements->Tag[sampleNo];
1933     break;
1934     case Points:
1935     out=mesh->Points->Tag[sampleNo];
1936     break;
1937     case ContactElementsZero:
1938     case ReducedContactElementsZero:
1939     case ContactElementsOne:
1940     case ReducedContactElementsOne:
1941     out=mesh->ContactElements->Tag[sampleNo];
1942     break;
1943     case DegreesOfFreedom:
1944     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1945     break;
1946     case ReducedDegreesOfFreedom:
1947     throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1948     break;
1949     default:
1950     stringstream temp;
1951     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1952     throw FinleyAdapterException(temp.str());
1953     break;
1954     }
1955     return out;
1956 jgs 110 }
1957    
1958 gross 1062
1959 gross 767 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1960     {
1961 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
1962 caltinay 4408 switch(functionSpaceType) {
1963     case Nodes:
1964 caltinay 4428 mesh->Nodes->setTags(newTag, mask);
1965 caltinay 4408 break;
1966     case ReducedNodes:
1967     throw FinleyAdapterException("Error - ReducedNodes does not support tags");
1968     case DegreesOfFreedom:
1969     throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1970     case ReducedDegreesOfFreedom:
1971     throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1972     case Elements:
1973     case ReducedElements:
1974 caltinay 4431 mesh->Elements->setTags(newTag, mask);
1975 caltinay 4408 break;
1976     case FaceElements:
1977     case ReducedFaceElements:
1978 caltinay 4431 mesh->FaceElements->setTags(newTag, mask);
1979 caltinay 4408 break;
1980     case Points:
1981 caltinay 4431 mesh->Points->setTags(newTag, mask);
1982 caltinay 4408 break;
1983     case ContactElementsZero:
1984     case ReducedContactElementsZero:
1985     case ContactElementsOne:
1986     case ReducedContactElementsOne:
1987 caltinay 4431 mesh->ContactElements->setTags(newTag, mask);
1988 caltinay 4408 break;
1989     default:
1990     stringstream temp;
1991     temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1992     throw FinleyAdapterException(temp.str());
1993     }
1994     checkFinleyError();
1995 gross 767 }
1996    
1997 caltinay 4408 void MeshAdapter::setTagMap(const string& name, int tag)
1998 gross 1044 {
1999 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
2000     mesh->addTagMap(name.c_str(), tag);
2001 caltinay 4408 checkFinleyError();
2002 gross 1044 }
2003 gross 767
2004 caltinay 2778 int MeshAdapter::getTag(const string& name) const
2005 gross 1044 {
2006 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
2007     int tag = mesh->getTag(name.c_str());
2008 caltinay 4408 checkFinleyError();
2009     return tag;
2010 gross 1044 }
2011    
2012 caltinay 2778 bool MeshAdapter::isValidTagName(const string& name) const
2013 gross 1044 {
2014 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
2015     return mesh->isValidTagName(name.c_str());
2016 gross 1044 }
2017    
2018 caltinay 2778 string MeshAdapter::showTagNames() const
2019 gross 1044 {
2020 caltinay 4408 stringstream temp;
2021 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
2022 caltinay 4437 TagMap::const_iterator it = mesh->tagMap.begin();
2023     while (it != mesh->tagMap.end()) {
2024     temp << it->first;
2025     ++it;
2026     if (it != mesh->tagMap.end())
2027     temp << ", ";
2028 caltinay 4408 }
2029     return temp.str();
2030 gross 1044 }
2031    
2032 gross 1716 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2033     {
2034 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
2035 caltinay 4408 switch(functionSpaceCode) {
2036     case Nodes:
2037 caltinay 4428 return mesh->Nodes->tagsInUse.size();
2038 caltinay 4408 case ReducedNodes:
2039     throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2040     case DegreesOfFreedom:
2041     throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2042     case ReducedDegreesOfFreedom:
2043     throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2044     case Elements:
2045     case ReducedElements:
2046 caltinay 4428 return mesh->Elements->tagsInUse.size();
2047 caltinay 4408 case FaceElements:
2048     case ReducedFaceElements:
2049 caltinay 4428 return mesh->FaceElements->tagsInUse.size();
2050 caltinay 4408 case Points:
2051 caltinay 4428 return mesh->Points->tagsInUse.size();
2052 caltinay 4408 case ContactElementsZero:
2053     case ReducedContactElementsZero:
2054     case ContactElementsOne:
2055     case ReducedContactElementsOne:
2056 caltinay 4428 return mesh->ContactElements->tagsInUse.size();
2057 caltinay 4408 default:
2058 caltinay 4428 stringstream ss;
2059     ss << "Finley does not know anything about function space type "
2060     << functionSpaceCode;
2061     throw FinleyAdapterException(ss.str());
2062 caltinay 4408 }
2063 caltinay 4428 return 0;
2064 gross 1716 }
2065 jfenwick 2487
2066     const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2067 gross 1716 {
2068 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
2069 caltinay 4408 switch(functionSpaceCode) {
2070     case Nodes:
2071 caltinay 4428 return &mesh->Nodes->tagsInUse[0];
2072 caltinay 4408 case ReducedNodes:
2073     throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2074     case DegreesOfFreedom:
2075     throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2076     case ReducedDegreesOfFreedom:
2077     throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2078     case Elements:
2079     case ReducedElements:
2080 caltinay 4428 return &mesh->Elements->tagsInUse[0];
2081 caltinay 4408 case FaceElements:
2082     case ReducedFaceElements:
2083 caltinay 4428 return &mesh->FaceElements->tagsInUse[0];
2084 caltinay 4408 case Points:
2085 caltinay 4428 return &mesh->Points->tagsInUse[0];
2086 caltinay 4408 case ContactElementsZero:
2087     case ReducedContactElementsZero:
2088     case ContactElementsOne:
2089     case ReducedContactElementsOne:
2090 caltinay 4428 return &mesh->ContactElements->tagsInUse[0];
2091 caltinay 4408 default:
2092     stringstream temp;
2093     temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2094     throw FinleyAdapterException(temp.str());
2095     }
2096 caltinay 4428 return NULL;
2097 gross 1716 }
2098    
2099    
2100 jfenwick 1802 bool MeshAdapter::canTag(int functionSpaceCode) const
2101     {
2102 caltinay 4408 switch(functionSpaceCode) {
2103     case Nodes:
2104     case Elements:
2105     case ReducedElements:
2106     case FaceElements:
2107     case ReducedFaceElements:
2108     case Points:
2109     case ContactElementsZero:
2110     case ReducedContactElementsZero:
2111     case ContactElementsOne:
2112     case ReducedContactElementsOne:
2113     return true;
2114     default:
2115     return false;
2116     }
2117 jfenwick 1802 }
2118    
2119 caltinay 4408 escript::AbstractDomain::StatusType MeshAdapter::getStatus() const
2120 gross 2533 {
2121 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
2122     return mesh->getStatus();
2123 gross 2533 }
2124 jfenwick 1802
2125 gross 2856 int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2126     {
2127 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
2128 caltinay 4408 int order =-1;
2129     switch(functionSpaceCode) {
2130     case Nodes:
2131     case DegreesOfFreedom:
2132     order=mesh->approximationOrder;
2133     break;
2134     case ReducedNodes:
2135     case ReducedDegreesOfFreedom:
2136     order=mesh->reducedApproximationOrder;
2137     break;
2138     case Elements:
2139     case FaceElements:
2140     case Points:
2141     case ContactElementsZero:
2142     case ContactElementsOne:
2143     order=mesh->integrationOrder;
2144     break;
2145     case ReducedElements:
2146     case ReducedFaceElements:
2147     case ReducedContactElementsZero:
2148     case ReducedContactElementsOne:
2149     order=mesh->reducedIntegrationOrder;
2150     break;
2151     default:
2152     stringstream temp;
2153     temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2154     throw FinleyAdapterException(temp.str());
2155     }
2156     return order;
2157 gross 2856 }
2158 gross 2533
2159 jfenwick 3259 bool MeshAdapter::supportsContactElements() const
2160 caltinay 2898 {
2161 caltinay 4408 return true;
2162 jfenwick 3259 }
2163    
2164 caltinay 4408 void MeshAdapter::addDiracPoints(const vector<double>& points,
2165     const vector<int>& tags) const
2166 gross 3515 {
2167 caltinay 4408 // points will be flattened
2168     const int dim = getDim();
2169     int numPoints=points.size()/dim;
2170     int numTags=tags.size();
2171 caltinay 4496 Mesh* mesh=m_finleyMesh.get();
2172 caltinay 4408
2173     if ( points.size() % dim != 0 ) {
2174     throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2175     }
2176    
2177     if ((numTags > 0) && (numPoints != numTags))
2178     throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2179    
2180 caltinay 4496 mesh->addPoints(numPoints, &points[0</