/[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 3090 - (hide annotations)
Wed Aug 11 00:51:28 2010 UTC (8 years, 8 months ago) by jfenwick
File size: 72859 byte(s)
Removed:
DUDLEY_CONTACT_ELEMENTS_1
DUDLEY_CONTACT_ELEMENTS_2
DUDLEY_REDUCED_CONTACT_ELEMENTS_1
DUDLEY_REDUCED_CONTACT_ELEMENTS_2

escript tests now query if the domain supports contact elements
before trying to use them.


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