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

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

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

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

Legend:
Removed from v.3344  
changed lines
  Added in v.4498

  ViewVC Help
Powered by ViewVC 1.1.26