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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4114 - (hide annotations)
Fri Dec 14 04:24:46 2012 UTC (6 years, 9 months ago) by caltinay
File size: 71448 byte(s)
Time to remove deprecated saveVTK/DX methods from Data and Domain.

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