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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3114 - (hide annotations)
Fri Aug 27 05:26:25 2010 UTC (8 years, 7 months ago) by jfenwick
File size: 72668 byte(s)
It doesn't pass all tests but this is major progress

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