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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4437 - (hide annotations)
Tue Jun 4 05:37:18 2013 UTC (5 years, 9 months ago) by caltinay
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 94260 byte(s)
finley:
-Removed obsolete reference files (OpenDX)
-Converted tag map to a std::map
-loosened tests to ignore tag order in file
-added some decorators for test skips

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