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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3490 - (hide annotations)
Wed Mar 30 02:24:33 2011 UTC (8 years, 1 month ago) by caltinay
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 85874 byte(s)
More gcc-4.6 fixes (mostly initialized-but-unused-var warnings)

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