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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.4346  
changed lines
  Added in v.4949

  ViewVC Help
Powered by ViewVC 1.1.26