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

Legend:
Removed from v.3998  
changed lines
  Added in v.4408

  ViewVC Help
Powered by ViewVC 1.1.26