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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4428 - (hide annotations)
Thu May 30 06:39:10 2013 UTC (6 years ago) by caltinay
File size: 94529 byte(s)
finley's NodeFile is now a class.
Associated changes:
- use of some C++ constructs/functions/types (1st pass)
- removal of obsolete pointer check
- merging of some duplicated code
- ...

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