/[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 4657 - (hide annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years, 1 month ago) by jfenwick
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 90691 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



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