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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4408 - (hide annotations)
Tue May 14 06:58:43 2013 UTC (6 years, 4 months ago) by caltinay
File size: 94960 byte(s)
Mostly no-op - trying to get some consistency in style while going through a
few files.
Also some minor changes related to C->C++ (e.g. no need to copy vector to pass
to functions etc.)

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