/[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 3086 - (hide annotations)
Thu Aug 5 05:07:58 2010 UTC (8 years, 8 months ago) by jfenwick
File size: 85811 byte(s)
Another pass at removing finley

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