/[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 3227 - (hide annotations)
Thu Sep 30 06:07:08 2010 UTC (8 years, 6 months ago) by jfenwick
File size: 72277 byte(s)
Pass1 or moving MPI stuff out of paso

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