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

Annotation of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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