/[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 4410 - (hide annotations)
Tue May 14 11:29:23 2013 UTC (5 years, 11 months ago) by caltinay
File size: 95049 byte(s)
Oops, I saw this coming.

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 4408 return Finley_NodeFile_getGlobalNumNodes(m_finleyMesh.get()->Nodes);
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     numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);
668     break;
669     case ReducedNodes:
670     numDataPointsPerSample=1;
671     numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
672     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     numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
731     }
732     break;
733     case ReducedDegreesOfFreedom:
734     if (mesh->Nodes!=NULL) {
735     numDataPointsPerSample=1;
736     numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
737     }
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     escriptDataC tmp;
1449     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1450     if (newDomain!=*this)
1451     throw FinleyAdapterException("Error - Illegal domain of new point locations");
1452     if (new_x.getFunctionSpace() == continuousFunction(*this)) {
1453     tmp = new_x.getDataC();
1454     Finley_Mesh_setCoordinates(mesh,&tmp);
1455     } else {
1456     throw FinleyAdapterException("As of escript version 3.3 SetX() only accepts ContinuousFunction arguments. Please interpolate.");
1457     //escript::Data new_x_inter=escript::Data( new_x, continuousFunction(*this) );
1458     //tmp = new_x_inter.getDataC();
1459     //Finley_Mesh_setCoordinates(mesh,&tmp);
1460     }
1461     checkFinleyError();
1462 jgs 82 }
1463 jgs 149
1464 jfenwick 2642 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1465     {
1466 caltinay 3993 if (getMPISize() > 1) {
1467 jfenwick 3259 #ifdef ESYS_MPI
1468 caltinay 3993 index_t myFirstNode=0, myLastNode=0, k=0;
1469     index_t* globalNodeIndex=0;
1470     Finley_Mesh* mesh_p=m_finleyMesh.get();
1471 caltinay 3998 /*
1472     * this method is only used by saveDataCSV which would use the returned
1473     * values for reduced nodes wrongly so this case is disabled for now
1474 caltinay 4408 if (fs_code == FINLEY_REDUCED_NODES) {
1475 caltinay 3993 myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1476     myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1477     globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1478 caltinay 4408 } else
1479 caltinay 3998 */
1480 caltinay 4408 if (fs_code == FINLEY_NODES) {
1481 caltinay 3993 myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1482     myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1483     globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1484 caltinay 4408 } else {
1485 caltinay 3993 throw FinleyAdapterException("Unsupported function space type for ownSample()");
1486     }
1487    
1488     k=globalNodeIndex[id];
1489     return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1490     #endif
1491 jfenwick 2642 }
1492 jfenwick 2644 return true;
1493 jfenwick 2642 }
1494    
1495    
1496 caltinay 2151 //
1497     // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1498     //
1499 caltinay 4408 escript::ASM_ptr MeshAdapter::newSystemMatrix(const int row_blocksize,
1500     const escript::FunctionSpace& row_functionspace,
1501     const int column_blocksize,
1502     const escript::FunctionSpace& column_functionspace,
1503     const int type) const
1504 jgs 82 {
1505 caltinay 4408 // is the domain right?
1506     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1507     if (row_domain!=*this)
1508     throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1509     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1510     if (col_domain!=*this)
1511     throw FinleyAdapterException("Error - domain of column function space does not match the domain of matrix generator.");
1512    
1513     int reduceRowOrder=0;
1514     int reduceColOrder=0;
1515     // is the function space type right?
1516     if (row_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1517     reduceRowOrder=1;
1518     } else if (row_functionspace.getTypeCode() != DegreesOfFreedom) {
1519     throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1520     }
1521     if (column_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1522     reduceColOrder=1;
1523     } else if (column_functionspace.getTypeCode() != DegreesOfFreedom) {
1524     throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1525     }
1526    
1527     // generate matrix:
1528     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(
1529     getFinley_Mesh(), reduceRowOrder, reduceColOrder);
1530     checkFinleyError();
1531     Paso_SystemMatrix* fsystemMatrix;
1532     const int trilinos = 0;
1533     if (trilinos) {
1534 ksteube 1312 #ifdef TRILINOS
1535 caltinay 4408 // FIXME: Allocation Epetra_VrbMatrix here...
1536 ksteube 1312 #endif
1537 caltinay 4408 } else {
1538     fsystemMatrix=Paso_SystemMatrix_alloc(type, fsystemMatrixPattern,
1539     row_blocksize, column_blocksize, FALSE);
1540     }
1541     checkPasoError();
1542     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1543     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1544     return escript::ASM_ptr(sma);
1545 jgs 82 }
1546 caltinay 2151
1547     //
1548 gross 1366 // creates a TransportProblemAdapter
1549 caltinay 2151 //
1550 caltinay 4408 escript::ATP_ptr MeshAdapter::newTransportProblem(const int blocksize,
1551     const escript::FunctionSpace& functionspace, const int type) const
1552 gross 1366 {
1553 caltinay 4408 // is the domain right?
1554     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1555     if (domain!=*this)
1556     throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1557    
1558     // is the function space type right?
1559     int reduceOrder=0;
1560     if (functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1561     reduceOrder=1;
1562     } else if (functionspace.getTypeCode() != DegreesOfFreedom) {
1563     throw FinleyAdapterException("Error - illegal function space type for transport problem.");
1564     }
1565    
1566     // generate transport problem:
1567     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(
1568     getFinley_Mesh(), reduceOrder, reduceOrder);
1569     checkFinleyError();
1570     Paso_TransportProblem* transportProblem;
1571     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern, blocksize);
1572     checkPasoError();
1573     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1574     TransportProblemAdapter* tpa=new TransportProblemAdapter(
1575     transportProblem, blocksize, functionspace);
1576     return escript::ATP_ptr(tpa);
1577 gross 1366 }
1578 jgs 149
1579 jgs 82 //
1580 caltinay 4408 // returns true if data on functionSpaceCode is considered as being cell centered
1581 jgs 82 //
1582     bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1583     {
1584 caltinay 4408 switch(functionSpaceCode) {
1585     case Nodes:
1586     case DegreesOfFreedom:
1587     case ReducedDegreesOfFreedom:
1588     return false;
1589     case Elements:
1590     case FaceElements:
1591     case Points:
1592     case ContactElementsZero:
1593     case ContactElementsOne:
1594     case ReducedElements:
1595     case ReducedFaceElements:
1596     case ReducedContactElementsZero:
1597     case ReducedContactElementsOne:
1598     return true;
1599     default:
1600     stringstream temp;
1601     temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1602     throw FinleyAdapterException(temp.str());
1603     break;
1604     }
1605     return false;
1606 jgs 82 }
1607 jgs 149
1608 jfenwick 2635 bool
1609 caltinay 2778 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1610 jfenwick 2635 {
1611 caltinay 4408 /* The idea is to use equivalence classes. [Types which can be interpolated
1612     back and forth]:
1613     class 1: DOF <-> Nodes
1614     class 2: ReducedDOF <-> ReducedNodes
1615     class 3: Points
1616     class 4: Elements
1617     class 5: ReducedElements
1618     class 6: FaceElements
1619     class 7: ReducedFaceElements
1620     class 8: ContactElementZero <-> ContactElementOne
1621     class 9: ReducedContactElementZero <-> ReducedContactElementOne
1622 jfenwick 2635
1623 caltinay 4408 There is also a set of lines. Interpolation is possible down a line but
1624     not between lines.
1625     class 1 and 2 belong to all lines so aren't considered.
1626     line 0: class 3
1627     line 1: class 4,5
1628     line 2: class 6,7
1629     line 3: class 8,9
1630 jfenwick 2635
1631 caltinay 4408 For classes with multiple members (eg class 2) we have vars to record if
1632     there is at least one instance.
1633     e.g. hasnodes is true if we have at least one instance of Nodes.
1634     */
1635 caltinay 2940 if (fs.empty())
1636 caltinay 2778 return false;
1637 caltinay 4408
1638 caltinay 2778 vector<int> hasclass(10);
1639 caltinay 4408 vector<int> hasline(4);
1640 jfenwick 2635 bool hasnodes=false;
1641     bool hasrednodes=false;
1642     bool hascez=false;
1643     bool hasrcez=false;
1644 caltinay 4408 for (int i=0;i<fs.size();++i) {
1645     switch(fs[i]) {
1646     case Nodes:
1647     hasnodes=true; // fall through
1648     case DegreesOfFreedom:
1649     hasclass[1]=1;
1650     break;
1651     case ReducedNodes:
1652     hasrednodes=true; // fall through
1653     case ReducedDegreesOfFreedom:
1654     hasclass[2]=1;
1655     break;
1656     case Points:
1657     hasline[0]=1;
1658     hasclass[3]=1;
1659     break;
1660     case Elements:
1661     hasclass[4]=1;
1662     hasline[1]=1;
1663     break;
1664     case ReducedElements:
1665     hasclass[5]=1;
1666     hasline[1]=1;
1667     break;
1668     case FaceElements:
1669     hasclass[6]=1;
1670     hasline[2]=1;
1671     break;
1672     case ReducedFaceElements:
1673     hasclass[7]=1;
1674     hasline[2]=1;
1675     break;
1676     case ContactElementsZero:
1677     hascez=true; // fall through
1678     case ContactElementsOne:
1679     hasclass[8]=1;
1680     hasline[3]=1;
1681     break;
1682     case ReducedContactElementsZero:
1683     hasrcez=true; // fall through
1684     case ReducedContactElementsOne:
1685     hasclass[9]=1;
1686     hasline[3]=1;
1687     break;
1688     default:
1689     return false;
1690     }
1691 jfenwick 2635 }
1692     int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1693 caltinay 4408
1694 jfenwick 2635 // fail if we have more than one leaf group
1695     if (totlines>1)
1696 caltinay 4408 return false; // there are at least two branches we can't interpolate between
1697     else if (totlines==1) {
1698     if (hasline[0]==1) // we have points
1699     resultcode=Points;
1700     else if (hasline[1]==1) {
1701     if (hasclass[5]==1)
1702     resultcode=ReducedElements;
1703     else
1704     resultcode=Elements;
1705     } else if (hasline[2]==1) {
1706     if (hasclass[7]==1)
1707     resultcode=ReducedFaceElements;
1708     else
1709     resultcode=FaceElements;
1710     } else { // so we must be in line3
1711     if (hasclass[9]==1) {
1712     // need something from class 9
1713     resultcode=(hasrcez ? ReducedContactElementsZero : ReducedContactElementsOne);
1714     } else {
1715     // something from class 8
1716     resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1717     }
1718     }
1719     } else { // totlines==0
1720     if (hasclass[2]==1) {
1721     // something from class 2
1722     resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
1723     } else {
1724     // something from class 1
1725     resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
1726     }
1727 jfenwick 2635 }
1728     return true;
1729     }
1730    
1731 caltinay 4408 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,
1732     int functionSpaceType_target) const
1733 jgs 82 {
1734 caltinay 4408 switch(functionSpaceType_source) {
1735     case Nodes:
1736     switch(functionSpaceType_target) {
1737     case Nodes:
1738     case ReducedNodes:
1739     case ReducedDegreesOfFreedom:
1740     case DegreesOfFreedom:
1741     case Elements:
1742     case ReducedElements:
1743     case FaceElements:
1744     case ReducedFaceElements:
1745     case Points:
1746     case ContactElementsZero:
1747     case ReducedContactElementsZero:
1748     case ContactElementsOne:
1749     case ReducedContactElementsOne:
1750     return true;
1751     default:
1752     stringstream temp;
1753     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1754     throw FinleyAdapterException(temp.str());
1755     }
1756     break;
1757     case ReducedNodes:
1758     switch(functionSpaceType_target) {
1759     case ReducedNodes:
1760     case ReducedDegreesOfFreedom:
1761     case Elements:
1762     case ReducedElements:
1763     case FaceElements:
1764     case ReducedFaceElements:
1765     case Points:
1766     case ContactElementsZero:
1767     case ReducedContactElementsZero:
1768     case ContactElementsOne:
1769     case ReducedContactElementsOne:
1770     return true;
1771     case Nodes:
1772     case DegreesOfFreedom:
1773     return false;
1774     default:
1775     stringstream temp;
1776     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1777     throw FinleyAdapterException(temp.str());
1778     }
1779     break;
1780     case Elements:
1781     if (functionSpaceType_target==Elements) {
1782     return true;
1783     } else if (functionSpaceType_target==ReducedElements) {
1784     return true;
1785     } else {
1786     return false;
1787     }
1788     case ReducedElements:
1789     if (functionSpaceType_target==ReducedElements) {
1790     return true;
1791     } else {
1792     return false;
1793     }
1794     case FaceElements:
1795     if (functionSpaceType_target==FaceElements) {
1796     return true;
1797     } else if (functionSpaceType_target==ReducedFaceElements) {
1798     return true;
1799     } else {
1800     return false;
1801     }
1802     case ReducedFaceElements:
1803     if (functionSpaceType_target==ReducedFaceElements) {
1804     return true;
1805     } else {
1806     return false;
1807     }
1808     case Points:
1809     if (functionSpaceType_target==Points) {
1810     return true;
1811     } else {
1812     return false;
1813     }
1814     case ContactElementsZero:
1815     case ContactElementsOne:
1816     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1817     return true;
1818     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1819     return true;
1820     } else {
1821     return false;
1822     }
1823     case ReducedContactElementsZero:
1824     case ReducedContactElementsOne:
1825     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1826     return true;
1827     } else {
1828     return false;
1829     }
1830     case DegreesOfFreedom:
1831     switch(functionSpaceType_target) {
1832     case ReducedDegreesOfFreedom:
1833     case DegreesOfFreedom:
1834     case Nodes:
1835     case ReducedNodes:
1836     case Elements:
1837     case ReducedElements:
1838     case Points:
1839     case FaceElements:
1840     case ReducedFaceElements:
1841     case ContactElementsZero:
1842     case ReducedContactElementsZero:
1843     case ContactElementsOne:
1844     case ReducedContactElementsOne:
1845     return true;
1846     default:
1847     stringstream temp;
1848     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1849     throw FinleyAdapterException(temp.str());
1850     }
1851     break;
1852     case ReducedDegreesOfFreedom:
1853     switch(functionSpaceType_target) {
1854     case ReducedDegreesOfFreedom:
1855     case ReducedNodes:
1856     case Elements:
1857     case ReducedElements:
1858     case FaceElements:
1859     case ReducedFaceElements:
1860     case Points:
1861     case ContactElementsZero:
1862     case ReducedContactElementsZero:
1863     case ContactElementsOne:
1864     case ReducedContactElementsOne:
1865     return true;
1866     case Nodes:
1867     case DegreesOfFreedom:
1868     return false;
1869     default:
1870     stringstream temp;
1871     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1872     throw FinleyAdapterException(temp.str());
1873     }
1874     break;
1875     default:
1876     stringstream temp;
1877     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1878     throw FinleyAdapterException(temp.str());
1879     break;
1880     }
1881     return false;
1882 jgs 82 }
1883 jgs 149
1884 caltinay 4408 signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source, int functionSpaceType_target) const
1885 jfenwick 4255 {
1886     if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1887     return 1;
1888 caltinay 4408
1889     if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1890 jfenwick 4255 return -1;
1891 caltinay 4408
1892 jfenwick 4255 return 0;
1893 caltinay 4408 }
1894    
1895     bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,
1896     const escript::AbstractDomain& targetDomain,
1897     int functionSpaceType_target) const
1898 jgs 82 {
1899 caltinay 4408 return false;
1900 jgs 82 }
1901 jgs 104
1902 caltinay 4408 bool MeshAdapter::operator==(const escript::AbstractDomain& other) const
1903 jgs 82 {
1904 caltinay 4408 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1905     if (temp) {
1906     return (m_finleyMesh==temp->m_finleyMesh);
1907     } else {
1908     return false;
1909     }
1910 jgs 82 }
1911    
1912 caltinay 4408 bool MeshAdapter::operator!=(const escript::AbstractDomain& other) const
1913 jgs 104 {
1914 caltinay 4408 return !(operator==(other));
1915 jgs 104 }
1916    
1917 gross 1859 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1918 jgs 102 {
1919 caltinay 4408 Finley_Mesh* mesh=m_finleyMesh.get();
1920     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
1921     package, symmetry, mesh->MPIInfo);
1922 jgs 102 }
1923 caltinay 3681
1924 gross 1859 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1925     {
1926 caltinay 4408 Finley_Mesh* mesh=m_finleyMesh.get();
1927     return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
1928     package, symmetry, mesh->MPIInfo);
1929 gross 1859 }
1930 jgs 149
1931 jgs 480 escript::Data MeshAdapter::getX() const
1932 jgs 102 {
1933 caltinay 4408 return continuousFunction(*this).getX();
1934 jgs 102 }
1935 jgs 149
1936 jgs 480 escript::Data MeshAdapter::getNormal() const
1937 jgs 102 {
1938 caltinay 4408 return functionOnBoundary(*this).getNormal();
1939 jgs 102 }
1940 jgs 149
1941 jgs 480 escript::Data MeshAdapter::getSize() const
1942 jgs 102 {
1943 caltinay 4408 return escript::function(*this).getSize();
1944 jgs 102 }
1945    
1946 jfenwick 2487 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1947 gross 1062 {
1948 caltinay 4408 int *out = NULL;
1949     Finley_Mesh* mesh=m_finleyMesh.get();
1950     switch (functionSpaceType) {
1951     case Nodes:
1952     out=mesh->Nodes->Id;
1953     break;
1954     case ReducedNodes:
1955     out=mesh->Nodes->reducedNodesId;
1956     break;
1957     case Elements:
1958     case ReducedElements:
1959     out=mesh->Elements->Id;
1960     break;
1961     case FaceElements:
1962     case ReducedFaceElements:
1963     out=mesh->FaceElements->Id;
1964     break;
1965     case Points:
1966     out=mesh->Points->Id;
1967     break;
1968     case ContactElementsZero:
1969     case ReducedContactElementsZero:
1970     case ContactElementsOne:
1971     case ReducedContactElementsOne:
1972     out=mesh->ContactElements->Id;
1973     break;
1974     case DegreesOfFreedom:
1975     out=mesh->Nodes->degreesOfFreedomId;
1976     break;
1977     case ReducedDegreesOfFreedom:
1978     out=mesh->Nodes->reducedDegreesOfFreedomId;
1979     break;
1980     default:
1981     stringstream temp;
1982     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1983     throw FinleyAdapterException(temp.str());
1984     break;
1985     }
1986     return out;
1987 gross 1062 }
1988 jgs 110 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1989     {
1990 caltinay 4408 int out=0;
1991     Finley_Mesh* mesh=m_finleyMesh.get();
1992     switch (functionSpaceType) {
1993     case Nodes:
1994     out=mesh->Nodes->Tag[sampleNo];
1995     break;
1996     case ReducedNodes:
1997     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1998     break;
1999     case Elements:
2000     case ReducedElements:
2001     out=mesh->Elements->Tag[sampleNo];
2002     break;
2003     case FaceElements:
2004     case ReducedFaceElements:
2005     out=mesh->FaceElements->Tag[sampleNo];
2006     break;
2007     case Points:
2008     out=mesh->Points->Tag[sampleNo];
2009     break;
2010     case ContactElementsZero:
2011     case ReducedContactElementsZero:
2012     case ContactElementsOne:
2013     case ReducedContactElementsOne:
2014     out=mesh->ContactElements->Tag[sampleNo];
2015     break;
2016     case DegreesOfFreedom:
2017     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2018     break;
2019     case ReducedDegreesOfFreedom:
2020     throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2021     break;
2022     default:
2023     stringstream temp;
2024     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2025     throw FinleyAdapterException(temp.str());
2026     break;
2027     }
2028     return out;
2029 jgs 110 }
2030    
2031 gross 1062
2032 gross 767 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2033     {
2034 caltinay 4408 Finley_Mesh* mesh=m_finleyMesh.get();
2035     escriptDataC tmp=mask.getDataC();
2036     switch(functionSpaceType) {
2037     case Nodes:
2038     Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2039     break;
2040     case ReducedNodes:
2041     throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2042     break;
2043     case DegreesOfFreedom:
2044     throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2045     break;
2046     case ReducedDegreesOfFreedom:
2047     throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2048     break;
2049     case Elements:
2050     case ReducedElements:
2051     Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2052     break;
2053     case FaceElements:
2054     case ReducedFaceElements:
2055     Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2056     break;
2057     case Points:
2058     Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2059     break;
2060     case ContactElementsZero:
2061     case ReducedContactElementsZero:
2062     case ContactElementsOne:
2063     case ReducedContactElementsOne:
2064     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2065     break;
2066     default:
2067     stringstream temp;
2068     temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2069     throw FinleyAdapterException(temp.str());
2070     }
2071     checkFinleyError();
2072 gross 767 }
2073    
2074 caltinay 4408 void MeshAdapter::setTagMap(const string& name, int tag)
2075 gross 1044 {
2076 caltinay 4408 Finley_Mesh* mesh=m_finleyMesh.get();
2077     Finley_Mesh_addTagMap(mesh, name.c_str(), tag);
2078     checkFinleyError();
2079 gross 1044 }
2080 gross 767
2081 caltinay 2778 int MeshAdapter::getTag(const string& name) const
2082 gross 1044 {
2083 caltinay 4408 Finley_Mesh* mesh=m_finleyMesh.get();
2084     int tag = Finley_Mesh_getTag(mesh, name.c_str());
2085     checkFinleyError();
2086     return tag;
2087 gross 1044 }
2088    
2089 caltinay 2778 bool MeshAdapter::isValidTagName(const string& name) const
2090 gross 1044 {
2091 caltinay 4408 Finley_Mesh* mesh=m_finleyMesh.get();
2092     return Finley_Mesh_isValidTagName(mesh, name.c_str());
2093 gross 1044 }
2094    
2095 caltinay 2778 string MeshAdapter::showTagNames() const
2096 gross 1044 {
2097 caltinay 4408 stringstream temp;
2098     Finley_Mesh* mesh=m_finleyMesh.get();
2099     Finley_TagMap* tag_map=mesh->TagMap;
2100     while (tag_map) {
2101     temp << tag_map->name;
2102     tag_map=tag_map->next;
2103     if (tag_map) temp << ", ";
2104     }
2105     return temp.str();
2106 gross 1044 }
2107    
2108 gross 1716 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2109     {
2110 caltinay 4408 Finley_Mesh* mesh=m_finleyMesh.get();
2111     dim_t numTags=0;
2112     switch(functionSpaceCode) {
2113     case Nodes:
2114     numTags=mesh->Nodes->numTagsInUse;
2115     break;
2116     case ReducedNodes:
2117     throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2118     break;
2119     case DegreesOfFreedom:
2120     throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2121     break;
2122     case ReducedDegreesOfFreedom:
2123     throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2124     break;
2125     case Elements:
2126     case ReducedElements:
2127     numTags=mesh->Elements->numTagsInUse;
2128     break;
2129     case FaceElements:
2130     case ReducedFaceElements:
2131     numTags=mesh->FaceElements->numTagsInUse;
2132     break;
2133     case Points:
2134     numTags=mesh->Points->numTagsInUse;
2135     break;
2136     case ContactElementsZero:
2137     case ReducedContactElementsZero:
2138     case ContactElementsOne:
2139     case ReducedContactElementsOne:
2140     numTags=mesh->ContactElements->numTagsInUse;
2141     break;
2142     default:
2143     stringstream temp;
2144     temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2145     throw FinleyAdapterException(temp.str());
2146     }
2147     return numTags;
2148 gross 1716 }
2149 jfenwick 2487
2150     const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2151 gross 1716 {
2152 caltinay 4408 Finley_Mesh* mesh=m_finleyMesh.get();
2153     index_t* tags=NULL;
2154     switch(functionSpaceCode) {
2155     case Nodes:
2156     tags=mesh->Nodes->tagsInUse;
2157     break;
2158     case ReducedNodes:
2159     throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2160     break;
2161     case DegreesOfFreedom:
2162     throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2163     break;
2164     case ReducedDegreesOfFreedom:
2165     throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2166     break;
2167     case Elements:
2168     case ReducedElements:
2169     tags=mesh->Elements->tagsInUse;
2170     break;
2171     case FaceElements:
2172     case ReducedFaceElements:
2173     tags=mesh->FaceElements->tagsInUse;
2174     break;
2175     case Points:
2176     tags=mesh->Points->tagsInUse;
2177     break;
2178     case ContactElementsZero:
2179     case ReducedContactElementsZero:
2180     case ContactElementsOne:
2181     case ReducedContactElementsOne:
2182     tags=mesh->ContactElements->tagsInUse;
2183     break;
2184     default:
2185     stringstream temp;
2186     temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2187     throw FinleyAdapterException(temp.str());
2188     }
2189     return tags;
2190 gross 1716 }
2191    
2192    
2193 jfenwick 1802 bool MeshAdapter::canTag(int functionSpaceCode) const
2194     {
2195 caltinay 4408 switch(functionSpaceCode) {
2196     case Nodes:
2197     case Elements:
2198     case ReducedElements:
2199     case FaceElements:
2200     case ReducedFaceElements:
2201     case Points:
2202     case ContactElementsZero:
2203     case ReducedContactElementsZero:
2204     case ContactElementsOne:
2205     case ReducedContactElementsOne:
2206     return true;
2207     default:
2208     return false;
2209     }
2210 jfenwick 1802 }
2211    
2212 caltinay 4408 escript::AbstractDomain::StatusType MeshAdapter::getStatus() const
2213 gross 2533 {
2214 caltinay 4408 Finley_Mesh* mesh=m_finleyMesh.get();
2215     return Finley_Mesh_getStatus(mesh);
2216 gross 2533 }
2217 jfenwick 1802
2218 gross 2856 int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2219     {
2220 caltinay 4408 Finley_Mesh* mesh=m_finleyMesh.get();
2221     int order =-1;
2222     switch(functionSpaceCode) {
2223     case Nodes:
2224     case DegreesOfFreedom:
2225     order=mesh->approximationOrder;
2226     break;
2227     case ReducedNodes:
2228     case ReducedDegreesOfFreedom:
2229     order=mesh->reducedApproximationOrder;
2230     break;
2231     case Elements:
2232     case FaceElements:
2233     case Points:
2234     case ContactElementsZero:
2235     case ContactElementsOne:
2236     order=mesh->integrationOrder;
2237     break;
2238     case ReducedElements:
2239     case ReducedFaceElements:
2240     case ReducedContactElementsZero:
2241     case ReducedContactElementsOne:
2242     order=mesh->reducedIntegrationOrder;
2243     break;
2244     default:
2245     stringstream temp;
2246     temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2247     throw FinleyAdapterException(temp.str());
2248     }
2249     return order;
2250 gross 2856 }
2251 gross 2533
2252 jfenwick 3259 bool MeshAdapter::supportsContactElements() const
2253 caltinay 2898 {
2254 caltinay 4408 return true;
2255 jfenwick 3259 }
2256    
2257 caltinay 4408 ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id,
2258     index_t order, index_t reducedOrder)
2259 jfenwick 3259 {
2260 caltinay 4408 m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2261 caltinay 2898 }
2262    
2263     ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2264     {
2265 caltinay 4408 Finley_ReferenceElementSet_dealloc(m_refSet);
2266 caltinay 2898 }