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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4441 - (hide annotations)
Fri Jun 7 02:23:49 2013 UTC (5 years, 10 months ago) by caltinay
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 91372 byte(s)
finley now uses Data objects directly instead of going through the C wrapper.

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