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

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

Parent Directory Parent Directory | Revision Log Revision Log


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