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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5136 - (hide annotations)
Tue Sep 9 07:13:55 2014 UTC (4 years, 6 months ago) by caltinay
File size: 72466 byte(s)
ripley now supports paso solvers again and returns an appropriate matrix type
id. Changed the getSystemMatrixTypeId() method to take a full SolverBuddy
instance and made some other simplifications.

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26