/[escript]/trunk/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Annotation of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3355 - (hide annotations)
Tue Nov 16 06:35:06 2010 UTC (8 years, 6 months ago) by caltinay
File size: 85795 byte(s)
Fixed building with boost < 1.40 by
- removing the saveVTK method from dudley's MeshAdapter completely
  (no point introducing it now when we are trying to get rid of it soon)
- adding a helper function to weipa's python layer which can be called from C++
  with older boost versions.


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 1312 extern "C" {
22 phornby 2078 #include "esysUtils/blocktimer.h"
23 ksteube 1312 }
24 jgs 480
25 caltinay 3344 #include <boost/python/import.hpp>
26    
27 jgs 82 using namespace std;
28     using namespace escript;
29    
30     namespace finley {
31 jgs 149
32 jgs 82 //
33 jgs 149 // define the static constants
34 jgs 82 MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
35     const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;
36     const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
37     const int MeshAdapter::Nodes=FINLEY_NODES;
38 gross 1062 const int MeshAdapter::ReducedNodes=FINLEY_REDUCED_NODES;
39 jgs 82 const int MeshAdapter::Elements=FINLEY_ELEMENTS;
40 gross 1059 const int MeshAdapter::ReducedElements=FINLEY_REDUCED_ELEMENTS;
41 jgs 82 const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;
42 gross 1059 const int MeshAdapter::ReducedFaceElements=FINLEY_REDUCED_FACE_ELEMENTS;
43 jgs 82 const int MeshAdapter::Points=FINLEY_POINTS;
44     const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;
45 gross 1059 const int MeshAdapter::ReducedContactElementsZero=FINLEY_REDUCED_CONTACT_ELEMENTS_1;
46 jgs 82 const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
47 gross 1059 const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;
48 jgs 82
49     MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
50     {
51 phornby 1628 setFunctionSpaceTypeNames();
52     //
53     // need to use a null_deleter as Finley_Mesh_free deletes the pointer
54     // for us.
55     m_finleyMesh.reset(finleyMesh,null_deleter());
56 jgs 82 }
57 jgs 149
58 jgs 82 //
59     // The copy constructor should just increment the use count
60     MeshAdapter::MeshAdapter(const MeshAdapter& in):
61     m_finleyMesh(in.m_finleyMesh)
62     {
63 phornby 1628 setFunctionSpaceTypeNames();
64 jgs 82 }
65    
66     MeshAdapter::~MeshAdapter()
67     {
68 phornby 1628 //
69     // I hope the case for the pointer being zero has been taken care of.
70     // cout << "In MeshAdapter destructor." << endl;
71     if (m_finleyMesh.unique()) {
72     Finley_Mesh_free(m_finleyMesh.get());
73     }
74 jgs 82 }
75    
76 ksteube 1312 int MeshAdapter::getMPISize() const
77     {
78     return m_finleyMesh.get()->MPIInfo->size;
79     }
80     int MeshAdapter::getMPIRank() const
81     {
82     return m_finleyMesh.get()->MPIInfo->rank;
83     }
84 ksteube 1877 void MeshAdapter::MPIBarrier() const
85     {
86 jfenwick 3259 #ifdef ESYS_MPI
87 ksteube 1877 MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
88     #endif
89     return;
90     }
91     bool MeshAdapter::onMasterProcessor() const
92     {
93     return m_finleyMesh.get()->MPIInfo->rank == 0;
94     }
95 ksteube 1312
96    
97 jfenwick 3259 #ifdef ESYS_MPI
98 jfenwick 2640 MPI_Comm
99     #else
100     unsigned int
101     #endif
102     MeshAdapter::getMPIComm() const
103     {
104 jfenwick 3259 #ifdef ESYS_MPI
105 jfenwick 2641 return m_finleyMesh->MPIInfo->comm;
106 jfenwick 2640 #else
107     return 0;
108     #endif
109     }
110    
111    
112 jgs 82 Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
113     return m_finleyMesh.get();
114     }
115    
116 caltinay 2778 void MeshAdapter::write(const string& fileName) const
117 jgs 82 {
118 phornby 1628 char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
119     strcpy(fName,fileName.c_str());
120     Finley_Mesh_write(m_finleyMesh.get(),fName);
121     checkFinleyError();
122     TMPMEMFREE(fName);
123 jgs 82 }
124    
125 jfenwick 2519 void MeshAdapter::Print_Mesh_Info(const bool full) const
126 ksteube 1326 {
127 phornby 1628 Finley_PrintMesh_Info(m_finleyMesh.get(), full);
128 ksteube 1326 }
129    
130 caltinay 2778 void MeshAdapter::dump(const string& fileName) const
131 ksteube 1312 {
132 ksteube 1343 #ifdef USE_NETCDF
133 ksteube 1990 const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
134 trankine 1632 NcVar *ids;
135 ksteube 1343 int *int_ptr;
136     Finley_Mesh *mesh = m_finleyMesh.get();
137 ksteube 1347 Finley_TagMap* tag_map;
138     int num_Tags = 0;
139 ksteube 1343 int mpi_size = mesh->MPIInfo->size;
140     int mpi_rank = mesh->MPIInfo->rank;
141     int numDim = mesh->Nodes->numDim;
142     int numNodes = mesh->Nodes->numNodes;
143     int num_Elements = mesh->Elements->numElements;
144     int num_FaceElements = mesh->FaceElements->numElements;
145     int num_ContactElements = mesh->ContactElements->numElements;
146     int num_Points = mesh->Points->numElements;
147     int num_Elements_numNodes = mesh->Elements->numNodes;
148     int num_FaceElements_numNodes = mesh->FaceElements->numNodes;
149     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
150 jfenwick 3259 #ifdef ESYS_MPI
151 ksteube 1800 MPI_Status status;
152     #endif
153 ksteube 1343
154 ksteube 1800 /* Incoming token indicates it's my turn to write */
155 jfenwick 3259 #ifdef ESYS_MPI
156 ksteube 1800 if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
157     #endif
158    
159 jfenwick 3259 char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
160 phornby 1628 mpi_size, mpi_rank);
161    
162 ksteube 1347 /* Figure out how much storage is required for tags */
163     tag_map = mesh->TagMap;
164 ksteube 1800 num_Tags = 0;
165 caltinay 2778 while (tag_map) {
166     num_Tags++;
167     tag_map=tag_map->next;
168 ksteube 1347 }
169    
170 ksteube 1343 // NetCDF error handler
171     NcError err(NcError::verbose_nonfatal);
172     // Create the file.
173     NcFile dataFile(newFileName, NcFile::Replace);
174 caltinay 2778 string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
175 ksteube 1343 // check if writing was successful
176     if (!dataFile.is_valid())
177 caltinay 2778 throw DataException(msgPrefix+"Open file for output");
178 ksteube 1343
179 caltinay 2778 // Define dimensions (num_Elements and dim_Elements are identical,
180     // dim_Elements only appears if > 0)
181 ksteube 1343 if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
182 caltinay 2778 throw DataException(msgPrefix+"add_dim(numNodes)");
183 ksteube 1343 if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
184 caltinay 2778 throw DataException(msgPrefix+"add_dim(numDim)");
185 ksteube 1343 if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )
186 caltinay 2778 throw DataException(msgPrefix+"add_dim(mpi_size)");
187 ksteube 1343 if (num_Elements>0)
188     if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
189 caltinay 2778 throw DataException(msgPrefix+"add_dim(dim_Elements)");
190 ksteube 1343 if (num_FaceElements>0)
191     if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
192 caltinay 2778 throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
193 ksteube 1343 if (num_ContactElements>0)
194     if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
195 caltinay 2778 throw DataException(msgPrefix+"add_dim(dim_ContactElements)");
196 ksteube 1343 if (num_Points>0)
197     if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
198 caltinay 2778 throw DataException(msgPrefix+"add_dim(dim_Points)");
199 ksteube 1343 if (num_Elements>0)
200     if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
201 caltinay 2778 throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
202 ksteube 1343 if (num_FaceElements>0)
203     if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
204 caltinay 2778 throw DataException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
205 ksteube 1343 if (num_ContactElements>0)
206     if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
207 caltinay 2778 throw DataException(msgPrefix+"add_dim(dim_ContactElements_numNodes)");
208 ksteube 1347 if (num_Tags>0)
209     if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
210 caltinay 2778 throw DataException(msgPrefix+"add_dim(dim_Tags)");
211 ksteube 1343
212 ksteube 1345 // Attributes: MPI size, MPI rank, Name, order, reduced_order
213 ksteube 1343 if (!dataFile.add_att("mpi_size", mpi_size) )
214 caltinay 2778 throw DataException(msgPrefix+"add_att(mpi_size)");
215 ksteube 1343 if (!dataFile.add_att("mpi_rank", mpi_rank) )
216 caltinay 2778 throw DataException(msgPrefix+"add_att(mpi_rank)");
217 ksteube 1343 if (!dataFile.add_att("Name",mesh->Name) )
218 caltinay 2778 throw DataException(msgPrefix+"add_att(Name)");
219 ksteube 1345 if (!dataFile.add_att("numDim",numDim) )
220 caltinay 2778 throw DataException(msgPrefix+"add_att(order)");
221 gross 2856 if (!dataFile.add_att("order",mesh->integrationOrder) )
222 caltinay 2778 throw DataException(msgPrefix+"add_att(order)");
223 gross 2856 if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
224 caltinay 2778 throw DataException(msgPrefix+"add_att(reduced_order)");
225 ksteube 1345 if (!dataFile.add_att("numNodes",numNodes) )
226 caltinay 2778 throw DataException(msgPrefix+"add_att(numNodes)");
227 ksteube 1343 if (!dataFile.add_att("num_Elements",num_Elements) )
228 caltinay 2778 throw DataException(msgPrefix+"add_att(num_Elements)");
229 ksteube 1343 if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
230 caltinay 2778 throw DataException(msgPrefix+"add_att(num_FaceElements)");
231 ksteube 1343 if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
232 caltinay 2778 throw DataException(msgPrefix+"add_att(num_ContactElements)");
233 ksteube 1343 if (!dataFile.add_att("num_Points",num_Points) )
234 caltinay 2778 throw DataException(msgPrefix+"add_att(num_Points)");
235 ksteube 1346 if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
236 caltinay 2778 throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
237 ksteube 1346 if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
238 caltinay 2778 throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
239 ksteube 1346 if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
240 caltinay 2778 throw DataException(msgPrefix+"add_att(num_ContactElements_numNodes)");
241 gross 2748 if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
242 caltinay 2778 throw DataException(msgPrefix+"add_att(Elements_TypeId)");
243 gross 2748 if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
244 caltinay 2778 throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
245 gross 2748 if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
246 caltinay 2778 throw DataException(msgPrefix+"add_att(ContactElements_TypeId)");
247 gross 2748 if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
248 caltinay 2778 throw DataException(msgPrefix+"add_att(Points_TypeId)");
249 ksteube 1347 if (!dataFile.add_att("num_Tags", num_Tags) )
250 caltinay 2778 throw DataException(msgPrefix+"add_att(num_Tags)");
251 ksteube 1343
252     // // // // // Nodes // // // // //
253    
254 caltinay 2842 // Nodes nodeDistribution
255     if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
256     throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
257     int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
258     if (! (ids->put(int_ptr, mpi_size+1)) )
259     throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
260    
261     // Nodes degreesOfFreedomDistribution
262     if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
263     throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
264     int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
265     if (! (ids->put(int_ptr, mpi_size+1)) )
266     throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
267    
268 caltinay 2778 // Only write nodes if non-empty because NetCDF doesn't like empty arrays
269     // (it treats them as NC_UNLIMITED)
270 ksteube 1345 if (numNodes>0) {
271    
272 phornby 1628 // Nodes Id
273     if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
274 caltinay 2778 throw DataException(msgPrefix+"add_var(Nodes_Id)");
275 phornby 1628 int_ptr = &mesh->Nodes->Id[0];
276     if (! (ids->put(int_ptr, numNodes)) )
277 caltinay 2778 throw DataException(msgPrefix+"put(Nodes_Id)");
278 ksteube 1343
279 phornby 1628 // Nodes Tag
280     if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
281 caltinay 2778 throw DataException(msgPrefix+"add_var(Nodes_Tag)");
282 phornby 1628 int_ptr = &mesh->Nodes->Tag[0];
283     if (! (ids->put(int_ptr, numNodes)) )
284 caltinay 2778 throw DataException(msgPrefix+"put(Nodes_Tag)");
285 ksteube 1343
286 phornby 1628 // Nodes gDOF
287     if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
288 caltinay 2778 throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
289 phornby 1628 int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
290     if (! (ids->put(int_ptr, numNodes)) )
291 caltinay 2778 throw DataException(msgPrefix+"put(Nodes_gDOF)");
292 ksteube 1343
293 phornby 1628 // Nodes global node index
294     if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
295 caltinay 2778 throw DataException(msgPrefix+"add_var(Nodes_gNI)");
296 phornby 1628 int_ptr = &mesh->Nodes->globalNodesIndex[0];
297     if (! (ids->put(int_ptr, numNodes)) )
298 caltinay 2778 throw DataException(msgPrefix+"put(Nodes_gNI)");
299 ksteube 1343
300 phornby 1628 // Nodes grDof
301     if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
302 caltinay 2778 throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
303 phornby 1628 int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
304     if (! (ids->put(int_ptr, numNodes)) )
305 caltinay 2778 throw DataException(msgPrefix+"put(Nodes_grDfI)");
306 ksteube 1343
307 phornby 1628 // Nodes grNI
308     if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
309 caltinay 2778 throw DataException(msgPrefix+"add_var(Nodes_grNI)");
310 phornby 1628 int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
311     if (! (ids->put(int_ptr, numNodes)) )
312 caltinay 2778 throw DataException(msgPrefix+"put(Nodes_grNI)");
313 ksteube 1343
314 phornby 1628 // Nodes Coordinates
315     if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
316 caltinay 2778 throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
317 phornby 1628 if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
318 caltinay 2778 throw DataException(msgPrefix+"put(Nodes_Coordinates)");
319 ksteube 1343
320 ksteube 1345 }
321    
322 ksteube 1343 // // // // // Elements // // // // //
323    
324     if (num_Elements>0) {
325    
326 phornby 1628 // Elements_Id
327     if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
328 caltinay 2778 throw DataException(msgPrefix+"add_var(Elements_Id)");
329 phornby 1628 int_ptr = &mesh->Elements->Id[0];
330     if (! (ids->put(int_ptr, num_Elements)) )
331 caltinay 2778 throw DataException(msgPrefix+"put(Elements_Id)");
332 ksteube 1343
333 phornby 1628 // Elements_Tag
334     if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
335 caltinay 2778 throw DataException(msgPrefix+"add_var(Elements_Tag)");
336 phornby 1628 int_ptr = &mesh->Elements->Tag[0];
337     if (! (ids->put(int_ptr, num_Elements)) )
338 caltinay 2778 throw DataException(msgPrefix+"put(Elements_Tag)");
339 ksteube 1343
340 phornby 1628 // Elements_Owner
341     if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
342 caltinay 2778 throw DataException(msgPrefix+"add_var(Elements_Owner)");
343 phornby 1628 int_ptr = &mesh->Elements->Owner[0];
344     if (! (ids->put(int_ptr, num_Elements)) )
345 caltinay 2778 throw DataException(msgPrefix+"put(Elements_Owner)");
346 ksteube 1343
347 phornby 1628 // Elements_Color
348     if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
349 caltinay 2778 throw DataException(msgPrefix+"add_var(Elements_Color)");
350 phornby 1628 int_ptr = &mesh->Elements->Color[0];
351     if (! (ids->put(int_ptr, num_Elements)) )
352 caltinay 2778 throw DataException(msgPrefix+"put(Elements_Color)");
353 ksteube 1343
354 phornby 1628 // Elements_Nodes
355     if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
356 caltinay 2778 throw DataException(msgPrefix+"add_var(Elements_Nodes)");
357 ksteube 1777 if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
358 caltinay 2778 throw DataException(msgPrefix+"put(Elements_Nodes)");
359 ksteube 1343
360     }
361    
362     // // // // // Face_Elements // // // // //
363    
364     if (num_FaceElements>0) {
365    
366 phornby 1628 // FaceElements_Id
367     if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
368 caltinay 2778 throw DataException(msgPrefix+"add_var(FaceElements_Id)");
369 phornby 1628 int_ptr = &mesh->FaceElements->Id[0];
370     if (! (ids->put(int_ptr, num_FaceElements)) )
371 caltinay 2778 throw DataException(msgPrefix+"put(FaceElements_Id)");
372 ksteube 1343
373 phornby 1628 // FaceElements_Tag
374     if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
375 caltinay 2778 throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
376 phornby 1628 int_ptr = &mesh->FaceElements->Tag[0];
377     if (! (ids->put(int_ptr, num_FaceElements)) )
378 caltinay 2778 throw DataException(msgPrefix+"put(FaceElements_Tag)");
379 ksteube 1343
380 phornby 1628 // FaceElements_Owner
381     if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
382 caltinay 2778 throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
383 phornby 1628 int_ptr = &mesh->FaceElements->Owner[0];
384     if (! (ids->put(int_ptr, num_FaceElements)) )
385 caltinay 2778 throw DataException(msgPrefix+"put(FaceElements_Owner)");
386 ksteube 1343
387 phornby 1628 // FaceElements_Color
388     if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
389 caltinay 2778 throw DataException(msgPrefix+"add_var(FaceElements_Color)");
390 phornby 1628 int_ptr = &mesh->FaceElements->Color[0];
391     if (! (ids->put(int_ptr, num_FaceElements)) )
392 caltinay 2778 throw DataException(msgPrefix+"put(FaceElements_Color)");
393 ksteube 1343
394 phornby 1628 // FaceElements_Nodes
395     if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
396 caltinay 2778 throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
397 ksteube 1777 if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
398 caltinay 2778 throw DataException(msgPrefix+"put(FaceElements_Nodes)");
399 ksteube 1343
400     }
401    
402     // // // // // Contact_Elements // // // // //
403    
404     if (num_ContactElements>0) {
405    
406 phornby 1628 // ContactElements_Id
407     if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
408 caltinay 2778 throw DataException(msgPrefix+"add_var(ContactElements_Id)");
409 phornby 1628 int_ptr = &mesh->ContactElements->Id[0];
410     if (! (ids->put(int_ptr, num_ContactElements)) )
411 caltinay 2778 throw DataException(msgPrefix+"put(ContactElements_Id)");
412 ksteube 1343
413 phornby 1628 // ContactElements_Tag
414     if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
415 caltinay 2778 throw DataException(msgPrefix+"add_var(ContactElements_Tag)");
416 phornby 1628 int_ptr = &mesh->ContactElements->Tag[0];
417     if (! (ids->put(int_ptr, num_ContactElements)) )
418 caltinay 2778 throw DataException(msgPrefix+"put(ContactElements_Tag)");
419 ksteube 1343
420 phornby 1628 // ContactElements_Owner
421     if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
422 caltinay 2778 throw DataException(msgPrefix+"add_var(ContactElements_Owner)");
423 phornby 1628 int_ptr = &mesh->ContactElements->Owner[0];
424     if (! (ids->put(int_ptr, num_ContactElements)) )
425 caltinay 2778 throw DataException(msgPrefix+"put(ContactElements_Owner)");
426 ksteube 1343
427 phornby 1628 // ContactElements_Color
428     if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
429 caltinay 2778 throw DataException(msgPrefix+"add_var(ContactElements_Color)");
430 phornby 1628 int_ptr = &mesh->ContactElements->Color[0];
431     if (! (ids->put(int_ptr, num_ContactElements)) )
432 caltinay 2778 throw DataException(msgPrefix+"put(ContactElements_Color)");
433 ksteube 1343
434 phornby 1628 // ContactElements_Nodes
435     if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
436 caltinay 2778 throw DataException(msgPrefix+"add_var(ContactElements_Nodes)");
437 ksteube 1777 if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
438 caltinay 2778 throw DataException(msgPrefix+"put(ContactElements_Nodes)");
439 ksteube 1343
440     }
441    
442     // // // // // Points // // // // //
443    
444     if (num_Points>0) {
445    
446 phornby 1628 fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");
447 ksteube 1343
448 phornby 1628 // Points_Id
449     if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
450 caltinay 2778 throw DataException(msgPrefix+"add_var(Points_Id)");
451 phornby 1628 int_ptr = &mesh->Points->Id[0];
452     if (! (ids->put(int_ptr, num_Points)) )
453 caltinay 2778 throw DataException(msgPrefix+"put(Points_Id)");
454 ksteube 1343
455 phornby 1628 // Points_Tag
456     if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
457 caltinay 2778 throw DataException(msgPrefix+"add_var(Points_Tag)");
458 phornby 1628 int_ptr = &mesh->Points->Tag[0];
459     if (! (ids->put(int_ptr, num_Points)) )
460 caltinay 2778 throw DataException(msgPrefix+"put(Points_Tag)");
461 ksteube 1343
462 phornby 1628 // Points_Owner
463     if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
464 caltinay 2778 throw DataException(msgPrefix+"add_var(Points_Owner)");
465 phornby 1628 int_ptr = &mesh->Points->Owner[0];
466     if (! (ids->put(int_ptr, num_Points)) )
467 caltinay 2778 throw DataException(msgPrefix+"put(Points_Owner)");
468 ksteube 1343
469 phornby 1628 // Points_Color
470     if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
471 caltinay 2778 throw DataException(msgPrefix+"add_var(Points_Color)");
472 phornby 1628 int_ptr = &mesh->Points->Color[0];
473     if (! (ids->put(int_ptr, num_Points)) )
474 caltinay 2778 throw DataException(msgPrefix+"put(Points_Color)");
475 ksteube 1343
476 phornby 1628 // Points_Nodes
477     // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
478     if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
479 caltinay 2778 throw DataException(msgPrefix+"add_var(Points_Nodes)");
480 ksteube 1777 if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
481 caltinay 2778 throw DataException(msgPrefix+"put(Points_Nodes)");
482 ksteube 1343
483     }
484    
485     // // // // // TagMap // // // // //
486    
487 ksteube 1347 if (num_Tags>0) {
488 ksteube 1343
489 phornby 1628 // Temp storage to gather node IDs
490     int *Tags_keys = TMPMEMALLOC(num_Tags, int);
491     char name_temp[4096];
492 ksteube 1347
493 phornby 1628 /* Copy tag data into temp arrays */
494     tag_map = mesh->TagMap;
495     if (tag_map) {
496     int i = 0;
497     while (tag_map) {
498     Tags_keys[i++] = tag_map->tag_key;
499     tag_map=tag_map->next;
500     }
501     }
502 ksteube 1347
503 phornby 1628 // Tags_keys
504     if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
505 caltinay 2778 throw DataException(msgPrefix+"add_var(Tags_keys)");
506 phornby 1628 int_ptr = &Tags_keys[0];
507     if (! (ids->put(int_ptr, num_Tags)) )
508 caltinay 2778 throw DataException(msgPrefix+"put(Tags_keys)");
509 ksteube 1347
510 phornby 1628 // Tags_names_*
511 caltinay 2778 // This is an array of strings, it should be stored as an array but
512     // instead I have hacked in one attribute per string because the NetCDF
513     // manual doesn't tell how to do an array of strings
514 phornby 1628 tag_map = mesh->TagMap;
515     if (tag_map) {
516     int i = 0;
517     while (tag_map) {
518     sprintf(name_temp, "Tags_name_%d", i);
519     if (!dataFile.add_att(name_temp, tag_map->name) )
520 caltinay 2778 throw DataException(msgPrefix+"add_att(Tags_names_XX)");
521 phornby 1628 tag_map=tag_map->next;
522     i++;
523     }
524     }
525 ksteube 1347
526 phornby 1628 TMPMEMFREE(Tags_keys);
527 ksteube 1347 }
528    
529 ksteube 1800 /* Send token to next MPI process so he can take his turn */
530 jfenwick 3259 #ifdef ESYS_MPI
531 ksteube 1800 if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
532     #endif
533 ksteube 1347
534 ksteube 1343 // NetCDF file is closed by destructor of NcFile object
535 ksteube 1800
536 ksteube 1343 #else
537     Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
538     #endif /* USE_NETCDF */
539     checkFinleyError();
540 ksteube 1312 }
541    
542 jgs 82 string MeshAdapter::getDescription() const
543     {
544 phornby 1628 return "FinleyMesh";
545 jgs 82 }
546    
547     string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
548     {
549 phornby 1628 FunctionSpaceNamesMapType::iterator loc;
550     loc=m_functionSpaceTypeNames.find(functionSpaceType);
551     if (loc==m_functionSpaceTypeNames.end()) {
552     return "Invalid function space type code.";
553     } else {
554     return loc->second;
555     }
556 jgs 82 }
557    
558     bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
559     {
560 phornby 1628 FunctionSpaceNamesMapType::iterator loc;
561     loc=m_functionSpaceTypeNames.find(functionSpaceType);
562     return (loc!=m_functionSpaceTypeNames.end());
563 jgs 82 }
564    
565     void MeshAdapter::setFunctionSpaceTypeNames()
566     {
567 phornby 1628 m_functionSpaceTypeNames.insert
568     (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));
569     m_functionSpaceTypeNames.insert
570     (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));
571     m_functionSpaceTypeNames.insert
572     (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
573     m_functionSpaceTypeNames.insert
574     (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));
575     m_functionSpaceTypeNames.insert
576     (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
577     m_functionSpaceTypeNames.insert
578     (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
579     m_functionSpaceTypeNames.insert
580     (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
581     m_functionSpaceTypeNames.insert
582     (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
583     m_functionSpaceTypeNames.insert
584     (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
585     m_functionSpaceTypeNames.insert
586     (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
587     m_functionSpaceTypeNames.insert
588     (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
589     m_functionSpaceTypeNames.insert
590     (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
591     m_functionSpaceTypeNames.insert
592     (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
593 jgs 82 }
594    
595     int MeshAdapter::getContinuousFunctionCode() const
596     {
597 phornby 1628 return Nodes;
598 jgs 82 }
599 gross 1062 int MeshAdapter::getReducedContinuousFunctionCode() const
600     {
601 phornby 1628 return ReducedNodes;
602 gross 1062 }
603 jgs 149
604 jgs 82 int MeshAdapter::getFunctionCode() const
605     {
606 phornby 1628 return Elements;
607 jgs 82 }
608 gross 1059 int MeshAdapter::getReducedFunctionCode() const
609     {
610 phornby 1628 return ReducedElements;
611 gross 1059 }
612 jgs 149
613 jgs 82 int MeshAdapter::getFunctionOnBoundaryCode() const
614     {
615 phornby 1628 return FaceElements;
616 jgs 82 }
617 gross 1059 int MeshAdapter::getReducedFunctionOnBoundaryCode() const
618     {
619 phornby 1628 return ReducedFaceElements;
620 gross 1059 }
621 jgs 149
622 jgs 82 int MeshAdapter::getFunctionOnContactZeroCode() const
623     {
624 phornby 1628 return ContactElementsZero;
625 jgs 82 }
626 gross 1059 int MeshAdapter::getReducedFunctionOnContactZeroCode() const
627     {
628 phornby 1628 return ReducedContactElementsZero;
629 gross 1059 }
630 jgs 149
631 jgs 82 int MeshAdapter::getFunctionOnContactOneCode() const
632     {
633 phornby 1628 return ContactElementsOne;
634 jgs 82 }
635 gross 1059 int MeshAdapter::getReducedFunctionOnContactOneCode() const
636     {
637 phornby 1628 return ReducedContactElementsOne;
638 gross 1059 }
639 jgs 82
640     int MeshAdapter::getSolutionCode() const
641     {
642 phornby 1628 return DegreesOfFreedom;
643 jgs 82 }
644 jgs 149
645 jgs 82 int MeshAdapter::getReducedSolutionCode() const
646     {
647 phornby 1628 return ReducedDegreesOfFreedom;
648 jgs 82 }
649 jgs 149
650 jgs 82 int MeshAdapter::getDiracDeltaFunctionCode() const
651     {
652 phornby 1628 return Points;
653 jgs 82 }
654 jgs 149
655 jgs 82 //
656     // return the spatial dimension of the Mesh:
657     //
658     int MeshAdapter::getDim() const
659     {
660 gross 2100 Finley_Mesh* mesh=m_finleyMesh.get();
661     int numDim=Finley_Mesh_getDim(mesh);
662 phornby 1628 checkFinleyError();
663     return numDim;
664 jgs 82 }
665 jgs 149
666 jgs 82 //
667 ksteube 1754 // Return the number of data points summed across all MPI processes
668     //
669     int MeshAdapter::getNumDataPointsGlobal() const
670     {
671     return Finley_NodeFile_getGlobalNumNodes(m_finleyMesh.get()->Nodes);
672     }
673    
674     //
675 jgs 82 // return the number of data points per sample and the number of samples
676     // needed to represent data on a parts of the mesh.
677     //
678     pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
679     {
680     int numDataPointsPerSample=0;
681     int numSamples=0;
682     Finley_Mesh* mesh=m_finleyMesh.get();
683     switch (functionSpaceCode) {
684 phornby 1628 case(Nodes):
685     numDataPointsPerSample=1;
686     numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);
687     break;
688     case(ReducedNodes):
689     numDataPointsPerSample=1;
690     numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
691     break;
692     case(Elements):
693     if (mesh->Elements!=NULL) {
694     numSamples=mesh->Elements->numElements;
695 gross 2748 numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
696 phornby 1628 }
697     break;
698     case(ReducedElements):
699     if (mesh->Elements!=NULL) {
700     numSamples=mesh->Elements->numElements;
701 gross 2748 numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
702 phornby 1628 }
703     break;
704     case(FaceElements):
705     if (mesh->FaceElements!=NULL) {
706 gross 2748 numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
707 phornby 1628 numSamples=mesh->FaceElements->numElements;
708     }
709     break;
710     case(ReducedFaceElements):
711     if (mesh->FaceElements!=NULL) {
712 gross 2748 numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
713 phornby 1628 numSamples=mesh->FaceElements->numElements;
714     }
715     break;
716     case(Points):
717     if (mesh->Points!=NULL) {
718     numDataPointsPerSample=1;
719     numSamples=mesh->Points->numElements;
720     }
721     break;
722     case(ContactElementsZero):
723     if (mesh->ContactElements!=NULL) {
724 gross 2748 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
725 phornby 1628 numSamples=mesh->ContactElements->numElements;
726     }
727     break;
728     case(ReducedContactElementsZero):
729     if (mesh->ContactElements!=NULL) {
730 gross 2748 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
731 phornby 1628 numSamples=mesh->ContactElements->numElements;
732     }
733     break;
734     case(ContactElementsOne):
735     if (mesh->ContactElements!=NULL) {
736 gross 2748 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
737 phornby 1628 numSamples=mesh->ContactElements->numElements;
738     }
739     break;
740     case(ReducedContactElementsOne):
741     if (mesh->ContactElements!=NULL) {
742 gross 2748 numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
743 phornby 1628 numSamples=mesh->ContactElements->numElements;
744     }
745     break;
746     case(DegreesOfFreedom):
747     if (mesh->Nodes!=NULL) {
748     numDataPointsPerSample=1;
749     numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
750     }
751     break;
752     case(ReducedDegreesOfFreedom):
753     if (mesh->Nodes!=NULL) {
754     numDataPointsPerSample=1;
755     numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
756     }
757     break;
758     default:
759     stringstream temp;
760     temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
761     throw FinleyAdapterException(temp.str());
762     break;
763     }
764     return pair<int,int>(numDataPointsPerSample,numSamples);
765 jgs 82 }
766 jgs 149
767 jgs 82 //
768     // adds linear PDE of second order into a given stiffness matrix and right hand side:
769     //
770     void MeshAdapter::addPDEToSystem(
771 jfenwick 3269 AbstractSystemMatrix& mat, escript::Data& rhs,
772 phornby 1628 const escript::Data& A, const escript::Data& B, const escript::Data& C,const escript::Data& D,const escript::Data& X,const escript::Data& Y,
773     const escript::Data& d, const escript::Data& y,
774     const escript::Data& d_contact,const escript::Data& y_contact) const
775 jgs 82 {
776 jfenwick 3269 SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
777     if (smat==0)
778     {
779     throw FinleyAdapterException("finley only supports its own system matrices.");
780     }
781 gross 798 escriptDataC _rhs=rhs.getDataC();
782     escriptDataC _A =A.getDataC();
783     escriptDataC _B=B.getDataC();
784     escriptDataC _C=C.getDataC();
785     escriptDataC _D=D.getDataC();
786     escriptDataC _X=X.getDataC();
787     escriptDataC _Y=Y.getDataC();
788     escriptDataC _d=d.getDataC();
789     escriptDataC _y=y.getDataC();
790     escriptDataC _d_contact=d_contact.getDataC();
791     escriptDataC _y_contact=y_contact.getDataC();
792 bcumming 751
793 jgs 82 Finley_Mesh* mesh=m_finleyMesh.get();
794 gross 798
795 jfenwick 3269 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
796 jgs 82 checkFinleyError();
797 bcumming 751
798 jfenwick 3269 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
799 jgs 82 checkFinleyError();
800 bcumming 751
801 jfenwick 3269 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
802 jgs 82 checkFinleyError();
803     }
804 jgs 149
805 gross 1204 void MeshAdapter::addPDEToLumpedSystem(
806 phornby 1628 escript::Data& mat,
807     const escript::Data& D,
808     const escript::Data& d) const
809 gross 1204 {
810     escriptDataC _mat=mat.getDataC();
811     escriptDataC _D=D.getDataC();
812     escriptDataC _d=d.getDataC();
813    
814     Finley_Mesh* mesh=m_finleyMesh.get();
815 ksteube 1312
816 gross 1204 Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
817 gross 2989 checkFinleyError();
818    
819 gross 1204 Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
820 gross 2989 checkFinleyError();
821 gross 1204
822     }
823    
824    
825 jgs 82 //
826 jgs 102 // adds linear PDE of second order into the right hand side only
827     //
828 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
829 jgs 102 {
830     Finley_Mesh* mesh=m_finleyMesh.get();
831 jgs 147
832 gross 798 escriptDataC _rhs=rhs.getDataC();
833     escriptDataC _X=X.getDataC();
834     escriptDataC _Y=Y.getDataC();
835     escriptDataC _y=y.getDataC();
836     escriptDataC _y_contact=y_contact.getDataC();
837    
838     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
839 jgs 102 checkFinleyError();
840 jgs 147
841 gross 798 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
842     checkFinleyError();
843 jgs 147
844 gross 798 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
845 jgs 102 checkFinleyError();
846     }
847 gross 1367 //
848     // adds PDE of second order into a transport problem
849     //
850     void MeshAdapter::addPDEToTransportProblem(
851 jfenwick 3269 AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
852 phornby 1628 const escript::Data& A, const escript::Data& B, const escript::Data& C,
853     const escript::Data& D,const escript::Data& X,const escript::Data& Y,
854     const escript::Data& d, const escript::Data& y,
855     const escript::Data& d_contact,const escript::Data& y_contact) const
856 gross 1367 {
857 jfenwick 3269 TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
858     if (tpa==0)
859     {
860     throw FinleyAdapterException("finley only supports its own transport problems.");
861     }
862    
863    
864 jfenwick 1796 DataTypes::ShapeType shape;
865 gross 1370 source.expand();
866 gross 1367 escriptDataC _source=source.getDataC();
867     escriptDataC _M=M.getDataC();
868     escriptDataC _A=A.getDataC();
869     escriptDataC _B=B.getDataC();
870     escriptDataC _C=C.getDataC();
871     escriptDataC _D=D.getDataC();
872     escriptDataC _X=X.getDataC();
873     escriptDataC _Y=Y.getDataC();
874     escriptDataC _d=d.getDataC();
875     escriptDataC _y=y.getDataC();
876     escriptDataC _d_contact=d_contact.getDataC();
877     escriptDataC _y_contact=y_contact.getDataC();
878 jgs 149
879 gross 1367 Finley_Mesh* mesh=m_finleyMesh.get();
880 jfenwick 3269 Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
881 phornby 1628
882 gross 1407 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
883 gross 1367 checkFinleyError();
884    
885 gross 1407 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
886 gross 1367 checkFinleyError();
887    
888     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
889     checkFinleyError();
890    
891     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
892     checkFinleyError();
893     }
894    
895 jgs 102 //
896 jgs 82 // interpolates data between different function spaces:
897     //
898 jgs 480 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
899 jgs 82 {
900 jfenwick 1872 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
901     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
902 phornby 1628 if (inDomain!=*this)
903     throw FinleyAdapterException("Error - Illegal domain of interpolant.");
904     if (targetDomain!=*this)
905     throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
906 jgs 82
907 phornby 1628 Finley_Mesh* mesh=m_finleyMesh.get();
908     escriptDataC _target=target.getDataC();
909     escriptDataC _in=in.getDataC();
910     switch(in.getFunctionSpace().getTypeCode()) {
911     case(Nodes):
912 gross 2385 switch(target.getFunctionSpace().getTypeCode()) {
913     case(Nodes):
914     case(ReducedNodes):
915     case(DegreesOfFreedom):
916     case(ReducedDegreesOfFreedom):
917     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
918 phornby 1628 break;
919 gross 2385 case(Elements):
920     case(ReducedElements):
921     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
922 phornby 1628 break;
923 gross 2385 case(FaceElements):
924     case(ReducedFaceElements):
925     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
926     break;
927     case(Points):
928     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
929     break;
930     case(ContactElementsZero):
931     case(ReducedContactElementsZero):
932     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
933     break;
934     case(ContactElementsOne):
935     case(ReducedContactElementsOne):
936     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
937     break;
938     default:
939     stringstream temp;
940     temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
941     throw FinleyAdapterException(temp.str());
942     break;
943     }
944     break;
945 phornby 1628 case(ReducedNodes):
946 gross 2385 switch(target.getFunctionSpace().getTypeCode()) {
947     case(Nodes):
948     case(ReducedNodes):
949     case(DegreesOfFreedom):
950     case(ReducedDegreesOfFreedom):
951 phornby 1628 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
952 gross 2385 break;
953     case(Elements):
954     case(ReducedElements):
955 phornby 1628 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
956 gross 2385 break;
957     case(FaceElements):
958     case(ReducedFaceElements):
959 phornby 1628 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
960 gross 2385 break;
961     case(Points):
962 phornby 1628 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
963 gross 2385 break;
964     case(ContactElementsZero):
965     case(ReducedContactElementsZero):
966 phornby 1628 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
967     break;
968 gross 2385 case(ContactElementsOne):
969     case(ReducedContactElementsOne):
970     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
971     break;
972     default:
973     stringstream temp;
974     temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
975     throw FinleyAdapterException(temp.str());
976     break;
977     }
978     break;
979 phornby 1628 case(Elements):
980 gross 2385 if (target.getFunctionSpace().getTypeCode()==Elements) {
981     Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
982     } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
983     Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
984     } else {
985     throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
986     }
987     break;
988 phornby 1628 case(ReducedElements):
989 gross 2385 if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
990     Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
991     } else {
992     throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
993     }
994     break;
995 phornby 1628 case(FaceElements):
996 gross 2385 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
997     Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
998     } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
999     Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1000     } else {
1001     throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1002     }
1003     break;
1004 phornby 1628 case(ReducedFaceElements):
1005 gross 2385 if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1006     Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1007     } else {
1008     throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
1009     }
1010     break;
1011 phornby 1628 case(Points):
1012 gross 2385 if (target.getFunctionSpace().getTypeCode()==Points) {
1013     Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1014     } else {
1015     throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1016     }
1017     break;
1018 phornby 1628 case(ContactElementsZero):
1019     case(ContactElementsOne):
1020 gross 2385 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1021     Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1022     } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1023     Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1024     } else {
1025     throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1026     }
1027     break;
1028 phornby 1628 case(ReducedContactElementsZero):
1029     case(ReducedContactElementsOne):
1030 gross 2385 if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1031     Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1032     } else {
1033     throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1034     }
1035 phornby 1628 break;
1036 gross 2385 case(DegreesOfFreedom):
1037     switch(target.getFunctionSpace().getTypeCode()) {
1038     case(ReducedDegreesOfFreedom):
1039     case(DegreesOfFreedom):
1040     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1041     break;
1042    
1043     case(Nodes):
1044     case(ReducedNodes):
1045     if (getMPISize()>1) {
1046     escript::Data temp=escript::Data(in);
1047     temp.expand();
1048     escriptDataC _in2 = temp.getDataC();
1049     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1050     } else {
1051     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1052     }
1053     break;
1054     case(Elements):
1055     case(ReducedElements):
1056     if (getMPISize()>1) {
1057 jfenwick 3259 escript::Data temp=escript::Data( in, continuousFunction(*this) );
1058 gross 2385 escriptDataC _in2 = temp.getDataC();
1059     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1060     } else {
1061     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1062     }
1063     break;
1064     case(FaceElements):
1065     case(ReducedFaceElements):
1066     if (getMPISize()>1) {
1067 jfenwick 3259 escript::Data temp=escript::Data( in, continuousFunction(*this) );
1068 gross 2385 escriptDataC _in2 = temp.getDataC();
1069     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1070    
1071     } else {
1072     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1073     }
1074     break;
1075     case(Points):
1076     if (getMPISize()>1) {
1077 jfenwick 3259 escript::Data temp=escript::Data( in, continuousFunction(*this) );
1078 gross 2385 escriptDataC _in2 = temp.getDataC();
1079     } else {
1080     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1081     }
1082     break;
1083     case(ContactElementsZero):
1084     case(ContactElementsOne):
1085     case(ReducedContactElementsZero):
1086     case(ReducedContactElementsOne):
1087     if (getMPISize()>1) {
1088 jfenwick 3259 escript::Data temp=escript::Data( in, continuousFunction(*this) );
1089 gross 2385 escriptDataC _in2 = temp.getDataC();
1090     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1091     } else {
1092     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1093     }
1094     break;
1095     default:
1096     stringstream temp;
1097     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1098     throw FinleyAdapterException(temp.str());
1099     break;
1100     }
1101     break;
1102     case(ReducedDegreesOfFreedom):
1103     switch(target.getFunctionSpace().getTypeCode()) {
1104     case(Nodes):
1105     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1106     break;
1107     case(ReducedNodes):
1108     if (getMPISize()>1) {
1109     escript::Data temp=escript::Data(in);
1110     temp.expand();
1111     escriptDataC _in2 = temp.getDataC();
1112     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1113     } else {
1114     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1115     }
1116     break;
1117     case(DegreesOfFreedom):
1118     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1119     break;
1120     case(ReducedDegreesOfFreedom):
1121     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1122     break;
1123     case(Elements):
1124     case(ReducedElements):
1125     if (getMPISize()>1) {
1126 jfenwick 3259 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
1127 gross 2385 escriptDataC _in2 = temp.getDataC();
1128     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1129     } else {
1130     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1131     }
1132     break;
1133     case(FaceElements):
1134     case(ReducedFaceElements):
1135     if (getMPISize()>1) {
1136 jfenwick 3259 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
1137 gross 2385 escriptDataC _in2 = temp.getDataC();
1138     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1139     } else {
1140     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1141     }
1142     break;
1143     case(Points):
1144     if (getMPISize()>1) {
1145 jfenwick 3259 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
1146 gross 2385 escriptDataC _in2 = temp.getDataC();
1147     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1148     } else {
1149     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1150     }
1151     break;
1152     case(ContactElementsZero):
1153     case(ContactElementsOne):
1154     case(ReducedContactElementsZero):
1155     case(ReducedContactElementsOne):
1156     if (getMPISize()>1) {
1157 jfenwick 3259 escript::Data temp=escript::Data( in, reducedContinuousFunction(*this) );
1158 gross 2385 escriptDataC _in2 = temp.getDataC();
1159     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1160     } else {
1161     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1162     }
1163     break;
1164     default:
1165     stringstream temp;
1166     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1167     throw FinleyAdapterException(temp.str());
1168     break;
1169     }
1170     break;
1171 phornby 1628 default:
1172     stringstream temp;
1173     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1174     throw FinleyAdapterException(temp.str());
1175     break;
1176     }
1177     checkFinleyError();
1178 jgs 82 }
1179    
1180     //
1181     // copies the locations of sample points into x:
1182     //
1183 jgs 480 void MeshAdapter::setToX(escript::Data& arg) const
1184 jgs 82 {
1185 jfenwick 1872 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1186 phornby 1628 if (argDomain!=*this)
1187     throw FinleyAdapterException("Error - Illegal domain of data point locations");
1188     Finley_Mesh* mesh=m_finleyMesh.get();
1189     // in case of values node coordinates we can do the job directly:
1190     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1191     escriptDataC _arg=arg.getDataC();
1192     Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1193     } else {
1194 jfenwick 3259 escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1195 phornby 1628 escriptDataC _tmp_data=tmp_data.getDataC();
1196     Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1197     // this is then interpolated onto arg:
1198     interpolateOnDomain(arg,tmp_data);
1199     }
1200     checkFinleyError();
1201 jgs 82 }
1202 jgs 149
1203 jgs 82 //
1204     // return the normal vectors at the location of data points as a Data object:
1205     //
1206 jgs 480 void MeshAdapter::setToNormal(escript::Data& normal) const
1207 jgs 82 {
1208 jfenwick 1872 /* const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1209     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1210 phornby 1628 if (normalDomain!=*this)
1211     throw FinleyAdapterException("Error - Illegal domain of normal locations");
1212     Finley_Mesh* mesh=m_finleyMesh.get();
1213     escriptDataC _normal=normal.getDataC();
1214     switch(normal.getFunctionSpace().getTypeCode()) {
1215     case(Nodes):
1216     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1217     break;
1218     case(ReducedNodes):
1219     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1220     break;
1221     case(Elements):
1222     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1223     break;
1224     case(ReducedElements):
1225     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1226     break;
1227     case (FaceElements):
1228     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1229     break;
1230     case (ReducedFaceElements):
1231     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1232     break;
1233     case(Points):
1234     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1235     break;
1236     case (ContactElementsOne):
1237     case (ContactElementsZero):
1238     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1239     break;
1240     case (ReducedContactElementsOne):
1241     case (ReducedContactElementsZero):
1242     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1243     break;
1244     case(DegreesOfFreedom):
1245     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1246     break;
1247     case(ReducedDegreesOfFreedom):
1248     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1249     break;
1250     default:
1251 jgs 150 stringstream temp;
1252     temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1253     throw FinleyAdapterException(temp.str());
1254 jgs 82 break;
1255 phornby 1628 }
1256     checkFinleyError();
1257 jgs 82 }
1258 jgs 149
1259 jgs 82 //
1260     // interpolates data to other domain:
1261     //
1262 jgs 480 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1263 jgs 82 {
1264 jfenwick 2087 const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1265     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1266     if (targetDomain!=this)
1267 phornby 1628 throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1268 jgs 82
1269 phornby 1628 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
1270 jgs 82 }
1271 jgs 149
1272 jgs 82 //
1273     // calculates the integral of a function defined of arg:
1274     //
1275 caltinay 2778 void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1276 jgs 82 {
1277 jfenwick 1872 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1278 phornby 1628 if (argDomain!=*this)
1279     throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1280 jgs 82
1281 phornby 1628 double blocktimer_start = blocktimer_time();
1282     Finley_Mesh* mesh=m_finleyMesh.get();
1283     escriptDataC _temp;
1284     escript::Data temp;
1285     escriptDataC _arg=arg.getDataC();
1286     switch(arg.getFunctionSpace().getTypeCode()) {
1287     case(Nodes):
1288 jfenwick 3259 temp=escript::Data( arg, escript::function(*this) );
1289 phornby 1628 _temp=temp.getDataC();
1290     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1291     break;
1292     case(ReducedNodes):
1293 jfenwick 3259 temp=escript::Data( arg, escript::function(*this) );
1294 phornby 1628 _temp=temp.getDataC();
1295     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1296     break;
1297     case(Elements):
1298     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1299     break;
1300     case(ReducedElements):
1301     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1302     break;
1303     case(FaceElements):
1304     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1305     break;
1306     case(ReducedFaceElements):
1307     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1308     break;
1309     case(Points):
1310     throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1311     break;
1312     case(ContactElementsZero):
1313     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1314     break;
1315     case(ReducedContactElementsZero):
1316     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1317     break;
1318     case(ContactElementsOne):
1319     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1320     break;
1321     case(ReducedContactElementsOne):
1322     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1323     break;
1324     case(DegreesOfFreedom):
1325 jfenwick 3259 temp=escript::Data( arg, escript::function(*this) );
1326 phornby 1628 _temp=temp.getDataC();
1327     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1328     break;
1329     case(ReducedDegreesOfFreedom):
1330 jfenwick 3259 temp=escript::Data( arg, escript::function(*this) );
1331 phornby 1628 _temp=temp.getDataC();
1332     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1333     break;
1334     default:
1335     stringstream temp;
1336     temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1337     throw FinleyAdapterException(temp.str());
1338     break;
1339     }
1340     checkFinleyError();
1341     blocktimer_increment("integrate()", blocktimer_start);
1342 jgs 82 }
1343 jgs 149
1344 jgs 82 //
1345     // calculates the gradient of arg:
1346     //
1347 jgs 480 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1348 jgs 82 {
1349 jfenwick 1872 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1350 phornby 1628 if (argDomain!=*this)
1351     throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1352 jfenwick 1872 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1353 phornby 1628 if (gradDomain!=*this)
1354     throw FinleyAdapterException("Error - Illegal domain of gradient");
1355 jgs 82
1356 phornby 1628 Finley_Mesh* mesh=m_finleyMesh.get();
1357     escriptDataC _grad=grad.getDataC();
1358     escriptDataC nodeDataC;
1359     escript::Data temp;
1360     if (getMPISize()>1) {
1361 ksteube 1312 if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1362 jfenwick 3259 temp=escript::Data( arg, continuousFunction(*this) );
1363 phornby 1628 nodeDataC = temp.getDataC();
1364 ksteube 1312 } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1365 jfenwick 3259 temp=escript::Data( arg, reducedContinuousFunction(*this) );
1366 phornby 1628 nodeDataC = temp.getDataC();
1367 ksteube 1312 } else {
1368 phornby 1628 nodeDataC = arg.getDataC();
1369 ksteube 1312 }
1370 phornby 1628 } else {
1371     nodeDataC = arg.getDataC();
1372     }
1373     switch(grad.getFunctionSpace().getTypeCode()) {
1374     case(Nodes):
1375     throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1376     break;
1377     case(ReducedNodes):
1378     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1379     break;
1380     case(Elements):
1381     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1382     break;
1383     case(ReducedElements):
1384     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1385     break;
1386     case(FaceElements):
1387     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1388     break;
1389     case(ReducedFaceElements):
1390     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1391     break;
1392     case(Points):
1393     throw FinleyAdapterException("Error - Gradient at points is not supported.");
1394     break;
1395     case(ContactElementsZero):
1396     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1397     break;
1398     case(ReducedContactElementsZero):
1399     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1400     break;
1401     case(ContactElementsOne):
1402     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1403     break;
1404     case(ReducedContactElementsOne):
1405     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1406     break;
1407     case(DegreesOfFreedom):
1408     throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1409     break;
1410     case(ReducedDegreesOfFreedom):
1411     throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1412     break;
1413     default:
1414     stringstream temp;
1415     temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1416     throw FinleyAdapterException(temp.str());
1417     break;
1418     }
1419     checkFinleyError();
1420 jgs 82 }
1421 jgs 149
1422 jgs 82 //
1423     // returns the size of elements:
1424     //
1425 jgs 480 void MeshAdapter::setToSize(escript::Data& size) const
1426 jgs 82 {
1427 phornby 1628 Finley_Mesh* mesh=m_finleyMesh.get();
1428     escriptDataC tmp=size.getDataC();
1429     switch(size.getFunctionSpace().getTypeCode()) {
1430     case(Nodes):
1431     throw FinleyAdapterException("Error - Size of nodes is not supported.");
1432     break;
1433     case(ReducedNodes):
1434     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1435     break;
1436     case(Elements):
1437     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1438     break;
1439     case(ReducedElements):
1440     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1441     break;
1442     case(FaceElements):
1443     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1444     break;
1445     case(ReducedFaceElements):
1446     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1447     break;
1448     case(Points):
1449     throw FinleyAdapterException("Error - Size of point elements is not supported.");
1450     break;
1451     case(ContactElementsZero):
1452     case(ContactElementsOne):
1453     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1454     break;
1455     case(ReducedContactElementsZero):
1456     case(ReducedContactElementsOne):
1457     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1458     break;
1459     case(DegreesOfFreedom):
1460     throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1461     break;
1462     case(ReducedDegreesOfFreedom):
1463     throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1464     break;
1465     default:
1466     stringstream temp;
1467     temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1468     throw FinleyAdapterException(temp.str());
1469     break;
1470     }
1471     checkFinleyError();
1472 jgs 82 }
1473 jgs 149
1474 caltinay 2151 //
1475     // sets the location of nodes
1476     //
1477 jgs 480 void MeshAdapter::setNewX(const escript::Data& new_x)
1478 jgs 82 {
1479 phornby 1628 Finley_Mesh* mesh=m_finleyMesh.get();
1480     escriptDataC tmp;
1481 jfenwick 1872 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1482 phornby 1628 if (newDomain!=*this)
1483     throw FinleyAdapterException("Error - Illegal domain of new point locations");
1484 jfenwick 3259 if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1485 gross 2518 tmp = new_x.getDataC();
1486     Finley_Mesh_setCoordinates(mesh,&tmp);
1487     } else {
1488 jfenwick 3259 escript::Data new_x_inter=escript::Data( new_x, continuousFunction(*this) );
1489 gross 2518 tmp = new_x_inter.getDataC();
1490     Finley_Mesh_setCoordinates(mesh,&tmp);
1491     }
1492 phornby 1628 checkFinleyError();
1493 jgs 82 }
1494 jgs 149
1495 caltinay 2151 //
1496     // Helper for the save* methods. Extracts optional data variable names and the
1497     // corresponding pointers from python dictionary. Caller must free arrays.
1498     //
1499     void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1500 jgs 153 {
1501 caltinay 2151 numData = boost::python::extract<int>(arg.attr("__len__")());
1502 phornby 1628 /* win32 refactor */
1503 caltinay 2151 names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1504     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1505     dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1506 jgs 153
1507 phornby 1628 boost::python::list keys=arg.keys();
1508 caltinay 2151 for (int i=0; i<numData; ++i) {
1509 caltinay 2778 string n=boost::python::extract<string>(keys[i]);
1510 phornby 1628 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1511 jfenwick 1872 if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1512 caltinay 2151 throw FinleyAdapterException("Error: Data must be defined on same Domain");
1513     data[i] = d.getDataC();
1514     dataPtr[i] = &(data[i]);
1515 caltinay 2150 names[i] = TMPMEMALLOC(n.length()+1, char);
1516     strcpy(names[i], n.c_str());
1517 phornby 1628 }
1518 caltinay 2151 }
1519    
1520     //
1521     // saves mesh and optionally data arrays in openDX format
1522     //
1523 caltinay 2778 void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1524 caltinay 2151 {
1525     int num_data;
1526     char **names;
1527     escriptDataC *data;
1528     escriptDataC **ptr_data;
1529    
1530     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1531     Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1532 phornby 1628 checkFinleyError();
1533    
1534     /* win32 refactor */
1535     TMPMEMFREE(data);
1536     TMPMEMFREE(ptr_data);
1537 caltinay 2150 for(int i=0; i<num_data; i++)
1538 phornby 1628 {
1539     TMPMEMFREE(names[i]);
1540     }
1541     TMPMEMFREE(names);
1542 woo409 757
1543 phornby 1628 return;
1544 jgs 82 }
1545 jgs 149
1546 caltinay 2151 //
1547     // saves mesh and optionally data arrays in VTK format
1548     //
1549 caltinay 2778 void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg, const string& metadata, const string& metadata_schema) const
1550 jgs 153 {
1551 caltinay 3355 boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1552     pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1553     metadata, metadata_schema, arg);
1554 caltinay 2151 }
1555 woo409 757
1556 jfenwick 2642 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1557     {
1558 jfenwick 3259 #ifdef ESYS_MPI
1559 jfenwick 2642 index_t myFirstNode=0, myLastNode=0, k=0;
1560     index_t* globalNodeIndex=0;
1561     Finley_Mesh* mesh_p=m_finleyMesh.get();
1562     if (fs_code == FINLEY_REDUCED_NODES)
1563     {
1564     myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1565     myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1566     globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1567     }
1568     else
1569     {
1570     myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1571     myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1572     globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1573     }
1574     k=globalNodeIndex[id];
1575     return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1576 jfenwick 2644 #endif
1577     return true;
1578 jfenwick 2642 }
1579    
1580    
1581    
1582 caltinay 2151 //
1583     // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1584     //
1585 jfenwick 3259 ASM_ptr MeshAdapter::newSystemMatrix(
1586 phornby 1628 const int row_blocksize,
1587     const escript::FunctionSpace& row_functionspace,
1588     const int column_blocksize,
1589     const escript::FunctionSpace& column_functionspace,
1590     const int type) const
1591 jgs 82 {
1592 phornby 1628 int reduceRowOrder=0;
1593     int reduceColOrder=0;
1594     // is the domain right?
1595 jfenwick 1872 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1596 phornby 1628 if (row_domain!=*this)
1597     throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1598 jfenwick 1872 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1599 phornby 1628 if (col_domain!=*this)
1600     throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1601     // is the function space type right
1602     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1603     reduceRowOrder=0;
1604     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1605     reduceRowOrder=1;
1606     } else {
1607     throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1608     }
1609     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1610     reduceColOrder=0;
1611     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1612     reduceColOrder=1;
1613     } else {
1614     throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1615     }
1616     // generate matrix:
1617    
1618     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1619     checkFinleyError();
1620     Paso_SystemMatrix* fsystemMatrix;
1621     int trilinos = 0;
1622     if (trilinos) {
1623 ksteube 1312 #ifdef TRILINOS
1624     /* Allocation Epetra_VrbMatrix here */
1625     #endif
1626 phornby 1628 }
1627     else {
1628 gross 2551 fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1629 phornby 1628 }
1630     checkPasoError();
1631     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1632 jfenwick 3259 SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1633     return ASM_ptr(sma);
1634     // return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1635 jgs 82 }
1636 caltinay 2151
1637     //
1638 gross 1366 // creates a TransportProblemAdapter
1639 caltinay 2151 //
1640 jfenwick 3259 ATP_ptr MeshAdapter::newTransportProblem(
1641 gross 2987 const bool useBackwardEuler,
1642 phornby 1628 const int blocksize,
1643     const escript::FunctionSpace& functionspace,
1644     const int type) const
1645 gross 1366 {
1646 phornby 1628 int reduceOrder=0;
1647     // is the domain right?
1648 jfenwick 1872 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1649 phornby 1628 if (domain!=*this)
1650     throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1651     // is the function space type right
1652     if (functionspace.getTypeCode()==DegreesOfFreedom) {
1653     reduceOrder=0;
1654     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1655     reduceOrder=1;
1656     } else {
1657     throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1658     }
1659     // generate matrix:
1660    
1661     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1662     checkFinleyError();
1663 gross 2987 Paso_TransportProblem* transportProblem;
1664     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1665 phornby 1628 checkPasoError();
1666     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1667 jfenwick 3259 TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1668     return ATP_ptr(tpa);
1669     // 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     temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1701     throw FinleyAdapterException(temp.str());
1702     break;
1703     }
1704     checkFinleyError();
1705     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     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1874     throw FinleyAdapterException(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     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1897     throw FinleyAdapterException(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     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1969     throw FinleyAdapterException(temp.str());
1970     }
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     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1992     throw FinleyAdapterException(temp.str());
1993     }
1994     break;
1995 phornby 1628 default:
1996     stringstream temp;
1997     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1998     throw FinleyAdapterException(temp.str());
1999     break;
2000     }
2001     checkFinleyError();
2002     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     return (m_finleyMesh==temp->m_finleyMesh);
2015     } 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 gross 2315 Finley_Mesh* mesh=m_finleyMesh.get();
2028     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 gross 2315 Finley_Mesh* mesh=m_finleyMesh.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 jfenwick 3259 return continuousFunction(*this).getX();
2043 jgs 102 }
2044 jgs 149
2045 jgs 480 escript::Data MeshAdapter::getNormal() const
2046 jgs 102 {
2047 jfenwick 3259 return functionOnBoundary(*this).getNormal();
2048 jgs 102 }
2049 jgs 149
2050 jgs 480 escript::Data MeshAdapter::getSize() const
2051 jgs 102 {
2052 jfenwick 3259 return escript::function(*this).getSize();
2053 jgs 102 }
2054    
2055 jfenwick 2487 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2056 gross 1062 {
2057 phornby 1628 int *out = NULL;
2058     Finley_Mesh* mesh=m_finleyMesh.get();
2059     switch (functionSpaceType) {
2060     case(Nodes):
2061     out=mesh->Nodes->Id;
2062     break;
2063     case(ReducedNodes):
2064     out=mesh->Nodes->reducedNodesId;
2065     break;
2066     case(Elements):
2067     out=mesh->Elements->Id;
2068     break;
2069     case(ReducedElements):
2070     out=mesh->Elements->Id;
2071     break;
2072     case(FaceElements):
2073     out=mesh->FaceElements->Id;
2074     break;
2075     case(ReducedFaceElements):
2076     out=mesh->FaceElements->Id;
2077     break;
2078     case(Points):
2079     out=mesh->Points->Id;
2080     break;
2081     case(ContactElementsZero):
2082     out=mesh->ContactElements->Id;
2083     break;
2084     case(ReducedContactElementsZero):
2085     out=mesh->ContactElements->Id;
2086     break;
2087     case(ContactElementsOne):
2088     out=mesh->ContactElements->Id;
2089     break;
2090     case(ReducedContactElementsOne):
2091     out=mesh->ContactElements->Id;
2092     break;
2093     case(DegreesOfFreedom):
2094     out=mesh->Nodes->degreesOfFreedomId;
2095     break;
2096     case(ReducedDegreesOfFreedom):
2097     out=mesh->Nodes->reducedDegreesOfFreedomId;
2098     break;
2099     default:
2100 gross 1062 stringstream temp;
2101     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2102     throw FinleyAdapterException(temp.str());
2103     break;
2104 phornby 1628 }
2105     return out;
2106 gross 1062 }
2107 jgs 110 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2108     {
2109 phornby 1628 int out=0;
2110     Finley_Mesh* mesh=m_finleyMesh.get();
2111     switch (functionSpaceType) {
2112     case(Nodes):
2113     out=mesh->Nodes->Tag[sampleNo];
2114     break;
2115     case(ReducedNodes):
2116     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2117     break;
2118     case(Elements):
2119     out=mesh->Elements->Tag[sampleNo];
2120     break;
2121     case(ReducedElements):
2122     out=mesh->Elements->Tag[sampleNo];
2123     break;
2124     case(FaceElements):
2125     out=mesh->FaceElements->Tag[sampleNo];
2126     break;
2127     case(ReducedFaceElements):
2128     out=mesh->FaceElements->Tag[sampleNo];
2129     break;
2130     case(Points):
2131     out=mesh->Points->Tag[sampleNo];
2132     break;
2133     case(ContactElementsZero):
2134     out=mesh->ContactElements->Tag[sampleNo];
2135     break;
2136     case(ReducedContactElementsZero):
2137     out=mesh->ContactElements->Tag[sampleNo];
2138     break;
2139     case(ContactElementsOne):
2140     out=mesh->ContactElements->Tag[sampleNo];
2141     break;
2142     case(ReducedContactElementsOne):
2143     out=mesh->ContactElements->Tag[sampleNo];
2144     break;
2145     case(DegreesOfFreedom):
2146     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2147     break;
2148     case(ReducedDegreesOfFreedom):
2149     throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2150     break;
2151     default:
2152     stringstream temp;
2153     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2154     throw FinleyAdapterException(temp.str());
2155     break;
2156     }
2157     return out;
2158 jgs 110 }
2159    
2160 gross 1062
2161 gross 767 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript