/[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 3242 - (hide annotations)
Tue Oct 5 06:08:28 2010 UTC (8 years, 8 months ago) by jfenwick
File size: 71783 byte(s)
Remove references to asContinuousDomain

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