/[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 2989 by gross, Thu Mar 18 01:45:55 2010 UTC revision 4428 by caltinay, Thu May 30 06:39:10 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2010 by University of Queensland  * Copyright (c) 2003-2013 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16    #include <pasowrap/PasoException.h>
17    #include <pasowrap/TransportProblemAdapter.h>
18  #include "MeshAdapter.h"  #include "MeshAdapter.h"
19  #include "escript/Data.h"  #include "escript/Data.h"
20  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
21    #include "esysUtils/blocktimer.h"
22    
23    #include <boost/python/import.hpp>
24  #ifdef USE_NETCDF  #ifdef USE_NETCDF
25  #include <netcdfcpp.h>  #include <netcdfcpp.h>
26  #endif  #endif
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #include "paso/Paso_MPI.h"  
 #endif  
 extern "C" {  
 #include "esysUtils/blocktimer.h"  
 }  
27    
28  using namespace std;  using namespace std;
29  using namespace escript;  using namespace paso;
30    
31  namespace finley {  namespace finley {
32    
 //  
33  // define the static constants  // define the static constants
34  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
35  const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;  const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;
# Line 50  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 PASO_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
96  #ifdef PASO_MPI  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 PASO_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 PASO_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 PASO_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 = Paso_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 PASO_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::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
632  {  {
633     return Points;      return Points;
634  }  }
635    
636  //  //
# Line 659  int MeshAdapter::getDiracDeltaFunctionCo Line 638  int MeshAdapter::getDiracDeltaFunctionCo
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 670  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 m_finleyMesh.get()->Nodes->getGlobalNumNodes();
653  }  }
654    
655  //  //
# Line 679  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=mesh->Nodes->getNumNodes();
668     break;              break;
669     case(ReducedNodes):          case ReducedNodes:
670     numDataPointsPerSample=1;              numDataPointsPerSample=1;
671     numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);              numSamples=mesh->Nodes->getNumReducedNodes();
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=mesh->Nodes->getNumDegreesOfFreedom();
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=mesh->Nodes->getNumReducedDegreesOfFreedom();
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                                   SystemMatrixAdapter& 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          const escript::Data& d, const escript::Data& y,
757  {          const escript::Data& d_contact, const escript::Data& y_contact,
758     escriptDataC _rhs=rhs.getDataC();          const escript::Data& d_dirac, const escript::Data& y_dirac) const
759     escriptDataC _A  =A.getDataC();  {
760     escriptDataC _B=B.getDataC();      SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
761     escriptDataC _C=C.getDataC();      if (!smat)
762     escriptDataC _D=D.getDataC();          throw FinleyAdapterException("finley only supports Paso system matrices.");
763     escriptDataC _X=X.getDataC();  
764     escriptDataC _Y=Y.getDataC();      escriptDataC _rhs=rhs.getDataC();
765     escriptDataC _d=d.getDataC();      escriptDataC _A  =A.getDataC();
766     escriptDataC _y=y.getDataC();      escriptDataC _B=B.getDataC();
767     escriptDataC _d_contact=d_contact.getDataC();      escriptDataC _C=C.getDataC();
768     escriptDataC _y_contact=y_contact.getDataC();      escriptDataC _D=D.getDataC();
769        escriptDataC _X=X.getDataC();
770     Finley_Mesh* mesh=m_finleyMesh.get();      escriptDataC _Y=Y.getDataC();
771        escriptDataC _d=d.getDataC();
772     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );      escriptDataC _y=y.getDataC();
773     checkFinleyError();      escriptDataC _d_contact=d_contact.getDataC();
774        escriptDataC _y_contact=y_contact.getDataC();
775     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );      escriptDataC _d_dirac=d_dirac.getDataC();
776     checkFinleyError();      escriptDataC _y_dirac=y_dirac.getDataC();
777    
778     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );      Finley_Mesh* mesh=m_finleyMesh.get();
779     checkFinleyError();      Paso_SystemMatrix* S = smat->getPaso_SystemMatrix();
780  }  
781        Finley_Assemble_PDE(mesh->Nodes, mesh->Elements, S, &_rhs, &_A, &_B, &_C,
782  void  MeshAdapter::addPDEToLumpedSystem(                          &_D, &_X, &_Y);
783                                          escript::Data& mat,      checkFinleyError();
784                                          const escript::Data& D,  
785                                          const escript::Data& d) const      Finley_Assemble_PDE(mesh->Nodes, mesh->FaceElements, S, &_rhs, 0, 0, 0,
786  {                          &_d, 0, &_y);
787     escriptDataC _mat=mat.getDataC();      checkFinleyError();
788     escriptDataC _D=D.getDataC();  
789     escriptDataC _d=d.getDataC();      Finley_Assemble_PDE(mesh->Nodes, mesh->ContactElements, S, &_rhs, 0, 0, 0,
790                &_d_contact, 0, &_y_contact);
791     Finley_Mesh* mesh=m_finleyMesh.get();      checkFinleyError();
792    
793     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);      Finley_Assemble_PDE(mesh->Nodes, mesh->Points, S, &_rhs, 0, 0, 0,
794     checkFinleyError();                          &_d_dirac, 0, &_y_dirac );
795          checkFinleyError();
796     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);  }
797     checkFinleyError();  
798    void MeshAdapter::addPDEToLumpedSystem(escript::Data& mat,
799            const escript::Data& D, const escript::Data& d,
800            const escript::Data& d_dirac, const bool useHRZ) const
801    {
802        escriptDataC _mat=mat.getDataC();
803        escriptDataC _D=D.getDataC();
804        escriptDataC _d=d.getDataC();
805        escriptDataC _d_dirac=d_dirac.getDataC();
806    
807        Finley_Mesh* mesh=m_finleyMesh.get();
808    
809  }      Finley_Assemble_LumpedSystem(mesh->Nodes, mesh->Elements, &_mat, &_D, useHRZ);
810        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  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        Finley_Mesh* mesh=m_finleyMesh.get();
827    
828        escriptDataC _rhs=rhs.getDataC();
829        escriptDataC _X=X.getDataC();
830        escriptDataC _Y=Y.getDataC();
831        escriptDataC _y=y.getDataC();
832        escriptDataC _y_contact=y_contact.getDataC();
833        escriptDataC _y_dirac=y_dirac.getDataC();
834    
835        Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y);
836        checkFinleyError();
837    
838     escriptDataC _rhs=rhs.getDataC();      Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y);
839     escriptDataC _X=X.getDataC();      checkFinleyError();
    escriptDataC _Y=Y.getDataC();  
    escriptDataC _y=y.getDataC();  
    escriptDataC _y_contact=y_contact.getDataC();  
840    
841     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );      Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y_contact);
842     checkFinleyError();      checkFinleyError();
843    
844     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );      Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs, 0, 0, 0, 0, 0, &_y_dirac);
845     checkFinleyError();      checkFinleyError();
   
    Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );  
    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                                             TransportProblemAdapter& 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          const escript::Data& d_contact, const escript::Data& y_contact,
857  {          const escript::Data& d_dirac, const escript::Data& y_dirac) const
858     DataTypes::ShapeType shape;  {
859     source.expand();      TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
860     escriptDataC _source=source.getDataC();      if (!tpa)
861     escriptDataC _M=M.getDataC();          throw FinleyAdapterException("finley only supports Paso transport problems.");
862     escriptDataC _A=A.getDataC();  
863     escriptDataC _B=B.getDataC();      source.expand();
864     escriptDataC _C=C.getDataC();      escriptDataC _source=source.getDataC();
865     escriptDataC _D=D.getDataC();      escriptDataC _M=M.getDataC();
866     escriptDataC _X=X.getDataC();      escriptDataC _A=A.getDataC();
867     escriptDataC _Y=Y.getDataC();      escriptDataC _B=B.getDataC();
868     escriptDataC _d=d.getDataC();      escriptDataC _C=C.getDataC();
869     escriptDataC _y=y.getDataC();      escriptDataC _D=D.getDataC();
870     escriptDataC _d_contact=d_contact.getDataC();      escriptDataC _X=X.getDataC();
871     escriptDataC _y_contact=y_contact.getDataC();      escriptDataC _Y=Y.getDataC();
872        escriptDataC _d=d.getDataC();
873     Finley_Mesh* mesh=m_finleyMesh.get();      escriptDataC _y=y.getDataC();
874     Paso_TransportProblem* _tp = tp.getPaso_TransportProblem();      escriptDataC _d_contact=d_contact.getDataC();
875        escriptDataC _y_contact=y_contact.getDataC();
876     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );      escriptDataC _d_dirac=d_dirac.getDataC();
877     checkFinleyError();      escriptDataC _y_dirac=y_dirac.getDataC();
878    
879     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );      Finley_Mesh* mesh=m_finleyMesh.get();
880     checkFinleyError();      Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
881    
882     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );      Finley_Assemble_PDE(mesh->Nodes, mesh->Elements, _tp->mass_matrix,
883     checkFinleyError();                          &_source, 0, 0, 0, &_M, 0, 0);
884        checkFinleyError();
885     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );  
886     checkFinleyError();      Finley_Assemble_PDE(mesh->Nodes, mesh->Elements, _tp->transport_matrix,
887  }                          &_source, &_A, &_B, &_C, &_D, &_X, &_Y);
888        checkFinleyError();
889  //  
890  // interpolates data between different function spaces:      Finley_Assemble_PDE(mesh->Nodes, mesh->FaceElements, _tp->transport_matrix,
891  //                          &_source, 0, 0, 0, &_d, 0, &_y);
892  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const      checkFinleyError();
893  {  
894     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));      Finley_Assemble_PDE(mesh->Nodes, mesh->ContactElements,
895     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));                          _tp->transport_matrix, &_source, 0, 0, 0, &_d_contact, 0, &_y_contact);
896     if (inDomain!=*this)        checkFinleyError();
897        throw FinleyAdapterException("Error - Illegal domain of interpolant.");  
898     if (targetDomain!=*this)      Finley_Assemble_PDE(mesh->Nodes, mesh->Points, _tp->transport_matrix,
899        throw FinleyAdapterException("Error - Illegal domain of interpolation target.");                          &_source, 0, 0, 0, &_d_dirac, 0, &_y_dirac);
900        checkFinleyError();
901     Finley_Mesh* mesh=m_finleyMesh.get();  }
902     escriptDataC _target=target.getDataC();  
903     escriptDataC _in=in.getDataC();  //
904     switch(in.getFunctionSpace().getTypeCode()) {  // interpolates data between different function spaces
905     case(Nodes):  //
906        switch(target.getFunctionSpace().getTypeCode()) {  void MeshAdapter::interpolateOnDomain(escript::Data& target, const escript::Data& in) const
907        case(Nodes):  {
908        case(ReducedNodes):      const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
909        case(DegreesOfFreedom):      const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
910        case(ReducedDegreesOfFreedom):      if (inDomain!=*this)
911        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);          throw FinleyAdapterException("Error - Illegal domain of interpolant.");
912        break;      if (targetDomain!=*this)
913        case(Elements):          throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
914        case(ReducedElements):  
915        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);      Finley_Mesh* mesh=m_finleyMesh.get();
916        break;      escriptDataC _target=target.getDataC();
917        case(FaceElements):      escriptDataC _in=in.getDataC();
918        case(ReducedFaceElements):      switch(in.getFunctionSpace().getTypeCode()) {
919        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);          case Nodes:
920        break;              switch(target.getFunctionSpace().getTypeCode()) {
921        case(Points):                  case Nodes:
922        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);                  case ReducedNodes:
923        break;                  case DegreesOfFreedom:
924        case(ContactElementsZero):                  case ReducedDegreesOfFreedom:
925        case(ReducedContactElementsZero):                      Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
926        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);                      break;
927        break;                  case Elements:
928        case(ContactElementsOne):                  case ReducedElements:
929        case(ReducedContactElementsOne):                      Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
930        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);                      break;
931        break;                  case FaceElements:
932        default:                  case ReducedFaceElements:
933           stringstream temp;                      Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
934           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();                      break;
935           throw FinleyAdapterException(temp.str());                  case Points:
936           break;                      Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
937        }                      break;
938        break;                  case ContactElementsZero:
939     case(ReducedNodes):                  case ReducedContactElementsZero:
940        switch(target.getFunctionSpace().getTypeCode()) {                  case ContactElementsOne:
941        case(Nodes):                  case ReducedContactElementsOne:
942        case(ReducedNodes):                      Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
943        case(DegreesOfFreedom):                      break;
944        case(ReducedDegreesOfFreedom):                  default:
945        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);                      stringstream temp;
946        break;                      temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
947        case(Elements):                      throw FinleyAdapterException(temp.str());
948        case(ReducedElements):                      break;
949        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);              }
950        break;              break;
951        case(FaceElements):          case ReducedNodes:
952        case(ReducedFaceElements):              switch(target.getFunctionSpace().getTypeCode()) {
953        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);                  case Nodes:
954        break;                  case ReducedNodes:
955        case(Points):                  case DegreesOfFreedom:
956        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);                  case ReducedDegreesOfFreedom:
957        break;                      Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
958        case(ContactElementsZero):                      break;
959        case(ReducedContactElementsZero):                  case Elements:
960        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);                  case ReducedElements:
961        break;                      Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
962        case(ContactElementsOne):                      break;
963        case(ReducedContactElementsOne):                  case FaceElements:
964        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);                  case ReducedFaceElements:
965        break;                      Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
966        default:                      break;
967           stringstream temp;                  case Points:
968           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();                      Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
969           throw FinleyAdapterException(temp.str());                      break;
970           break;                  case ContactElementsZero:
971        }                  case ReducedContactElementsZero:
972        break;                  case ContactElementsOne:
973     case(Elements):                  case ReducedContactElementsOne:
974        if (target.getFunctionSpace().getTypeCode()==Elements) {                      Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
975           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);                      break;
976        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {                  default:
977           Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);                      stringstream temp;
978        } else {                      temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
979           throw FinleyAdapterException("Error - No interpolation with data on elements possible.");                      throw FinleyAdapterException(temp.str());
980        }                      break;
981        break;              }
982     case(ReducedElements):              break;
983        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {          case Elements:
984           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);              if (target.getFunctionSpace().getTypeCode()==Elements) {
985        } else {                  Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
986           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");              } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
987        }                  Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
988        break;              } else {
989     case(FaceElements):                  throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
990        if (target.getFunctionSpace().getTypeCode()==FaceElements) {              }
991           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);              break;
992        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {          case ReducedElements:
993           Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);              if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
994        } else {                  Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
995           throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");              } else {
996        }                  throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
997        break;              }
998     case(ReducedFaceElements):              break;
999        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {          case FaceElements:
1000           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);              if (target.getFunctionSpace().getTypeCode()==FaceElements) {
1001        } else {                  Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1002           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");              } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1003        }                  Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1004        break;              } else {
1005     case(Points):                  throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1006        if (target.getFunctionSpace().getTypeCode()==Points) {              }
1007           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);              break;
1008        } else {          case ReducedFaceElements:
1009           throw FinleyAdapterException("Error - No interpolation with data on points possible.");              if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1010        }                  Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1011        break;              } else {
1012     case(ContactElementsZero):                  throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
1013     case(ContactElementsOne):              }
1014        if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {              break;
1015           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);          case Points:
1016        } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {              if (target.getFunctionSpace().getTypeCode()==Points) {
1017           Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);                  Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1018        } else {              } else {
1019           throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");                  throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1020        }              }
1021        break;              break;
1022     case(ReducedContactElementsZero):          case ContactElementsZero:
1023     case(ReducedContactElementsOne):          case ContactElementsOne:
1024        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {              if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1025           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);                  Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1026        } else {              } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1027           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");                  Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1028        }              } else {
1029        break;                  throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1030     case(DegreesOfFreedom):                    }
1031        switch(target.getFunctionSpace().getTypeCode()) {              break;
1032        case(ReducedDegreesOfFreedom):          case ReducedContactElementsZero:
1033        case(DegreesOfFreedom):          case ReducedContactElementsOne:
1034        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);              if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1035        break;                  Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1036                  } else {
1037        case(Nodes):                  throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1038        case(ReducedNodes):              }
1039        if (getMPISize()>1) {              break;
1040           escript::Data temp=escript::Data(in);          case DegreesOfFreedom:
1041           temp.expand();              switch(target.getFunctionSpace().getTypeCode()) {
1042           escriptDataC _in2 = temp.getDataC();                  case ReducedDegreesOfFreedom:
1043           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);                  case DegreesOfFreedom:
1044        } else {                      Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1045           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);                      break;
1046        }  
1047        break;                  case Nodes:
1048        case(Elements):                  case ReducedNodes:
1049        case(ReducedElements):                      if (getMPISize()>1) {
1050        if (getMPISize()>1) {                          escript::Data temp=escript::Data(in);
1051           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );                          temp.expand();
1052           escriptDataC _in2 = temp.getDataC();                          escriptDataC _in2 = temp.getDataC();
1053           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);                          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1054        } else {                      } else {
1055           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);                          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1056        }                      }
1057        break;                      break;
1058        case(FaceElements):                  case Elements:
1059        case(ReducedFaceElements):                  case ReducedElements:
1060        if (getMPISize()>1) {                      if (getMPISize()>1) {
1061           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );                          escript::Data temp=escript::Data(in, continuousFunction(*this));
1062           escriptDataC _in2 = temp.getDataC();                          escriptDataC _in2 = temp.getDataC();
1063           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);                          Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1064                          } else {
1065        } else {                          Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1066           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);                      }
1067        }                      break;
1068        break;                  case FaceElements:
1069        case(Points):                  case ReducedFaceElements:
1070        if (getMPISize()>1) {                      if (getMPISize()>1) {
1071           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );                          escript::Data temp=escript::Data(in, continuousFunction(*this));
1072           escriptDataC _in2 = temp.getDataC();                          escriptDataC _in2 = temp.getDataC();
1073        } else {                          Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1074           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);                      } else {
1075        }                          Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1076        break;                      }
1077        case(ContactElementsZero):                      break;
1078        case(ContactElementsOne):                  case Points:
1079        case(ReducedContactElementsZero):                      if (getMPISize()>1) {
1080        case(ReducedContactElementsOne):                          //escript::Data temp=escript::Data(in, continuousFunction(*this) );
1081        if (getMPISize()>1) {                          //escriptDataC _in2 = temp.getDataC();
1082           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );                      } else {
1083           escriptDataC _in2 = temp.getDataC();                          Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1084           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);                      }
1085        } else {                      break;
1086           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);                  case ContactElementsZero:
1087        }                  case ContactElementsOne:
1088        break;                  case ReducedContactElementsZero:
1089        default:                  case ReducedContactElementsOne:
1090           stringstream temp;                      if (getMPISize()>1) {
1091           temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();                          escript::Data temp=escript::Data(in, continuousFunction(*this));
1092           throw FinleyAdapterException(temp.str());                          escriptDataC _in2 = temp.getDataC();
1093           break;                          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1094        }                      } else {
1095        break;                          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1096     case(ReducedDegreesOfFreedom):                      }
1097        switch(target.getFunctionSpace().getTypeCode()) {                      break;
1098        case(Nodes):                  default:
1099        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");                      stringstream temp;
1100        break;                      temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1101        case(ReducedNodes):                      throw FinleyAdapterException(temp.str());
1102        if (getMPISize()>1) {                      break;
1103           escript::Data temp=escript::Data(in);              }
1104           temp.expand();              break;
1105           escriptDataC _in2 = temp.getDataC();          case ReducedDegreesOfFreedom:
1106           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);              switch(target.getFunctionSpace().getTypeCode()) {
1107        } else {                  case Nodes:
1108           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);                      throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1109        }                      break;
1110        break;                  case ReducedNodes:
1111        case(DegreesOfFreedom):                      if (getMPISize()>1) {
1112        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");                          escript::Data temp=escript::Data(in);
1113        break;                          temp.expand();
1114        case(ReducedDegreesOfFreedom):                          escriptDataC _in2 = temp.getDataC();
1115        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);                          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1116        break;                      } else {
1117        case(Elements):                          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1118        case(ReducedElements):                      }
1119        if (getMPISize()>1) {                      break;
1120           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );                  case DegreesOfFreedom:
1121           escriptDataC _in2 = temp.getDataC();                      throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1122           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);                      break;
1123        } else {                  case ReducedDegreesOfFreedom:
1124           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);                      Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1125        }                      break;
1126        break;                  case Elements:
1127        case(FaceElements):                  case ReducedElements:
1128        case(ReducedFaceElements):                      if (getMPISize()>1) {
1129        if (getMPISize()>1) {                          escript::Data temp=escript::Data(in, reducedContinuousFunction(*this) );
1130           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );                          escriptDataC _in2 = temp.getDataC();
1131           escriptDataC _in2 = temp.getDataC();                          Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1132           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);                      } else {
1133        } else {                          Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1134           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);                      }
1135        }                      break;
1136        break;                  case FaceElements:
1137        case(Points):                  case ReducedFaceElements:
1138        if (getMPISize()>1) {                      if (getMPISize()>1) {
1139           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );                          escript::Data temp=escript::Data(in, reducedContinuousFunction(*this) );
1140           escriptDataC _in2 = temp.getDataC();                          escriptDataC _in2 = temp.getDataC();
1141           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);                          Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1142        } else {                      } else {
1143           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);                          Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1144        }                      }
1145        break;                      break;
1146        case(ContactElementsZero):                  case Points:
1147        case(ContactElementsOne):                      if (getMPISize()>1) {
1148        case(ReducedContactElementsZero):                          escript::Data temp=escript::Data(in, reducedContinuousFunction(*this));
1149        case(ReducedContactElementsOne):                          escriptDataC _in2 = temp.getDataC();
1150        if (getMPISize()>1) {                          Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1151           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );                      } else {
1152           escriptDataC _in2 = temp.getDataC();                          Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1153           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);                      }
1154        } else {                      break;
1155           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);                  case ContactElementsZero:
1156        }                  case ContactElementsOne:
1157        break;                  case ReducedContactElementsZero:
1158        default:                  case ReducedContactElementsOne:
1159           stringstream temp;                      if (getMPISize()>1) {
1160           temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();                          escript::Data temp=escript::Data(in, reducedContinuousFunction(*this));
1161           throw FinleyAdapterException(temp.str());                          escriptDataC _in2 = temp.getDataC();
1162           break;                          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1163        }                      } else {
1164        break;                          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1165     default:                      }
1166        stringstream temp;                      break;
1167        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();                  default:
1168        throw FinleyAdapterException(temp.str());                      stringstream temp;
1169        break;                      temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1170     }                      throw FinleyAdapterException(temp.str());
1171     checkFinleyError();                      break;
1172                }
1173                break;
1174            default:
1175                stringstream temp;
1176                temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1177                throw FinleyAdapterException(temp.str());
1178                break;
1179        }
1180        checkFinleyError();
1181  }  }
1182    
1183  //  //
1184  // copies the locations of sample points into x:  // copies the locations of sample points into x
1185  //  //
1186  void MeshAdapter::setToX(escript::Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
1187  {  {
1188     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));      const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1189     if (argDomain!=*this)      if (argDomain!=*this)
1190        throw FinleyAdapterException("Error - Illegal domain of data point locations");          throw FinleyAdapterException("Error - Illegal domain of data point locations");
1191     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
1192     // in case of values node coordinates we can do the job directly:      // in case of appropriate function space we can do the job directly:
1193     if (arg.getFunctionSpace().getTypeCode()==Nodes) {      if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1194        escriptDataC _arg=arg.getDataC();          escriptDataC _arg=arg.getDataC();
1195        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);          Finley_Assemble_NodeCoordinates(mesh->Nodes, &_arg);
1196     } else {      } else {
1197        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);          escript::Data tmp_data=Vector(0., continuousFunction(*this), true);
1198        escriptDataC _tmp_data=tmp_data.getDataC();          escriptDataC _tmp_data=tmp_data.getDataC();
1199        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);          Finley_Assemble_NodeCoordinates(mesh->Nodes, &_tmp_data);
1200        // this is then interpolated onto arg:          // this is then interpolated onto arg:
1201        interpolateOnDomain(arg,tmp_data);          interpolateOnDomain(arg, tmp_data);
1202     }      }
1203     checkFinleyError();      checkFinleyError();
1204  }  }
1205    
1206  //  //
1207  // 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
1208  //  //
1209  void MeshAdapter::setToNormal(escript::Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
1210  {  {
1211  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/      const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1212     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));      if (normalDomain!=*this)
1213     if (normalDomain!=*this)          throw FinleyAdapterException("Error - Illegal domain of normal locations");
1214        throw FinleyAdapterException("Error - Illegal domain of normal locations");      Finley_Mesh* mesh=m_finleyMesh.get();
1215     Finley_Mesh* mesh=m_finleyMesh.get();      escriptDataC _normal=normal.getDataC();
1216     escriptDataC _normal=normal.getDataC();      switch(normal.getFunctionSpace().getTypeCode()) {
1217     switch(normal.getFunctionSpace().getTypeCode()) {          case Nodes:
1218     case(Nodes):              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1219     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");              break;
1220     break;          case ReducedNodes:
1221     case(ReducedNodes):              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1222     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");              break;
1223     break;          case Elements:
1224     case(Elements):              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1225     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");              break;
1226     break;          case ReducedElements:
1227     case(ReducedElements):              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1228     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");              break;
1229     break;          case FaceElements:
1230     case (FaceElements):          case ReducedFaceElements:
1231     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);              Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1232     break;              break;
1233     case (ReducedFaceElements):          case Points:
1234     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1235     break;              break;
1236     case(Points):          case ContactElementsOne:
1237     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");          case ContactElementsZero:
1238     break;          case ReducedContactElementsOne:
1239     case (ContactElementsOne):          case ReducedContactElementsZero:
1240     case (ContactElementsZero):              Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1241     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);              break;
1242     break;          case DegreesOfFreedom:
1243     case (ReducedContactElementsOne):              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1244     case (ReducedContactElementsZero):              break;
1245     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);          case ReducedDegreesOfFreedom:
1246     break;              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1247     case(DegreesOfFreedom):              break;
1248     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");          default:
1249     break;              stringstream temp;
1250     case(ReducedDegreesOfFreedom):              temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1251     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");              throw FinleyAdapterException(temp.str());
1252     break;              break;
1253     default:      }
1254        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();  
1255  }  }
1256    
1257  //  //
1258  // interpolates data to other domain:  // interpolates data to other domain
1259  //  //
1260  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1261  {  {
1262     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();      escript::const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1263     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());      const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1264     if (targetDomain!=this)      if (targetDomain!=this)
1265        throw FinleyAdapterException("Error - Illegal domain of interpolation target");          throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1266    
1267     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");      throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
1268  }  }
1269    
1270  //  //
1271  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined on arg
1272  //  //
1273  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
1274  {  {
1275     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));      const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1276     if (argDomain!=*this)      if (argDomain!=*this)
1277        throw FinleyAdapterException("Error - Illegal domain of integration kernel");          throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1278    
1279     double blocktimer_start = blocktimer_time();      double blocktimer_start = blocktimer_time();
1280     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
1281     escriptDataC _temp;      escriptDataC _temp;
1282     escript::Data temp;      escript::Data temp;
1283     escriptDataC _arg=arg.getDataC();      escriptDataC _arg=arg.getDataC();
1284     switch(arg.getFunctionSpace().getTypeCode()) {      switch(arg.getFunctionSpace().getTypeCode()) {
1285     case(Nodes):          case Nodes:
1286     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );              temp=escript::Data(arg, escript::function(*this));
1287     _temp=temp.getDataC();              _temp=temp.getDataC();
1288     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);              Finley_Assemble_integrate(mesh->Nodes, mesh->Elements, &_temp, &integrals[0]);
1289     break;              break;
1290     case(ReducedNodes):          case ReducedNodes:
1291     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );              temp=escript::Data( arg, escript::function(*this) );
1292     _temp=temp.getDataC();              _temp=temp.getDataC();
1293     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);              Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1294     break;              break;
1295     case(Elements):          case Elements:
1296     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);          case ReducedElements:
1297     break;              Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1298     case(ReducedElements):              break;
1299     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);          case FaceElements:
1300     break;          case ReducedFaceElements:
1301     case(FaceElements):              Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1302     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);              break;
1303     break;          case Points:
1304     case(ReducedFaceElements):              throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1305     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);              break;
1306     break;          case ContactElementsZero:
1307     case(Points):          case ReducedContactElementsZero:
1308     throw FinleyAdapterException("Error - Integral of data on points is not supported.");          case ContactElementsOne:
1309     break;          case ReducedContactElementsOne:
1310     case(ContactElementsZero):              Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1311     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);              break;
1312     break;          case DegreesOfFreedom:
1313     case(ReducedContactElementsZero):          case ReducedDegreesOfFreedom:
1314     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);              temp=escript::Data(arg, escript::function(*this));
1315     break;              _temp=temp.getDataC();
1316     case(ContactElementsOne):              Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1317     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);              break;
1318     break;          default:
1319     case(ReducedContactElementsOne):              stringstream temp;
1320     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();
1321     break;              throw FinleyAdapterException(temp.str());
1322     case(DegreesOfFreedom):              break;
1323     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );      }
1324     _temp=temp.getDataC();      checkFinleyError();
1325     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);      blocktimer_increment("integrate()", blocktimer_start);
1326     break;  }
1327     case(ReducedDegreesOfFreedom):  
1328     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );  //
1329     _temp=temp.getDataC();  // calculates the gradient of arg
1330     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);  //
1331     break;  void MeshAdapter::setToGradient(escript::Data& grad, const escript::Data& arg) const
1332     default:  {
1333        stringstream temp;      const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1334        temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();      if (argDomain!=*this)
1335        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1336        break;      const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1337     }      if (gradDomain!=*this)
1338     checkFinleyError();          throw FinleyAdapterException("Error - Illegal domain of gradient");
1339     blocktimer_increment("integrate()", blocktimer_start);  
1340  }      Finley_Mesh* mesh=m_finleyMesh.get();
1341        escriptDataC _grad=grad.getDataC();
1342  //      escriptDataC nodeDataC;
1343  // calculates the gradient of arg:      escript::Data temp;
1344  //      if (getMPISize()>1) {
1345  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const          if (arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom) {
1346  {              temp=escript::Data(arg, continuousFunction(*this));
1347     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));              nodeDataC = temp.getDataC();
1348     if (argDomain!=*this)          } else if(arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom) {
1349        throw FinleyAdapterException("Error - Illegal domain of gradient argument");              temp=escript::Data(arg, reducedContinuousFunction(*this));
1350     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));              nodeDataC = temp.getDataC();
1351     if (gradDomain!=*this)          } else {
1352        throw FinleyAdapterException("Error - Illegal domain of gradient");              nodeDataC = arg.getDataC();
1353            }
1354     Finley_Mesh* mesh=m_finleyMesh.get();      } else {
1355     escriptDataC _grad=grad.getDataC();          nodeDataC = arg.getDataC();
1356     escriptDataC nodeDataC;      }
1357     escript::Data temp;      switch(grad.getFunctionSpace().getTypeCode()) {
1358     if (getMPISize()>1) {          case Nodes:
1359        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {              throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1360           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );              break;
1361           nodeDataC = temp.getDataC();          case ReducedNodes:
1362        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {              throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1363           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );              break;
1364           nodeDataC = temp.getDataC();          case Elements:
1365        } else {          case ReducedElements:
1366           nodeDataC = arg.getDataC();              Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1367        }              break;
1368     } else {          case FaceElements:
1369        nodeDataC = arg.getDataC();          case ReducedFaceElements:
1370     }              Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1371     switch(grad.getFunctionSpace().getTypeCode()) {              break;
1372     case(Nodes):          case Points:
1373     throw FinleyAdapterException("Error - Gradient at nodes is not supported.");              throw FinleyAdapterException("Error - Gradient at points is not supported.");
1374     break;              break;
1375     case(ReducedNodes):          case ContactElementsZero:
1376     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");          case ReducedContactElementsZero:
1377     break;          case ContactElementsOne:
1378     case(Elements):          case ReducedContactElementsOne:
1379     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);              Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1380     break;              break;
1381     case(ReducedElements):          case DegreesOfFreedom:
1382     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);              throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1383     break;              break;
1384     case(FaceElements):          case ReducedDegreesOfFreedom:
1385     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);              throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1386     break;              break;
1387     case(ReducedFaceElements):          default:
1388     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);              stringstream temp;
1389     break;              temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1390     case(Points):              throw FinleyAdapterException(temp.str());
1391     throw FinleyAdapterException("Error - Gradient at points is not supported.");              break;
1392     break;      }
1393     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();  
1394  }  }
1395    
1396  //  //
1397  // returns the size of elements:  // returns the size of elements
1398  //  //
1399  void MeshAdapter::setToSize(escript::Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
1400  {  {
1401     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
1402     escriptDataC tmp=size.getDataC();      escriptDataC tmp=size.getDataC();
1403     switch(size.getFunctionSpace().getTypeCode()) {      switch(size.getFunctionSpace().getTypeCode()) {
1404     case(Nodes):          case Nodes:
1405     throw FinleyAdapterException("Error - Size of nodes is not supported.");              throw FinleyAdapterException("Error - Size of nodes is not supported.");
1406     break;              break;
1407     case(ReducedNodes):          case ReducedNodes:
1408     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");              throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1409     break;              break;
1410     case(Elements):          case Elements:
1411     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);          case ReducedElements:
1412     break;              Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1413     case(ReducedElements):              break;
1414     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);          case FaceElements:
1415     break;          case ReducedFaceElements:
1416     case(FaceElements):              Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1417     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);              break;
1418     break;          case Points:
1419     case(ReducedFaceElements):              throw FinleyAdapterException("Error - Size of point elements is not supported.");
1420     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);              break;
1421     break;          case ContactElementsZero:
1422     case(Points):          case ContactElementsOne:
1423     throw FinleyAdapterException("Error - Size of point elements is not supported.");          case ReducedContactElementsZero:
1424     break;          case ReducedContactElementsOne:
1425     case(ContactElementsZero):              Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1426     case(ContactElementsOne):              break;
1427     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);          case DegreesOfFreedom:
1428     break;              throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1429     case(ReducedContactElementsZero):              break;
1430     case(ReducedContactElementsOne):          case ReducedDegreesOfFreedom:
1431     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);              throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1432     break;              break;
1433     case(DegreesOfFreedom):          default:
1434     throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");              stringstream temp;
1435     break;              temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1436     case(ReducedDegreesOfFreedom):              throw FinleyAdapterException(temp.str());
1437     throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");              break;
1438     break;      }
1439     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();  
1440  }  }
1441    
1442  //  //
# Line 1466  void MeshAdapter::setToSize(escript::Dat Line 1444  void MeshAdapter::setToSize(escript::Dat
1444  //  //
1445  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1446  {  {
1447     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
1448     escriptDataC tmp;      const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1449     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));      if (newDomain!=*this)
1450     if (newDomain!=*this)          throw FinleyAdapterException("Error - Illegal domain of new point locations");
1451        throw FinleyAdapterException("Error - Illegal domain of new point locations");      if (new_x.getFunctionSpace() == continuousFunction(*this)) {
1452     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {          Finley_Mesh_setCoordinates(mesh, new_x);
1453         tmp = new_x.getDataC();      } else {
1454         Finley_Mesh_setCoordinates(mesh,&tmp);          throw FinleyAdapterException("As of escript version 3.3 SetX() only accepts ContinuousFunction arguments. Please interpolate.");
1455     } else {          //escript::Data new_x_inter=escript::Data(new_x, continuousFunction(*this));
1456         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );          //Finley_Mesh_setCoordinates(mesh, new_x_inter);
1457         tmp = new_x_inter.getDataC();      }
1458         Finley_Mesh_setCoordinates(mesh,&tmp);      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);  
   
    return;  
 }  
   
 //  
 // saves mesh and optionally data arrays in VTK format  
 //  
 void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const  
 {  
    int num_data;  
    char **names;  
    escriptDataC *data;  
    escriptDataC **ptr_data;  
   
    extractArgsFromDict(arg, num_data, names, data, ptr_data);  
    Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());  
    checkFinleyError();  
   
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0; i<num_data; i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
1459  }  }
1460    
1461  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1462  {  {
1463  #ifdef PASO_MPI      if (getMPISize() > 1) {
1464      index_t myFirstNode=0, myLastNode=0, k=0;  #ifdef ESYS_MPI
1465      index_t* globalNodeIndex=0;          index_t myFirstNode=0, myLastNode=0, k=0;
1466      Finley_Mesh* mesh_p=m_finleyMesh.get();          index_t* globalNodeIndex=0;
1467      if (fs_code == FINLEY_REDUCED_NODES)          Finley_Mesh* mesh_p=m_finleyMesh.get();
1468      {          /*
1469      myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);           * this method is only used by saveDataCSV which would use the returned
1470      myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);           * values for reduced nodes wrongly so this case is disabled for now
1471      globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);          if (fs_code == FINLEY_REDUCED_NODES) {
1472      }              myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1473      else              myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1474      {              globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1475      myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);          } else
1476      myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);          */
1477      globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);          if (fs_code == FINLEY_NODES) {
1478      }              myFirstNode = mesh_p->Nodes->getFirstNode();
1479      k=globalNodeIndex[id];              myLastNode = mesh_p->Nodes->getLastNode();
1480      return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );              globalNodeIndex = mesh_p->Nodes->borrowGlobalNodesIndex();
1481            } else {
1482                throw FinleyAdapterException("Unsupported function space type for ownSample()");
1483            }
1484    
1485            k=globalNodeIndex[id];
1486            return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1487  #endif  #endif
1488        }
1489      return true;      return true;
1490  }  }
1491    
1492    
   
1493  //  //
1494  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1495  //  //
1496  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  escript::ASM_ptr MeshAdapter::newSystemMatrix(const int row_blocksize,
1497                                                   const int row_blocksize,                              const escript::FunctionSpace& row_functionspace,
1498                                                   const escript::FunctionSpace& row_functionspace,                              const int column_blocksize,
1499                                                   const int column_blocksize,                              const escript::FunctionSpace& column_functionspace,
1500                                                   const escript::FunctionSpace& column_functionspace,                              const int type) const
1501                                                   const int type) const  {
1502  {      // is the domain right?
1503     int reduceRowOrder=0;      const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1504     int reduceColOrder=0;      if (row_domain!=*this)
1505     // is the domain right?          throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1506     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));      const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1507     if (row_domain!=*this)      if (col_domain!=*this)
1508        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.");
1509     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));  
1510     if (col_domain!=*this)      int reduceRowOrder=0;
1511        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");      int reduceColOrder=0;
1512     // is the function space type right      // is the function space type right?
1513     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {      if (row_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1514        reduceRowOrder=0;          reduceRowOrder=1;
1515     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      } else if (row_functionspace.getTypeCode() != DegreesOfFreedom) {
1516        reduceRowOrder=1;          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1517     } else {      }
1518        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");      if (column_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1519     }          reduceColOrder=1;
1520     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {      } else if (column_functionspace.getTypeCode() != DegreesOfFreedom) {
1521        reduceColOrder=0;          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1522     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      }
1523        reduceColOrder=1;  
1524     } else {      // generate matrix:
1525        throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(
1526     }              getFinley_Mesh(), reduceRowOrder, reduceColOrder);
1527     // generate matrix:      checkFinleyError();
1528        Paso_SystemMatrix* fsystemMatrix;
1529     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);      const int trilinos = 0;
1530     checkFinleyError();      if (trilinos) {
    Paso_SystemMatrix* fsystemMatrix;  
    int trilinos = 0;  
    if (trilinos) {  
1531  #ifdef TRILINOS  #ifdef TRILINOS
1532        /* Allocation Epetra_VrbMatrix here */          // FIXME: Allocation Epetra_VrbMatrix here...
1533  #endif  #endif
1534     }      } else {
1535     else {          fsystemMatrix=Paso_SystemMatrix_alloc(type, fsystemMatrixPattern,
1536        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);                  row_blocksize, column_blocksize, FALSE);
1537     }      }
1538     checkPasoError();      checkPasoError();
1539     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1540     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);      SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1541        return escript::ASM_ptr(sma);
1542  }  }
1543    
1544  //  //
1545  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1546  //  //
1547  TransportProblemAdapter MeshAdapter::newTransportProblem(  escript::ATP_ptr MeshAdapter::newTransportProblem(const int blocksize,
1548                                                           const bool useBackwardEuler,          const escript::FunctionSpace& functionspace, const int type) const
1549                                                           const int blocksize,  {
1550                                                           const escript::FunctionSpace& functionspace,      // is the domain right?
1551                                                           const int type) const      const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1552  {      if (domain!=*this)
1553     int reduceOrder=0;          throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1554     // is the domain right?  
1555     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));      // is the function space type right?
1556     if (domain!=*this)      int reduceOrder=0;
1557        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");      if (functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1558     // is the function space type right          reduceOrder=1;
1559     if (functionspace.getTypeCode()==DegreesOfFreedom) {      } else if (functionspace.getTypeCode() != DegreesOfFreedom) {
1560        reduceOrder=0;          throw FinleyAdapterException("Error - illegal function space type for transport problem.");
1561     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      }
1562        reduceOrder=1;  
1563     } else {      // generate transport problem:
1564        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(
1565     }              getFinley_Mesh(), reduceOrder, reduceOrder);
1566     // generate matrix:      checkFinleyError();
1567        Paso_TransportProblem* transportProblem;
1568     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);      transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern, blocksize);
1569     checkFinleyError();      checkPasoError();
1570     Paso_TransportProblem* transportProblem;      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1571     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);      TransportProblemAdapter* tpa=new TransportProblemAdapter(
1572     checkPasoError();              transportProblem, blocksize, functionspace);
1573     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);      return escript::ATP_ptr(tpa);
    return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);  
1574  }  }
1575    
1576  //  //
1577  // vtkObject MeshAdapter::createVtkObject() const  // returns true if data on functionSpaceCode is considered as being cell centered
 // TODO:  
1578  //  //
   
 //  
 // returns true if data at the atom_type is considered as being cell centered:  
1579  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1580  {  {
1581     switch(functionSpaceCode) {      switch(functionSpaceCode) {
1582     case(Nodes):          case Nodes:
1583     case(DegreesOfFreedom):          case DegreesOfFreedom:
1584     case(ReducedDegreesOfFreedom):          case ReducedDegreesOfFreedom:
1585     return false;              return false;
1586     break;          case Elements:
1587     case(Elements):          case FaceElements:
1588     case(FaceElements):          case Points:
1589     case(Points):          case ContactElementsZero:
1590     case(ContactElementsZero):          case ContactElementsOne:
1591     case(ContactElementsOne):          case ReducedElements:
1592     case(ReducedElements):          case ReducedFaceElements:
1593     case(ReducedFaceElements):          case ReducedContactElementsZero:
1594     case(ReducedContactElementsZero):          case ReducedContactElementsOne:
1595     case(ReducedContactElementsOne):              return true;
1596     return true;          default:
1597     break;              stringstream temp;
1598     default:              temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1599        stringstream temp;              throw FinleyAdapterException(temp.str());
1600        temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;              break;
1601        throw FinleyAdapterException(temp.str());      }
1602        break;      return false;
    }  
    checkFinleyError();  
    return false;  
1603  }  }
1604    
1605  bool  bool
1606  MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const  MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1607  {  {
1608     /* 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
1609      class 1: DOF <-> Nodes         back and forth]:
1610      class 2: ReducedDOF <-> ReducedNodes          class 1: DOF <-> Nodes
1611      class 3: Points          class 2: ReducedDOF <-> ReducedNodes
1612      class 4: Elements          class 3: Points
1613      class 5: ReducedElements          class 4: Elements
1614      class 6: FaceElements          class 5: ReducedElements
1615      class 7: ReducedFaceElements          class 6: FaceElements
1616      class 8: ContactElementZero <-> ContactElementOne          class 7: ReducedFaceElements
1617      class 9: ReducedContactElementZero <-> ReducedContactElementOne          class 8: ContactElementZero <-> ContactElementOne
1618            class 9: ReducedContactElementZero <-> ReducedContactElementOne
1619     There is also a set of lines. Interpolation is possible down a line but not between lines.  
1620     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
1621      line 0: class 3      not between lines.
1622      line 1: class 4,5      class 1 and 2 belong to all lines so aren't considered.
1623      line 2: class 6,7          line 0: class 3
1624      line 3: class 8,9          line 1: class 4,5
1625            line 2: class 6,7
1626     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
1627     eg hasnodes is true if we have at least one instance of Nodes.  
1628     */      For classes with multiple members (eg class 2) we have vars to record if
1629        there is at least one instance.
1630        e.g. hasnodes is true if we have at least one instance of Nodes.
1631        */
1632      if (fs.empty())      if (fs.empty())
     {  
1633          return false;          return false;
1634      }  
1635      vector<int> hasclass(10);      vector<int> hasclass(10);
1636      vector<int> hasline(4);      vector<int> hasline(4);    
1637      bool hasnodes=false;      bool hasnodes=false;
1638      bool hasrednodes=false;      bool hasrednodes=false;
1639      bool hascez=false;      bool hascez=false;
1640      bool hasrcez=false;      bool hasrcez=false;
1641      for (int i=0;i<fs.size();++i)      for (int i=0;i<fs.size();++i) {
1642      {          switch(fs[i]) {
1643      switch(fs[i])              case Nodes:
1644      {                  hasnodes=true; // fall through
1645      case(Nodes):   hasnodes=true;   // no break is deliberate              case DegreesOfFreedom:
1646      case(DegreesOfFreedom):                  hasclass[1]=1;
1647          hasclass[1]=1;                  break;
1648          break;              case ReducedNodes:
1649      case(ReducedNodes):    hasrednodes=true;    // no break is deliberate                  hasrednodes=true; // fall through
1650      case(ReducedDegreesOfFreedom):              case ReducedDegreesOfFreedom:
1651          hasclass[2]=1;                  hasclass[2]=1;
1652          break;                  break;
1653      case(Points):              case Points:
1654          hasline[0]=1;                  hasline[0]=1;
1655          hasclass[3]=1;                  hasclass[3]=1;
1656          break;                  break;
1657      case(Elements):              case Elements:
1658          hasclass[4]=1;                  hasclass[4]=1;
1659          hasline[1]=1;                  hasline[1]=1;
1660          break;                  break;
1661      case(ReducedElements):              case ReducedElements:
1662          hasclass[5]=1;                  hasclass[5]=1;
1663          hasline[1]=1;                  hasline[1]=1;
1664          break;                  break;
1665      case(FaceElements):              case FaceElements:
1666          hasclass[6]=1;                  hasclass[6]=1;
1667          hasline[2]=1;                  hasline[2]=1;
1668          break;                  break;
1669      case(ReducedFaceElements):              case ReducedFaceElements:
1670          hasclass[7]=1;                  hasclass[7]=1;
1671          hasline[2]=1;                  hasline[2]=1;
1672          break;                  break;
1673      case(ContactElementsZero):  hascez=true;    // no break is deliberate              case ContactElementsZero:
1674      case(ContactElementsOne):                  hascez=true; // fall through
1675          hasclass[8]=1;              case ContactElementsOne:
1676          hasline[3]=1;                  hasclass[8]=1;
1677          break;                  hasline[3]=1;
1678      case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate                  break;
1679      case(ReducedContactElementsOne):              case ReducedContactElementsZero:
1680          hasclass[9]=1;                  hasrcez=true; // fall through
1681          hasline[3]=1;              case ReducedContactElementsOne:
1682          break;                  hasclass[9]=1;
1683      default:                  hasline[3]=1;
1684          return false;                  break;
1685      }              default:
1686                    return false;
1687            }
1688      }      }
1689      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  
1690    
1691        // fail if we have more than one leaf group
1692      if (totlines>1)      if (totlines>1)
1693      {          return false; // there are at least two branches we can't interpolate between
1694      return false;   // there are at least two branches we can't interpolate between      else if (totlines==1) {
1695      }          if (hasline[0]==1)              // we have points
1696      else if (totlines==1)              resultcode=Points;
1697      {          else if (hasline[1]==1) {
1698      if (hasline[0]==1)      // we have points              if (hasclass[5]==1)
1699      {                  resultcode=ReducedElements;
1700          resultcode=Points;              else
1701      }                  resultcode=Elements;
1702      else if (hasline[1]==1)          } else if (hasline[2]==1) {
1703      {              if (hasclass[7]==1)
1704          if (hasclass[5]==1)                  resultcode=ReducedFaceElements;
1705          {              else
1706          resultcode=ReducedElements;                  resultcode=FaceElements;
1707          }          } else {   // so we must be in line3
1708          else              if (hasclass[9]==1) {
1709          {                  // need something from class 9
1710          resultcode=Elements;                  resultcode=(hasrcez ? ReducedContactElementsZero : ReducedContactElementsOne);
1711          }              } else {
1712      }                  // something from class 8
1713      else if (hasline[2]==1)                  resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1714      {              }
1715          if (hasclass[7]==1)          }
1716          {      } else { // totlines==0
1717          resultcode=ReducedFaceElements;          if (hasclass[2]==1) {
1718          }              // something from class 2
1719          else              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
1720          {          } else {
1721          resultcode=FaceElements;              // something from class 1
1722          }              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
1723      }          }
     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);  
     }  
1724      }      }
1725      return true;      return true;
1726  }  }
1727    
1728  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,
1729                                                 int functionSpaceType_target) const
1730  {  {
1731     switch(functionSpaceType_source) {      switch(functionSpaceType_source) {
1732     case(Nodes):          case Nodes:
1733      switch(functionSpaceType_target) {              switch(functionSpaceType_target) {
1734      case(Nodes):                  case Nodes:
1735      case(ReducedNodes):                  case ReducedNodes:
1736      case(ReducedDegreesOfFreedom):                  case ReducedDegreesOfFreedom:
1737      case(DegreesOfFreedom):                  case DegreesOfFreedom:
1738      case(Elements):                  case Elements:
1739      case(ReducedElements):                  case ReducedElements:
1740      case(FaceElements):                  case FaceElements:
1741      case(ReducedFaceElements):                  case ReducedFaceElements:
1742      case(Points):                  case Points:
1743      case(ContactElementsZero):                  case ContactElementsZero:
1744      case(ReducedContactElementsZero):                  case ReducedContactElementsZero:
1745      case(ContactElementsOne):                  case ContactElementsOne:
1746      case(ReducedContactElementsOne):                  case ReducedContactElementsOne:
1747      return true;                      return true;
1748      default:                  default:
1749            stringstream temp;                      stringstream temp;
1750            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;
1751            throw FinleyAdapterException(temp.str());                      throw FinleyAdapterException(temp.str());
1752     }              }
1753     break;              break;
1754     case(ReducedNodes):          case ReducedNodes:
1755      switch(functionSpaceType_target) {              switch(functionSpaceType_target) {
1756      case(ReducedNodes):                  case ReducedNodes:
1757      case(ReducedDegreesOfFreedom):                  case ReducedDegreesOfFreedom:
1758      case(Elements):                  case Elements:
1759      case(ReducedElements):                  case ReducedElements:
1760      case(FaceElements):                  case FaceElements:
1761      case(ReducedFaceElements):                  case ReducedFaceElements:
1762      case(Points):                  case Points:
1763      case(ContactElementsZero):                  case ContactElementsZero:
1764      case(ReducedContactElementsZero):                  case ReducedContactElementsZero:
1765      case(ContactElementsOne):                  case ContactElementsOne:
1766      case(ReducedContactElementsOne):                  case ReducedContactElementsOne:
1767      return true;                      return true;
1768      case(Nodes):                  case Nodes:
1769      case(DegreesOfFreedom):                  case DegreesOfFreedom:
1770      return false;                      return false;
1771      default:                  default:
1772          stringstream temp;                      stringstream temp;
1773          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;
1774          throw FinleyAdapterException(temp.str());                      throw FinleyAdapterException(temp.str());
1775     }              }
1776     break;              break;
1777     case(Elements):          case Elements:
1778      if (functionSpaceType_target==Elements) {              if (functionSpaceType_target==Elements) {
1779        return true;                  return true;
1780      } else if (functionSpaceType_target==ReducedElements) {              } else if (functionSpaceType_target==ReducedElements) {
1781        return true;                  return true;
1782          } else {              } else {
1783            return false;                  return false;
1784          }              }
1785     case(ReducedElements):          case ReducedElements:
1786      if (functionSpaceType_target==ReducedElements) {              if (functionSpaceType_target==ReducedElements) {
1787        return true;                  return true;
1788      } else {              } else {
1789            return false;                  return false;
1790      }              }
1791     case(FaceElements):          case FaceElements:
1792      if (functionSpaceType_target==FaceElements) {              if (functionSpaceType_target==FaceElements) {
1793              return true;                  return true;
1794      } else if (functionSpaceType_target==ReducedFaceElements) {              } else if (functionSpaceType_target==ReducedFaceElements) {
1795              return true;                  return true;
1796      } else {              } else {
1797              return false;                  return false;
1798      }              }
1799     case(ReducedFaceElements):          case ReducedFaceElements:
1800      if (functionSpaceType_target==ReducedFaceElements) {              if (functionSpaceType_target==ReducedFaceElements) {
1801              return true;                  return true;
1802      } else {              } else {
1803          return false;                  return false;
1804      }              }
1805     case(Points):          case Points:
1806      if (functionSpaceType_target==Points) {              if (functionSpaceType_target==Points) {
1807              return true;                  return true;
1808      } else {              } else {
1809              return false;                  return false;
1810      }              }
1811     case(ContactElementsZero):          case ContactElementsZero:
1812     case(ContactElementsOne):          case ContactElementsOne:
1813      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {              if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1814              return true;                  return true;
1815      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {              } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1816              return true;                  return true;
1817      } else {              } else {
1818              return false;                  return false;
1819      }              }
1820     case(ReducedContactElementsZero):          case ReducedContactElementsZero:
1821     case(ReducedContactElementsOne):          case ReducedContactElementsOne:
1822      if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {              if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1823              return true;                  return true;
1824      } else {              } else {
1825              return false;                  return false;
1826      }              }
1827     case(DegreesOfFreedom):          case DegreesOfFreedom:
1828      switch(functionSpaceType_target) {              switch(functionSpaceType_target) {
1829      case(ReducedDegreesOfFreedom):                  case ReducedDegreesOfFreedom:
1830      case(DegreesOfFreedom):                  case DegreesOfFreedom:
1831      case(Nodes):                  case Nodes:
1832      case(ReducedNodes):                  case ReducedNodes:
1833      case(Elements):                  case Elements:
1834      case(ReducedElements):                  case ReducedElements:
1835      case(Points):                  case Points:
1836      case(FaceElements):                  case FaceElements:
1837      case(ReducedFaceElements):                  case ReducedFaceElements:
1838      case(ContactElementsZero):                  case ContactElementsZero:
1839      case(ReducedContactElementsZero):                  case ReducedContactElementsZero:
1840      case(ContactElementsOne):                  case ContactElementsOne:
1841      case(ReducedContactElementsOne):                  case ReducedContactElementsOne:
1842      return true;                      return true;
1843      default:                  default:
1844          stringstream temp;                      stringstream temp;
1845          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;
1846          throw FinleyAdapterException(temp.str());                      throw FinleyAdapterException(temp.str());
1847      }              }
1848      break;              break;
1849     case(ReducedDegreesOfFreedom):          case ReducedDegreesOfFreedom:
1850     switch(functionSpaceType_target) {              switch(functionSpaceType_target) {
1851      case(ReducedDegreesOfFreedom):                  case ReducedDegreesOfFreedom:
1852      case(ReducedNodes):                  case ReducedNodes:
1853      case(Elements):                  case Elements:
1854      case(ReducedElements):                  case ReducedElements:
1855      case(FaceElements):                  case FaceElements:
1856      case(ReducedFaceElements):                  case ReducedFaceElements:
1857      case(Points):                  case Points:
1858      case(ContactElementsZero):                  case ContactElementsZero:
1859      case(ReducedContactElementsZero):                  case ReducedContactElementsZero:
1860      case(ContactElementsOne):                  case ContactElementsOne:
1861      case(ReducedContactElementsOne):                  case ReducedContactElementsOne:
1862      return true;                      return true;
1863      case(Nodes):                  case Nodes:
1864      case(DegreesOfFreedom):                  case DegreesOfFreedom:
1865      return false;                      return false;
1866      default:                  default:
1867          stringstream temp;                      stringstream temp;
1868          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;
1869          throw FinleyAdapterException(temp.str());                      throw FinleyAdapterException(temp.str());
1870      }              }
1871      break;              break;
1872     default:          default:
1873        stringstream temp;              stringstream temp;
1874        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;
1875        throw FinleyAdapterException(temp.str());              throw FinleyAdapterException(temp.str());
1876        break;              break;
1877     }      }
1878     checkFinleyError();      return false;
1879     return false;  }
1880  }  
1881    signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source, int functionSpaceType_target) const
1882  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const  {
1883  {      if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1884     return false;          return 1;
1885  }  
1886        if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1887  bool MeshAdapter::operator==(const AbstractDomain& other) const          return -1;
1888  {  
1889     const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);      return 0;
    if (temp!=0) {  
       return (m_finleyMesh==temp->m_finleyMesh);  
    } else {  
       return false;  
    }  
1890  }  }
1891    
1892  bool MeshAdapter::operator!=(const AbstractDomain& other) const  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,
1893            const escript::AbstractDomain& targetDomain,
1894            int functionSpaceType_target) const
1895  {  {
1896     return !(operator==(other));      return false;
1897    }
1898    
1899    bool MeshAdapter::operator==(const escript::AbstractDomain& other) const
1900    {
1901        const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1902        if (temp) {
1903            return (m_finleyMesh==temp->m_finleyMesh);
1904        } else {
1905            return false;
1906        }
1907    }
1908    
1909    bool MeshAdapter::operator!=(const escript::AbstractDomain& other) const
1910    {
1911        return !(operator==(other));
1912  }  }
1913    
1914  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
1915  {  {
1916     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
1917     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);      return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
1918     checkPasoError();                 package, symmetry, mesh->MPIInfo);
    return out;  
1919  }  }
1920    
1921  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
1922  {  {
1923     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
1924     int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);      return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
1925     checkPasoError();                  package, symmetry, mesh->MPIInfo);
    return out;  
1926  }  }
1927    
1928  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
1929  {  {
1930     return continuousFunction(asAbstractContinuousDomain()).getX();      return continuousFunction(*this).getX();
1931  }  }
1932    
1933  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
1934  {  {
1935     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();      return functionOnBoundary(*this).getNormal();
1936  }  }
1937    
1938  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
1939  {  {
1940     return escript::function(asAbstractContinuousDomain()).getSize();      return escript::function(*this).getSize();
1941  }  }
1942    
1943  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1944  {  {
1945     int *out = NULL;      int *out = NULL;
1946     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
1947     switch (functionSpaceType) {      switch (functionSpaceType) {
1948     case(Nodes):          case Nodes:
1949     out=mesh->Nodes->Id;              out=mesh->Nodes->Id;
1950     break;              break;
1951     case(ReducedNodes):          case ReducedNodes:
1952     out=mesh->Nodes->reducedNodesId;              out=mesh->Nodes->reducedNodesId;
1953     break;              break;
1954     case(Elements):          case Elements:
1955     out=mesh->Elements->Id;          case ReducedElements:
1956     break;              out=mesh->Elements->Id;
1957     case(ReducedElements):              break;
1958     out=mesh->Elements->Id;          case FaceElements:
1959     break;          case ReducedFaceElements:
1960     case(FaceElements):              out=mesh->FaceElements->Id;
1961     out=mesh->FaceElements->Id;              break;
1962     break;          case Points:
1963     case(ReducedFaceElements):              out=mesh->Points->Id;
1964     out=mesh->FaceElements->Id;              break;
1965     break;          case ContactElementsZero:
1966     case(Points):          case ReducedContactElementsZero:
1967     out=mesh->Points->Id;          case ContactElementsOne:
1968     break;          case ReducedContactElementsOne:
1969     case(ContactElementsZero):              out=mesh->ContactElements->Id;
1970     out=mesh->ContactElements->Id;              break;
1971     break;          case DegreesOfFreedom:
1972     case(ReducedContactElementsZero):              out=mesh->Nodes->degreesOfFreedomId;
1973     out=mesh->ContactElements->Id;              break;
1974     break;          case ReducedDegreesOfFreedom:
1975     case(ContactElementsOne):              out=mesh->Nodes->reducedDegreesOfFreedomId;
1976     out=mesh->ContactElements->Id;              break;
1977     break;          default:
1978     case(ReducedContactElementsOne):              stringstream temp;
1979     out=mesh->ContactElements->Id;              temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1980     break;              throw FinleyAdapterException(temp.str());
1981     case(DegreesOfFreedom):              break;
1982     out=mesh->Nodes->degreesOfFreedomId;      }
1983     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;  
1984  }  }
1985  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1986  {  {
1987     int out=0;      int out=0;
1988     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
1989     switch (functionSpaceType) {      switch (functionSpaceType) {
1990     case(Nodes):          case Nodes:
1991     out=mesh->Nodes->Tag[sampleNo];              out=mesh->Nodes->Tag[sampleNo];
1992     break;              break;
1993     case(ReducedNodes):          case ReducedNodes:
1994     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");              throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1995     break;              break;
1996     case(Elements):          case Elements:
1997     out=mesh->Elements->Tag[sampleNo];          case ReducedElements:
1998     break;              out=mesh->Elements->Tag[sampleNo];
1999     case(ReducedElements):              break;
2000     out=mesh->Elements->Tag[sampleNo];          case FaceElements:
2001     break;          case ReducedFaceElements:
2002     case(FaceElements):              out=mesh->FaceElements->Tag[sampleNo];
2003     out=mesh->FaceElements->Tag[sampleNo];              break;
2004     break;          case Points:
2005     case(ReducedFaceElements):              out=mesh->Points->Tag[sampleNo];
2006     out=mesh->FaceElements->Tag[sampleNo];              break;
2007     break;          case ContactElementsZero:
2008     case(Points):          case ReducedContactElementsZero:
2009     out=mesh->Points->Tag[sampleNo];          case ContactElementsOne:
2010     break;          case ReducedContactElementsOne:
2011     case(ContactElementsZero):              out=mesh->ContactElements->Tag[sampleNo];
2012     out=mesh->ContactElements->Tag[sampleNo];              break;
2013     break;          case DegreesOfFreedom:
2014     case(ReducedContactElementsZero):              throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2015     out=mesh->ContactElements->Tag[sampleNo];              break;
2016     break;          case ReducedDegreesOfFreedom:
2017     case(ContactElementsOne):              throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2018     out=mesh->ContactElements->Tag[sampleNo];              break;
2019     break;          default:
2020     case(ReducedContactElementsOne):              stringstream temp;
2021     out=mesh->ContactElements->Tag[sampleNo];              temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2022     break;              throw FinleyAdapterException(temp.str());
2023     case(DegreesOfFreedom):              break;
2024     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");      }
2025     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;  
2026  }  }
2027    
2028    
2029  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
2030  {  {
2031     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2032     escriptDataC tmp=mask.getDataC();      escriptDataC tmp=mask.getDataC();
2033     switch(functionSpaceType) {      switch(functionSpaceType) {
2034     case(Nodes):          case Nodes:
2035     Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);              mesh->Nodes->setTags(newTag, mask);
2036     break;              break;
2037     case(ReducedNodes):          case ReducedNodes:
2038     throw FinleyAdapterException("Error - ReducedNodes does not support tags");              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2039     break;              break;
2040     case(DegreesOfFreedom):          case DegreesOfFreedom:
2041     throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2042     break;              break;
2043     case(ReducedDegreesOfFreedom):          case ReducedDegreesOfFreedom:
2044     throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2045     break;              break;
2046     case(Elements):          case Elements:
2047     Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);          case ReducedElements:
2048     break;              Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2049     case(ReducedElements):              break;
2050     Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);          case FaceElements:
2051     break;          case ReducedFaceElements:
2052     case(FaceElements):              Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2053     Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);              break;
2054     break;          case Points:
2055     case(ReducedFaceElements):              Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2056     Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);              break;
2057     break;          case ContactElementsZero:
2058     case(Points):          case ReducedContactElementsZero:
2059     Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);          case ContactElementsOne:
2060     break;          case ReducedContactElementsOne:
2061     case(ContactElementsZero):              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2062     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);              break;
2063     break;          default:
2064     case(ReducedContactElementsZero):              stringstream temp;
2065     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);              temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2066     break;              throw FinleyAdapterException(temp.str());
2067     case(ContactElementsOne):      }
2068     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);      checkFinleyError();
2069     break;  }
2070     case(ReducedContactElementsOne):  
2071     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);  void MeshAdapter::setTagMap(const string& name, int tag)
2072     break;  {
2073     default:      Finley_Mesh* mesh=m_finleyMesh.get();
2074        stringstream temp;      Finley_Mesh_addTagMap(mesh, name.c_str(), tag);
2075        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);  
    checkPasoError();  
    // throwStandardException("MeshAdapter::set TagMap is not implemented.");  
2076  }  }
2077    
2078  int MeshAdapter::getTag(const string& name) const  int MeshAdapter::getTag(const string& name) const
2079  {  {
2080     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2081     int tag=0;      int tag = Finley_Mesh_getTag(mesh, name.c_str());
2082     tag=Finley_Mesh_getTag(mesh, name.c_str());      checkFinleyError();
2083     checkPasoError();      return tag;
    // throwStandardException("MeshAdapter::getTag is not implemented.");  
    return tag;  
2084  }  }
2085    
2086  bool MeshAdapter::isValidTagName(const string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2087  {  {
2088     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2089     return Finley_Mesh_isValidTagName(mesh,name.c_str());      return Finley_Mesh_isValidTagName(mesh, name.c_str());
2090  }  }
2091    
2092  string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2093  {  {
2094     stringstream temp;      stringstream temp;
2095     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2096     Finley_TagMap* tag_map=mesh->TagMap;      Finley_TagMap* tag_map=mesh->TagMap;
2097     while (tag_map) {      while (tag_map) {
2098        temp << tag_map->name;          temp << tag_map->name;
2099        tag_map=tag_map->next;          tag_map=tag_map->next;
2100        if (tag_map) temp << ", ";          if (tag_map) temp << ", ";
2101     }      }
2102     return temp.str();      return temp.str();
2103  }  }
2104    
2105  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2106  {  {
2107    Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2108    dim_t numTags=0;      switch(functionSpaceCode) {
2109    switch(functionSpaceCode) {          case Nodes:
2110     case(Nodes):              return mesh->Nodes->tagsInUse.size();
2111            numTags=mesh->Nodes->numTagsInUse;          case ReducedNodes:
2112            break;              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2113     case(ReducedNodes):          case DegreesOfFreedom:
2114            throw FinleyAdapterException("Error - ReducedNodes does not support tags");              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2115            break;          case ReducedDegreesOfFreedom:
2116     case(DegreesOfFreedom):              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2117            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");          case Elements:
2118            break;          case ReducedElements:
2119     case(ReducedDegreesOfFreedom):              return mesh->Elements->tagsInUse.size();
2120            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");          case FaceElements:
2121            break;          case ReducedFaceElements:
2122     case(Elements):              return mesh->FaceElements->tagsInUse.size();
2123     case(ReducedElements):          case Points:
2124            numTags=mesh->Elements->numTagsInUse;              return mesh->Points->tagsInUse.size();
2125            break;          case ContactElementsZero:
2126     case(FaceElements):          case ReducedContactElementsZero:
2127     case(ReducedFaceElements):          case ContactElementsOne:
2128            numTags=mesh->FaceElements->numTagsInUse;          case ReducedContactElementsOne:
2129            break;              return mesh->ContactElements->tagsInUse.size();
2130     case(Points):          default:
2131            numTags=mesh->Points->numTagsInUse;              stringstream ss;
2132            break;              ss << "Finley does not know anything about function space type "
2133     case(ContactElementsZero):                   << functionSpaceCode;
2134     case(ReducedContactElementsZero):              throw FinleyAdapterException(ss.str());
2135     case(ContactElementsOne):      }
2136     case(ReducedContactElementsOne):      return 0;
           numTags=mesh->ContactElements->numTagsInUse;  
           break;  
    default:  
       stringstream temp;  
       temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;  
       throw FinleyAdapterException(temp.str());  
   }  
   return numTags;  
2137  }  }
2138    
2139  const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2140  {  {
2141    Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2142    index_t* tags=NULL;      switch(functionSpaceCode) {
2143    switch(functionSpaceCode) {          case Nodes:
2144     case(Nodes):              return &mesh->Nodes->tagsInUse[0];
2145            tags=mesh->Nodes->tagsInUse;          case ReducedNodes:
2146            break;              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2147     case(ReducedNodes):          case DegreesOfFreedom:
2148            throw FinleyAdapterException("Error - ReducedNodes does not support tags");              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2149            break;          case ReducedDegreesOfFreedom:
2150     case(DegreesOfFreedom):              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2151            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");          case Elements:
2152            break;          case ReducedElements:
2153     case(ReducedDegreesOfFreedom):              return &mesh->Elements->tagsInUse[0];
2154            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");          case FaceElements:
2155            break;          case ReducedFaceElements:
2156     case(Elements):              return &mesh->FaceElements->tagsInUse[0];
2157     case(ReducedElements):          case Points:
2158            tags=mesh->Elements->tagsInUse;              return &mesh->Points->tagsInUse[0];
2159            break;          case ContactElementsZero:
2160     case(FaceElements):          case ReducedContactElementsZero:
2161     case(ReducedFaceElements):          case ContactElementsOne:
2162            tags=mesh->FaceElements->tagsInUse;          case ReducedContactElementsOne:
2163            break;              return &mesh->ContactElements->tagsInUse[0];
2164     case(Points):          default:
2165            tags=mesh->Points->tagsInUse;              stringstream temp;
2166            break;              temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2167     case(ContactElementsZero):              throw FinleyAdapterException(temp.str());
2168     case(ReducedContactElementsZero):      }
2169     case(ContactElementsOne):      return NULL;
    case(ReducedContactElementsOne):  
           tags=mesh->ContactElements->tagsInUse;  
           break;  
    default:  
       stringstream temp;  
       temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;  
       throw FinleyAdapterException(temp.str());  
   }  
   return tags;  
2170  }  }
2171    
2172    
2173  bool MeshAdapter::canTag(int functionSpaceCode) const  bool MeshAdapter::canTag(int functionSpaceCode) const
2174  {  {
2175    switch(functionSpaceCode) {      switch(functionSpaceCode) {
2176     case(Nodes):          case Nodes:
2177     case(Elements):          case Elements:
2178     case(ReducedElements):          case ReducedElements:
2179     case(FaceElements):          case FaceElements:
2180     case(ReducedFaceElements):          case ReducedFaceElements:
2181     case(Points):          case Points:
2182     case(ContactElementsZero):          case ContactElementsZero:
2183     case(ReducedContactElementsZero):          case ReducedContactElementsZero:
2184     case(ContactElementsOne):          case ContactElementsOne:
2185     case(ReducedContactElementsOne):          case ReducedContactElementsOne:
2186            return true;              return true;
2187     case(ReducedNodes):          default:
2188     case(DegreesOfFreedom):              return false;
2189     case(ReducedDegreesOfFreedom):      }
       return false;  
    default:  
     return false;  
   }  
2190  }  }
2191    
2192  AbstractDomain::StatusType MeshAdapter::getStatus() const  escript::AbstractDomain::StatusType MeshAdapter::getStatus() const
2193  {  {
2194    Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2195    return Finley_Mesh_getStatus(mesh);      return Finley_Mesh_getStatus(mesh);
2196  }  }
2197    
2198  int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const  int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2199  {  {
2200          Finley_Mesh* mesh=m_finleyMesh.get();
2201    Finley_Mesh* mesh=m_finleyMesh.get();      int order =-1;
2202    int order =-1;      switch(functionSpaceCode) {
2203    switch(functionSpaceCode) {          case Nodes:
2204     case(Nodes):          case DegreesOfFreedom:
2205     case(DegreesOfFreedom):              order=mesh->approximationOrder;
2206            order=mesh->approximationOrder;              break;
2207            break;          case ReducedNodes:
2208     case(ReducedNodes):          case ReducedDegreesOfFreedom:
2209     case(ReducedDegreesOfFreedom):              order=mesh->reducedApproximationOrder;
2210            order=mesh->reducedApproximationOrder;              break;
2211            break;          case Elements:
2212     case(Elements):          case FaceElements:
2213     case(FaceElements):          case Points:
2214     case(Points):          case ContactElementsZero:
2215     case(ContactElementsZero):          case ContactElementsOne:
2216     case(ContactElementsOne):              order=mesh->integrationOrder;
2217            order=mesh->integrationOrder;              break;
2218            break;          case ReducedElements:
2219     case(ReducedElements):          case ReducedFaceElements:
2220     case(ReducedFaceElements):          case ReducedContactElementsZero:
2221     case(ReducedContactElementsZero):          case ReducedContactElementsOne:
2222     case(ReducedContactElementsOne):              order=mesh->reducedIntegrationOrder;
2223            order=mesh->reducedIntegrationOrder;              break;
2224            break;          default:
2225     default:              stringstream temp;
2226        stringstream temp;              temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2227        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;              throw FinleyAdapterException(temp.str());
2228        throw FinleyAdapterException(temp.str());      }
2229    }      return order;
2230    return order;  }
2231    
2232    bool MeshAdapter::supportsContactElements() const
2233    {
2234        return true;
2235  }  }
2236    
2237  ReferenceElementSetWrapper::ReferenceElementSetWrapper(ElementTypeId id, index_t order, index_t reducedOrder)  ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id,
2238            index_t order, index_t reducedOrder)
2239  {  {
2240    m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2241  }  }
2242    
2243  ReferenceElementSetWrapper::~ReferenceElementSetWrapper()  ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2244  {  {
2245    Finley_ReferenceElementSet_dealloc(m_refSet);      Finley_ReferenceElementSet_dealloc(m_refSet);
2246    }
2247    
2248    void MeshAdapter::addDiracPoints(const vector<double>& points,
2249                                     const vector<int>& tags) const
2250    {
2251        // points will be flattened
2252        const int dim = getDim();
2253        int numPoints=points.size()/dim;
2254        int numTags=tags.size();
2255        Finley_Mesh* mesh=m_finleyMesh.get();
2256    
2257        if ( points.size() % dim != 0 ) {
2258            throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2259        }
2260    
2261        if ((numTags > 0) && (numPoints != numTags))
2262            throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2263    
2264        Finley_Mesh_addPoints(mesh, numPoints, &points[0], &tags[0]);
2265        checkFinleyError();
2266  }  }
2267    
2268    // void MeshAdapter::addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2269    // {
2270    //       const int dim = getDim();
2271    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2272    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2273    //       Finley_Mesh* mesh=m_finleyMesh.get();
2274    //
2275    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2276    //       throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2277    //
2278    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2279    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2280    //
2281    //       for (int i=0;i<numPoints;++i) {
2282    //         int tag_id=-1;
2283    //         int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2284    //         if  ( numComps !=   dim ) {
2285    //                stringstream temp;
2286    //                temp << "Error - illegal number of components " << numComps << " for point " << i;
2287    //                throw FinleyAdapterException(temp.str());
2288    //         }
2289    //         points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2290    //         if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2291    //         if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2292    //
2293    //         if ( numTags > 0) {
2294    //                boost::python::extract<string> ex_str(tags[i]);
2295    //                if  ( ex_str.check() ) {
2296    //                    tag_id=getTag( ex_str());
2297    //                } else {
2298    //                     boost::python::extract<int> ex_int(tags[i]);
2299    //                     if ( ex_int.check() ) {
2300    //                         tag_id=ex_int();
2301    //                     } else {
2302    //                          stringstream temp;
2303    //                          temp << "Error - unable to extract tag for point " << i;
2304    //                          throw FinleyAdapterException(temp.str());
2305    //                    }
2306    //                }
2307    //         }
2308    //            tags_ptr[i]=tag_id;
2309    //       }
2310    //
2311    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2312    //       checkPasoError();
2313    //
2314    //       TMPMEMFREE(points_ptr);
2315    //       TMPMEMFREE(tags_ptr);
2316    // }
2317    
2318    /*
2319    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2320    {
2321        boost::python::list points =  boost::python::list();
2322        boost::python::list tags =  boost::python::list();
2323        points.append(point);
2324        tags.append(tag);
2325        addDiracPoints(points, tags);
2326    }
2327    */
2328    
2329    /*
2330    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const string& tag) const
2331    {
2332        boost::python::list points =   boost::python::list();
2333        boost::python::list tags =   boost::python::list();
2334        points.append(point);
2335        tags.append(tag);
2336        addDiracPoints(points, tags);
2337    }
2338    */
2339  }  // end of namespace  }  // end of namespace
2340    

Legend:
Removed from v.2989  
changed lines
  Added in v.4428

  ViewVC Help
Powered by ViewVC 1.1.26