/[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 4410 by caltinay, Tue May 14 11:29:23 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 Finley_NodeFile_getGlobalNumNodes(m_finleyMesh.get()->Nodes);
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=Finley_NodeFile_getNumNodes(mesh->Nodes);
668     break;              break;
669     case(ReducedNodes):          case ReducedNodes:
670     numDataPointsPerSample=1;              numDataPointsPerSample=1;
671     numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);              numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
672     break;              break;
673     case(Elements):          case Elements:
674     if (mesh->Elements!=NULL) {              if (mesh->Elements!=NULL) {
675        numSamples=mesh->Elements->numElements;                  numSamples=mesh->Elements->numElements;
676        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;                  numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
677     }              }
678     break;              break;
679     case(ReducedElements):          case ReducedElements:
680     if (mesh->Elements!=NULL) {              if (mesh->Elements!=NULL) {
681        numSamples=mesh->Elements->numElements;                  numSamples=mesh->Elements->numElements;
682        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;                  numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
683     }              }
684     break;              break;
685     case(FaceElements):          case FaceElements:
686     if (mesh->FaceElements!=NULL) {              if (mesh->FaceElements!=NULL) {
687        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;                  numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
688        numSamples=mesh->FaceElements->numElements;                  numSamples=mesh->FaceElements->numElements;
689     }              }
690     break;              break;
691     case(ReducedFaceElements):          case ReducedFaceElements:
692     if (mesh->FaceElements!=NULL) {              if (mesh->FaceElements!=NULL) {
693        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;                  numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
694        numSamples=mesh->FaceElements->numElements;                  numSamples=mesh->FaceElements->numElements;
695     }              }
696     break;              break;
697     case(Points):          case Points:
698     if (mesh->Points!=NULL) {              if (mesh->Points!=NULL) {
699        numDataPointsPerSample=1;                  numDataPointsPerSample=1;
700        numSamples=mesh->Points->numElements;                  numSamples=mesh->Points->numElements;
701     }              }
702     break;              break;
703     case(ContactElementsZero):          case ContactElementsZero:
704     if (mesh->ContactElements!=NULL) {              if (mesh->ContactElements!=NULL) {
705        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;                  numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
706        numSamples=mesh->ContactElements->numElements;                  numSamples=mesh->ContactElements->numElements;
707     }              }
708     break;              break;
709     case(ReducedContactElementsZero):          case ReducedContactElementsZero:
710     if (mesh->ContactElements!=NULL) {              if (mesh->ContactElements!=NULL) {
711        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;                  numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
712        numSamples=mesh->ContactElements->numElements;                  numSamples=mesh->ContactElements->numElements;
713     }              }
714     break;              break;
715     case(ContactElementsOne):          case ContactElementsOne:
716     if (mesh->ContactElements!=NULL) {              if (mesh->ContactElements!=NULL) {
717        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;                  numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
718        numSamples=mesh->ContactElements->numElements;                  numSamples=mesh->ContactElements->numElements;
719     }              }
720     break;              break;
721     case(ReducedContactElementsOne):          case ReducedContactElementsOne:
722     if (mesh->ContactElements!=NULL) {              if (mesh->ContactElements!=NULL) {
723        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;                  numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
724        numSamples=mesh->ContactElements->numElements;                  numSamples=mesh->ContactElements->numElements;
725     }              }
726     break;              break;
727     case(DegreesOfFreedom):          case DegreesOfFreedom:
728     if (mesh->Nodes!=NULL) {              if (mesh->Nodes!=NULL) {
729        numDataPointsPerSample=1;                  numDataPointsPerSample=1;
730        numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);                  numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
731     }              }
732     break;              break;
733     case(ReducedDegreesOfFreedom):          case ReducedDegreesOfFreedom:
734     if (mesh->Nodes!=NULL) {              if (mesh->Nodes!=NULL) {
735        numDataPointsPerSample=1;                  numDataPointsPerSample=1;
736        numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);                  numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
737     }              }
738     break;              break;
739     default:          default:
740        stringstream temp;              stringstream temp;
741        temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();              temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
742        throw FinleyAdapterException(temp.str());              throw FinleyAdapterException(temp.str());
743        break;              break;
744     }      }
745     return pair<int,int>(numDataPointsPerSample,numSamples);      return pair<int,int>(numDataPointsPerSample,numSamples);
746  }  }
747    
748  //  //
749  // adds linear PDE of second order into a given stiffness matrix and right hand side:  // adds linear PDE of second order into a given stiffness matrix and right
750    // hand side:
751  //  //
752  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
753                                   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;      escriptDataC tmp;
1449     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));      const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1450     if (newDomain!=*this)      if (newDomain!=*this)
1451        throw FinleyAdapterException("Error - Illegal domain of new point locations");          throw FinleyAdapterException("Error - Illegal domain of new point locations");
1452     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {      if (new_x.getFunctionSpace() == continuousFunction(*this)) {
1453         tmp = new_x.getDataC();          tmp = new_x.getDataC();
1454         Finley_Mesh_setCoordinates(mesh,&tmp);          Finley_Mesh_setCoordinates(mesh,&tmp);
1455     } else {      } else {
1456         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );          throw FinleyAdapterException("As of escript version 3.3 SetX() only accepts ContinuousFunction arguments. Please interpolate.");
1457         tmp = new_x_inter.getDataC();          //escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1458         Finley_Mesh_setCoordinates(mesh,&tmp);          //tmp = new_x_inter.getDataC();
1459     }          //Finley_Mesh_setCoordinates(mesh,&tmp);
1460     checkFinleyError();      }
1461  }      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);  
1462  }  }
1463    
1464  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1465  {  {
1466  #ifdef PASO_MPI      if (getMPISize() > 1) {
1467      index_t myFirstNode=0, myLastNode=0, k=0;  #ifdef ESYS_MPI
1468      index_t* globalNodeIndex=0;          index_t myFirstNode=0, myLastNode=0, k=0;
1469      Finley_Mesh* mesh_p=m_finleyMesh.get();          index_t* globalNodeIndex=0;
1470      if (fs_code == FINLEY_REDUCED_NODES)          Finley_Mesh* mesh_p=m_finleyMesh.get();
1471      {          /*
1472      myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);           * this method is only used by saveDataCSV which would use the returned
1473      myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);           * values for reduced nodes wrongly so this case is disabled for now
1474      globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);          if (fs_code == FINLEY_REDUCED_NODES) {
1475      }              myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1476      else              myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1477      {              globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1478      myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);          } else
1479      myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);          */
1480      globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);          if (fs_code == FINLEY_NODES) {
1481      }              myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1482      k=globalNodeIndex[id];              myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1483      return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );              globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1484            } else {
1485                throw FinleyAdapterException("Unsupported function space type for ownSample()");
1486            }
1487    
1488            k=globalNodeIndex[id];
1489            return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1490  #endif  #endif
1491        }
1492      return true;      return true;
1493  }  }
1494    
1495    
   
1496  //  //
1497  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1498  //  //
1499  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  escript::ASM_ptr MeshAdapter::newSystemMatrix(const int row_blocksize,
1500                                                   const int row_blocksize,                              const escript::FunctionSpace& row_functionspace,
1501                                                   const escript::FunctionSpace& row_functionspace,                              const int column_blocksize,
1502                                                   const int column_blocksize,                              const escript::FunctionSpace& column_functionspace,
1503                                                   const escript::FunctionSpace& column_functionspace,                              const int type) const
1504                                                   const int type) const  {
1505  {      // is the domain right?
1506     int reduceRowOrder=0;      const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1507     int reduceColOrder=0;      if (row_domain!=*this)
1508     // is the domain right?          throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1509     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));      const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1510     if (row_domain!=*this)      if (col_domain!=*this)
1511        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.");
1512     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));  
1513     if (col_domain!=*this)      int reduceRowOrder=0;
1514        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");      int reduceColOrder=0;
1515     // is the function space type right      // is the function space type right?
1516     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {      if (row_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1517        reduceRowOrder=0;          reduceRowOrder=1;
1518     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      } else if (row_functionspace.getTypeCode() != DegreesOfFreedom) {
1519        reduceRowOrder=1;          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1520     } else {      }
1521        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");      if (column_functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1522     }          reduceColOrder=1;
1523     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {      } else if (column_functionspace.getTypeCode() != DegreesOfFreedom) {
1524        reduceColOrder=0;          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1525     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      }
1526        reduceColOrder=1;  
1527     } else {      // generate matrix:
1528        throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(
1529     }              getFinley_Mesh(), reduceRowOrder, reduceColOrder);
1530     // generate matrix:      checkFinleyError();
1531        Paso_SystemMatrix* fsystemMatrix;
1532     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);      const int trilinos = 0;
1533     checkFinleyError();      if (trilinos) {
    Paso_SystemMatrix* fsystemMatrix;  
    int trilinos = 0;  
    if (trilinos) {  
1534  #ifdef TRILINOS  #ifdef TRILINOS
1535        /* Allocation Epetra_VrbMatrix here */          // FIXME: Allocation Epetra_VrbMatrix here...
1536  #endif  #endif
1537     }      } else {
1538     else {          fsystemMatrix=Paso_SystemMatrix_alloc(type, fsystemMatrixPattern,
1539        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);                  row_blocksize, column_blocksize, FALSE);
1540     }      }
1541     checkPasoError();      checkPasoError();
1542     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1543     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);      SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1544        return escript::ASM_ptr(sma);
1545  }  }
1546    
1547  //  //
1548  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1549  //  //
1550  TransportProblemAdapter MeshAdapter::newTransportProblem(  escript::ATP_ptr MeshAdapter::newTransportProblem(const int blocksize,
1551                                                           const bool useBackwardEuler,          const escript::FunctionSpace& functionspace, const int type) const
1552                                                           const int blocksize,  {
1553                                                           const escript::FunctionSpace& functionspace,      // is the domain right?
1554                                                           const int type) const      const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1555  {      if (domain!=*this)
1556     int reduceOrder=0;          throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1557     // is the domain right?  
1558     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));      // is the function space type right?
1559     if (domain!=*this)      int reduceOrder=0;
1560        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");      if (functionspace.getTypeCode() == ReducedDegreesOfFreedom) {
1561     // is the function space type right          reduceOrder=1;
1562     if (functionspace.getTypeCode()==DegreesOfFreedom) {      } else if (functionspace.getTypeCode() != DegreesOfFreedom) {
1563        reduceOrder=0;          throw FinleyAdapterException("Error - illegal function space type for transport problem.");
1564     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      }
1565        reduceOrder=1;  
1566     } else {      // generate transport problem:
1567        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(
1568     }              getFinley_Mesh(), reduceOrder, reduceOrder);
1569     // generate matrix:      checkFinleyError();
1570        Paso_TransportProblem* transportProblem;
1571     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);      transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern, blocksize);
1572     checkFinleyError();      checkPasoError();
1573     Paso_TransportProblem* transportProblem;      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1574     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);      TransportProblemAdapter* tpa=new TransportProblemAdapter(
1575     checkPasoError();              transportProblem, blocksize, functionspace);
1576     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);      return escript::ATP_ptr(tpa);
    return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);  
1577  }  }
1578    
1579  //  //
1580  // vtkObject MeshAdapter::createVtkObject() const  // returns true if data on functionSpaceCode is considered as being cell centered
 // TODO:  
1581  //  //
   
 //  
 // returns true if data at the atom_type is considered as being cell centered:  
1582  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1583  {  {
1584     switch(functionSpaceCode) {      switch(functionSpaceCode) {
1585     case(Nodes):          case Nodes:
1586     case(DegreesOfFreedom):          case DegreesOfFreedom:
1587     case(ReducedDegreesOfFreedom):          case ReducedDegreesOfFreedom:
1588     return false;              return false;
1589     break;          case Elements:
1590     case(Elements):          case FaceElements:
1591     case(FaceElements):          case Points:
1592     case(Points):          case ContactElementsZero:
1593     case(ContactElementsZero):          case ContactElementsOne:
1594     case(ContactElementsOne):          case ReducedElements:
1595     case(ReducedElements):          case ReducedFaceElements:
1596     case(ReducedFaceElements):          case ReducedContactElementsZero:
1597     case(ReducedContactElementsZero):          case ReducedContactElementsOne:
1598     case(ReducedContactElementsOne):              return true;
1599     return true;          default:
1600     break;              stringstream temp;
1601     default:              temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1602        stringstream temp;              throw FinleyAdapterException(temp.str());
1603        temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;              break;
1604        throw FinleyAdapterException(temp.str());      }
1605        break;      return false;
    }  
    checkFinleyError();  
    return false;  
1606  }  }
1607    
1608  bool  bool
1609  MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const  MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1610  {  {
1611     /* 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
1612      class 1: DOF <-> Nodes         back and forth]:
1613      class 2: ReducedDOF <-> ReducedNodes          class 1: DOF <-> Nodes
1614      class 3: Points          class 2: ReducedDOF <-> ReducedNodes
1615      class 4: Elements          class 3: Points
1616      class 5: ReducedElements          class 4: Elements
1617      class 6: FaceElements          class 5: ReducedElements
1618      class 7: ReducedFaceElements          class 6: FaceElements
1619      class 8: ContactElementZero <-> ContactElementOne          class 7: ReducedFaceElements
1620      class 9: ReducedContactElementZero <-> ReducedContactElementOne          class 8: ContactElementZero <-> ContactElementOne
1621            class 9: ReducedContactElementZero <-> ReducedContactElementOne
1622     There is also a set of lines. Interpolation is possible down a line but not between lines.  
1623     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
1624      line 0: class 3      not between lines.
1625      line 1: class 4,5      class 1 and 2 belong to all lines so aren't considered.
1626      line 2: class 6,7          line 0: class 3
1627      line 3: class 8,9          line 1: class 4,5
1628            line 2: class 6,7
1629     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
1630     eg hasnodes is true if we have at least one instance of Nodes.  
1631     */      For classes with multiple members (eg class 2) we have vars to record if
1632        there is at least one instance.
1633        e.g. hasnodes is true if we have at least one instance of Nodes.
1634        */
1635      if (fs.empty())      if (fs.empty())
     {  
1636          return false;          return false;
1637      }  
1638      vector<int> hasclass(10);      vector<int> hasclass(10);
1639      vector<int> hasline(4);      vector<int> hasline(4);    
1640      bool hasnodes=false;      bool hasnodes=false;
1641      bool hasrednodes=false;      bool hasrednodes=false;
1642      bool hascez=false;      bool hascez=false;
1643      bool hasrcez=false;      bool hasrcez=false;
1644      for (int i=0;i<fs.size();++i)      for (int i=0;i<fs.size();++i) {
1645      {          switch(fs[i]) {
1646      switch(fs[i])              case Nodes:
1647      {                  hasnodes=true; // fall through
1648      case(Nodes):   hasnodes=true;   // no break is deliberate              case DegreesOfFreedom:
1649      case(DegreesOfFreedom):                  hasclass[1]=1;
1650          hasclass[1]=1;                  break;
1651          break;              case ReducedNodes:
1652      case(ReducedNodes):    hasrednodes=true;    // no break is deliberate                  hasrednodes=true; // fall through
1653      case(ReducedDegreesOfFreedom):              case ReducedDegreesOfFreedom:
1654          hasclass[2]=1;                  hasclass[2]=1;
1655          break;                  break;
1656      case(Points):              case Points:
1657          hasline[0]=1;                  hasline[0]=1;
1658          hasclass[3]=1;                  hasclass[3]=1;
1659          break;                  break;
1660      case(Elements):              case Elements:
1661          hasclass[4]=1;                  hasclass[4]=1;
1662          hasline[1]=1;                  hasline[1]=1;
1663          break;                  break;
1664      case(ReducedElements):              case ReducedElements:
1665          hasclass[5]=1;                  hasclass[5]=1;
1666          hasline[1]=1;                  hasline[1]=1;
1667          break;                  break;
1668      case(FaceElements):              case FaceElements:
1669          hasclass[6]=1;                  hasclass[6]=1;
1670          hasline[2]=1;                  hasline[2]=1;
1671          break;                  break;
1672      case(ReducedFaceElements):              case ReducedFaceElements:
1673          hasclass[7]=1;                  hasclass[7]=1;
1674          hasline[2]=1;                  hasline[2]=1;
1675          break;                  break;
1676      case(ContactElementsZero):  hascez=true;    // no break is deliberate              case ContactElementsZero:
1677      case(ContactElementsOne):                  hascez=true; // fall through
1678          hasclass[8]=1;              case ContactElementsOne:
1679          hasline[3]=1;                  hasclass[8]=1;
1680          break;                  hasline[3]=1;
1681      case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate                  break;
1682      case(ReducedContactElementsOne):              case ReducedContactElementsZero:
1683          hasclass[9]=1;                  hasrcez=true; // fall through
1684          hasline[3]=1;              case ReducedContactElementsOne:
1685          break;                  hasclass[9]=1;
1686      default:                  hasline[3]=1;
1687          return false;                  break;
1688      }              default:
1689                    return false;
1690            }
1691      }      }
1692      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  
1693    
1694        // fail if we have more than one leaf group
1695      if (totlines>1)      if (totlines>1)
1696      {          return false; // there are at least two branches we can't interpolate between
1697      return false;   // there are at least two branches we can't interpolate between      else if (totlines==1) {
1698      }          if (hasline[0]==1)              // we have points
1699      else if (totlines==1)              resultcode=Points;
1700      {          else if (hasline[1]==1) {
1701      if (hasline[0]==1)      // we have points              if (hasclass[5]==1)
1702      {                  resultcode=ReducedElements;
1703          resultcode=Points;              else
1704      }                  resultcode=Elements;
1705      else if (hasline[1]==1)          } else if (hasline[2]==1) {
1706      {              if (hasclass[7]==1)
1707          if (hasclass[5]==1)                  resultcode=ReducedFaceElements;
1708          {              else
1709          resultcode=ReducedElements;                  resultcode=FaceElements;
1710          }          } else {   // so we must be in line3
1711          else              if (hasclass[9]==1) {
1712          {                  // need something from class 9
1713          resultcode=Elements;                  resultcode=(hasrcez ? ReducedContactElementsZero : ReducedContactElementsOne);
1714          }              } else {
1715      }                  // something from class 8
1716      else if (hasline[2]==1)                  resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1717      {              }
1718          if (hasclass[7]==1)          }
1719          {      } else { // totlines==0
1720          resultcode=ReducedFaceElements;          if (hasclass[2]==1) {
1721          }              // something from class 2
1722          else              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
1723          {          } else {
1724          resultcode=FaceElements;              // something from class 1
1725          }              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
1726      }          }
     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);  
     }  
1727      }      }
1728      return true;      return true;
1729  }  }
1730    
1731  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,
1732                                                 int functionSpaceType_target) const
1733  {  {
1734     switch(functionSpaceType_source) {      switch(functionSpaceType_source) {
1735     case(Nodes):          case Nodes:
1736      switch(functionSpaceType_target) {              switch(functionSpaceType_target) {
1737      case(Nodes):                  case Nodes:
1738      case(ReducedNodes):                  case ReducedNodes:
1739      case(ReducedDegreesOfFreedom):                  case ReducedDegreesOfFreedom:
1740      case(DegreesOfFreedom):                  case DegreesOfFreedom:
1741      case(Elements):                  case Elements:
1742      case(ReducedElements):                  case ReducedElements:
1743      case(FaceElements):                  case FaceElements:
1744      case(ReducedFaceElements):                  case ReducedFaceElements:
1745      case(Points):                  case Points:
1746      case(ContactElementsZero):                  case ContactElementsZero:
1747      case(ReducedContactElementsZero):                  case ReducedContactElementsZero:
1748      case(ContactElementsOne):                  case ContactElementsOne:
1749      case(ReducedContactElementsOne):                  case ReducedContactElementsOne:
1750      return true;                      return true;
1751      default:                  default:
1752            stringstream temp;                      stringstream temp;
1753            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;
1754            throw FinleyAdapterException(temp.str());                      throw FinleyAdapterException(temp.str());
1755     }              }
1756     break;              break;
1757     case(ReducedNodes):          case ReducedNodes:
1758      switch(functionSpaceType_target) {              switch(functionSpaceType_target) {
1759      case(ReducedNodes):                  case ReducedNodes:
1760      case(ReducedDegreesOfFreedom):                  case ReducedDegreesOfFreedom:
1761      case(Elements):                  case Elements:
1762      case(ReducedElements):                  case ReducedElements:
1763      case(FaceElements):                  case FaceElements:
1764      case(ReducedFaceElements):                  case ReducedFaceElements:
1765      case(Points):                  case Points:
1766      case(ContactElementsZero):                  case ContactElementsZero:
1767      case(ReducedContactElementsZero):                  case ReducedContactElementsZero:
1768      case(ContactElementsOne):                  case ContactElementsOne:
1769      case(ReducedContactElementsOne):                  case ReducedContactElementsOne:
1770      return true;                      return true;
1771      case(Nodes):                  case Nodes:
1772      case(DegreesOfFreedom):                  case DegreesOfFreedom:
1773      return false;                      return false;
1774      default:                  default:
1775          stringstream temp;                      stringstream temp;
1776          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;                      temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1777          throw FinleyAdapterException(temp.str());                      throw FinleyAdapterException(temp.str());
1778     }              }
1779     break;              break;
1780     case(Elements):          case Elements:
1781      if (functionSpaceType_target==Elements) {              if (functionSpaceType_target==Elements) {
1782        return true;                  return true;
1783      } else if (functionSpaceType_target==ReducedElements) {              } else if (functionSpaceType_target==ReducedElements) {
1784        return true;                  return true;
1785          } else {              } else {
1786            return false;                  return false;
1787          }              }
1788     case(ReducedElements):          case ReducedElements:
1789      if (functionSpaceType_target==ReducedElements) {              if (functionSpaceType_target==ReducedElements) {
1790        return true;                  return true;
1791      } else {              } else {
1792            return false;                  return false;
1793      }              }
1794     case(FaceElements):          case FaceElements:
1795      if (functionSpaceType_target==FaceElements) {              if (functionSpaceType_target==FaceElements) {
1796              return true;                  return true;
1797      } else if (functionSpaceType_target==ReducedFaceElements) {              } else if (functionSpaceType_target==ReducedFaceElements) {
1798              return true;                  return true;
1799      } else {              } else {
1800              return false;                  return false;
1801      }              }
1802     case(ReducedFaceElements):          case ReducedFaceElements:
1803      if (functionSpaceType_target==ReducedFaceElements) {              if (functionSpaceType_target==ReducedFaceElements) {
1804              return true;                  return true;
1805      } else {              } else {
1806          return false;                  return false;
1807      }              }
1808     case(Points):          case Points:
1809      if (functionSpaceType_target==Points) {              if (functionSpaceType_target==Points) {
1810              return true;                  return true;
1811      } else {              } else {
1812              return false;                  return false;
1813      }              }
1814     case(ContactElementsZero):          case ContactElementsZero:
1815     case(ContactElementsOne):          case ContactElementsOne:
1816      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {              if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1817              return true;                  return true;
1818      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {              } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1819              return true;                  return true;
1820      } else {              } else {
1821              return false;                  return false;
1822      }              }
1823     case(ReducedContactElementsZero):          case ReducedContactElementsZero:
1824     case(ReducedContactElementsOne):          case ReducedContactElementsOne:
1825      if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {              if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1826              return true;                  return true;
1827      } else {              } else {
1828              return false;                  return false;
1829      }              }
1830     case(DegreesOfFreedom):          case DegreesOfFreedom:
1831      switch(functionSpaceType_target) {              switch(functionSpaceType_target) {
1832      case(ReducedDegreesOfFreedom):                  case ReducedDegreesOfFreedom:
1833      case(DegreesOfFreedom):                  case DegreesOfFreedom:
1834      case(Nodes):                  case Nodes:
1835      case(ReducedNodes):                  case ReducedNodes:
1836      case(Elements):                  case Elements:
1837      case(ReducedElements):                  case ReducedElements:
1838      case(Points):                  case Points:
1839      case(FaceElements):                  case FaceElements:
1840      case(ReducedFaceElements):                  case ReducedFaceElements:
1841      case(ContactElementsZero):                  case ContactElementsZero:
1842      case(ReducedContactElementsZero):                  case ReducedContactElementsZero:
1843      case(ContactElementsOne):                  case ContactElementsOne:
1844      case(ReducedContactElementsOne):                  case ReducedContactElementsOne:
1845      return true;                      return true;
1846      default:                  default:
1847          stringstream temp;                      stringstream temp;
1848          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;
1849          throw FinleyAdapterException(temp.str());                      throw FinleyAdapterException(temp.str());
1850      }              }
1851      break;              break;
1852     case(ReducedDegreesOfFreedom):          case ReducedDegreesOfFreedom:
1853     switch(functionSpaceType_target) {              switch(functionSpaceType_target) {
1854      case(ReducedDegreesOfFreedom):                  case ReducedDegreesOfFreedom:
1855      case(ReducedNodes):                  case ReducedNodes:
1856      case(Elements):                  case Elements:
1857      case(ReducedElements):                  case ReducedElements:
1858      case(FaceElements):                  case FaceElements:
1859      case(ReducedFaceElements):                  case ReducedFaceElements:
1860      case(Points):                  case Points:
1861      case(ContactElementsZero):                  case ContactElementsZero:
1862      case(ReducedContactElementsZero):                  case ReducedContactElementsZero:
1863      case(ContactElementsOne):                  case ContactElementsOne:
1864      case(ReducedContactElementsOne):                  case ReducedContactElementsOne:
1865      return true;                      return true;
1866      case(Nodes):                  case Nodes:
1867      case(DegreesOfFreedom):                  case DegreesOfFreedom:
1868      return false;                      return false;
1869      default:                  default:
1870          stringstream temp;                      stringstream temp;
1871          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;
1872          throw FinleyAdapterException(temp.str());                      throw FinleyAdapterException(temp.str());
1873      }              }
1874      break;              break;
1875     default:          default:
1876        stringstream temp;              stringstream temp;
1877        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;
1878        throw FinleyAdapterException(temp.str());              throw FinleyAdapterException(temp.str());
1879        break;              break;
1880     }      }
1881     checkFinleyError();      return false;
1882     return false;  }
1883  }  
1884    signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source, int functionSpaceType_target) const
1885  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const  {
1886  {      if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1887     return false;          return 1;
1888  }  
1889        if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1890  bool MeshAdapter::operator==(const AbstractDomain& other) const          return -1;
1891  {  
1892     const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);      return 0;
    if (temp!=0) {  
       return (m_finleyMesh==temp->m_finleyMesh);  
    } else {  
       return false;  
    }  
1893  }  }
1894    
1895  bool MeshAdapter::operator!=(const AbstractDomain& other) const  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,
1896            const escript::AbstractDomain& targetDomain,
1897            int functionSpaceType_target) const
1898  {  {
1899     return !(operator==(other));      return false;
1900    }
1901    
1902    bool MeshAdapter::operator==(const escript::AbstractDomain& other) const
1903    {
1904        const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1905        if (temp) {
1906            return (m_finleyMesh==temp->m_finleyMesh);
1907        } else {
1908            return false;
1909        }
1910    }
1911    
1912    bool MeshAdapter::operator!=(const escript::AbstractDomain& other) const
1913    {
1914        return !(operator==(other));
1915  }  }
1916    
1917  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
1918  {  {
1919     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
1920     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);      return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
1921     checkPasoError();                 package, symmetry, mesh->MPIInfo);
    return out;  
1922  }  }
1923    
1924  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
1925  {  {
1926     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
1927     int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);      return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
1928     checkPasoError();                  package, symmetry, mesh->MPIInfo);
    return out;  
1929  }  }
1930    
1931  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
1932  {  {
1933     return continuousFunction(asAbstractContinuousDomain()).getX();      return continuousFunction(*this).getX();
1934  }  }
1935    
1936  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
1937  {  {
1938     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();      return functionOnBoundary(*this).getNormal();
1939  }  }
1940    
1941  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
1942  {  {
1943     return escript::function(asAbstractContinuousDomain()).getSize();      return escript::function(*this).getSize();
1944  }  }
1945    
1946  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1947  {  {
1948     int *out = NULL;      int *out = NULL;
1949     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
1950     switch (functionSpaceType) {      switch (functionSpaceType) {
1951     case(Nodes):          case Nodes:
1952     out=mesh->Nodes->Id;              out=mesh->Nodes->Id;
1953     break;              break;
1954     case(ReducedNodes):          case ReducedNodes:
1955     out=mesh->Nodes->reducedNodesId;              out=mesh->Nodes->reducedNodesId;
1956     break;              break;
1957     case(Elements):          case Elements:
1958     out=mesh->Elements->Id;          case ReducedElements:
1959     break;              out=mesh->Elements->Id;
1960     case(ReducedElements):              break;
1961     out=mesh->Elements->Id;          case FaceElements:
1962     break;          case ReducedFaceElements:
1963     case(FaceElements):              out=mesh->FaceElements->Id;
1964     out=mesh->FaceElements->Id;              break;
1965     break;          case Points:
1966     case(ReducedFaceElements):              out=mesh->Points->Id;
1967     out=mesh->FaceElements->Id;              break;
1968     break;          case ContactElementsZero:
1969     case(Points):          case ReducedContactElementsZero:
1970     out=mesh->Points->Id;          case ContactElementsOne:
1971     break;          case ReducedContactElementsOne:
1972     case(ContactElementsZero):              out=mesh->ContactElements->Id;
1973     out=mesh->ContactElements->Id;              break;
1974     break;          case DegreesOfFreedom:
1975     case(ReducedContactElementsZero):              out=mesh->Nodes->degreesOfFreedomId;
1976     out=mesh->ContactElements->Id;              break;
1977     break;          case ReducedDegreesOfFreedom:
1978     case(ContactElementsOne):              out=mesh->Nodes->reducedDegreesOfFreedomId;
1979     out=mesh->ContactElements->Id;              break;
1980     break;          default:
1981     case(ReducedContactElementsOne):              stringstream temp;
1982     out=mesh->ContactElements->Id;              temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1983     break;              throw FinleyAdapterException(temp.str());
1984     case(DegreesOfFreedom):              break;
1985     out=mesh->Nodes->degreesOfFreedomId;      }
1986     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;  
1987  }  }
1988  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1989  {  {
1990     int out=0;      int out=0;
1991     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
1992     switch (functionSpaceType) {      switch (functionSpaceType) {
1993     case(Nodes):          case Nodes:
1994     out=mesh->Nodes->Tag[sampleNo];              out=mesh->Nodes->Tag[sampleNo];
1995     break;              break;
1996     case(ReducedNodes):          case ReducedNodes:
1997     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");              throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1998     break;              break;
1999     case(Elements):          case Elements:
2000     out=mesh->Elements->Tag[sampleNo];          case ReducedElements:
2001     break;              out=mesh->Elements->Tag[sampleNo];
2002     case(ReducedElements):              break;
2003     out=mesh->Elements->Tag[sampleNo];          case FaceElements:
2004     break;          case ReducedFaceElements:
2005     case(FaceElements):              out=mesh->FaceElements->Tag[sampleNo];
2006     out=mesh->FaceElements->Tag[sampleNo];              break;
2007     break;          case Points:
2008     case(ReducedFaceElements):              out=mesh->Points->Tag[sampleNo];
2009     out=mesh->FaceElements->Tag[sampleNo];              break;
2010     break;          case ContactElementsZero:
2011     case(Points):          case ReducedContactElementsZero:
2012     out=mesh->Points->Tag[sampleNo];          case ContactElementsOne:
2013     break;          case ReducedContactElementsOne:
2014     case(ContactElementsZero):              out=mesh->ContactElements->Tag[sampleNo];
2015     out=mesh->ContactElements->Tag[sampleNo];              break;
2016     break;          case DegreesOfFreedom:
2017     case(ReducedContactElementsZero):              throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2018     out=mesh->ContactElements->Tag[sampleNo];              break;
2019     break;          case ReducedDegreesOfFreedom:
2020     case(ContactElementsOne):              throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2021     out=mesh->ContactElements->Tag[sampleNo];              break;
2022     break;          default:
2023     case(ReducedContactElementsOne):              stringstream temp;
2024     out=mesh->ContactElements->Tag[sampleNo];              temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2025     break;              throw FinleyAdapterException(temp.str());
2026     case(DegreesOfFreedom):              break;
2027     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");      }
2028     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;  
2029  }  }
2030    
2031    
2032  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
2033  {  {
2034     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2035     escriptDataC tmp=mask.getDataC();      escriptDataC tmp=mask.getDataC();
2036     switch(functionSpaceType) {      switch(functionSpaceType) {
2037     case(Nodes):          case Nodes:
2038     Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);              Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2039     break;              break;
2040     case(ReducedNodes):          case ReducedNodes:
2041     throw FinleyAdapterException("Error - ReducedNodes does not support tags");              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2042     break;              break;
2043     case(DegreesOfFreedom):          case DegreesOfFreedom:
2044     throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2045     break;              break;
2046     case(ReducedDegreesOfFreedom):          case ReducedDegreesOfFreedom:
2047     throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2048     break;              break;
2049     case(Elements):          case Elements:
2050     Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);          case ReducedElements:
2051     break;              Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2052     case(ReducedElements):              break;
2053     Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);          case FaceElements:
2054     break;          case ReducedFaceElements:
2055     case(FaceElements):              Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2056     Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);              break;
2057     break;          case Points:
2058     case(ReducedFaceElements):              Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2059     Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);              break;
2060     break;          case ContactElementsZero:
2061     case(Points):          case ReducedContactElementsZero:
2062     Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);          case ContactElementsOne:
2063     break;          case ReducedContactElementsOne:
2064     case(ContactElementsZero):              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2065     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);              break;
2066     break;          default:
2067     case(ReducedContactElementsZero):              stringstream temp;
2068     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);              temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2069     break;              throw FinleyAdapterException(temp.str());
2070     case(ContactElementsOne):      }
2071     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);      checkFinleyError();
2072     break;  }
2073     case(ReducedContactElementsOne):  
2074     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);  void MeshAdapter::setTagMap(const string& name, int tag)
2075     break;  {
2076     default:      Finley_Mesh* mesh=m_finleyMesh.get();
2077        stringstream temp;      Finley_Mesh_addTagMap(mesh, name.c_str(), tag);
2078        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.");  
2079  }  }
2080    
2081  int MeshAdapter::getTag(const string& name) const  int MeshAdapter::getTag(const string& name) const
2082  {  {
2083     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2084     int tag=0;      int tag = Finley_Mesh_getTag(mesh, name.c_str());
2085     tag=Finley_Mesh_getTag(mesh, name.c_str());      checkFinleyError();
2086     checkPasoError();      return tag;
    // throwStandardException("MeshAdapter::getTag is not implemented.");  
    return tag;  
2087  }  }
2088    
2089  bool MeshAdapter::isValidTagName(const string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2090  {  {
2091     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2092     return Finley_Mesh_isValidTagName(mesh,name.c_str());      return Finley_Mesh_isValidTagName(mesh, name.c_str());
2093  }  }
2094    
2095  string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2096  {  {
2097     stringstream temp;      stringstream temp;
2098     Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2099     Finley_TagMap* tag_map=mesh->TagMap;      Finley_TagMap* tag_map=mesh->TagMap;
2100     while (tag_map) {      while (tag_map) {
2101        temp << tag_map->name;          temp << tag_map->name;
2102        tag_map=tag_map->next;          tag_map=tag_map->next;
2103        if (tag_map) temp << ", ";          if (tag_map) temp << ", ";
2104     }      }
2105     return temp.str();      return temp.str();
2106  }  }
2107    
2108  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2109  {  {
2110    Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2111    dim_t numTags=0;      dim_t numTags=0;
2112    switch(functionSpaceCode) {      switch(functionSpaceCode) {
2113     case(Nodes):          case Nodes:
2114            numTags=mesh->Nodes->numTagsInUse;              numTags=mesh->Nodes->numTagsInUse;
2115            break;              break;
2116     case(ReducedNodes):          case ReducedNodes:
2117            throw FinleyAdapterException("Error - ReducedNodes does not support tags");              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2118            break;              break;
2119     case(DegreesOfFreedom):          case DegreesOfFreedom:
2120            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2121            break;              break;
2122     case(ReducedDegreesOfFreedom):          case ReducedDegreesOfFreedom:
2123            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2124            break;              break;
2125     case(Elements):          case Elements:
2126     case(ReducedElements):          case ReducedElements:
2127            numTags=mesh->Elements->numTagsInUse;              numTags=mesh->Elements->numTagsInUse;
2128            break;              break;
2129     case(FaceElements):          case FaceElements:
2130     case(ReducedFaceElements):          case ReducedFaceElements:
2131            numTags=mesh->FaceElements->numTagsInUse;              numTags=mesh->FaceElements->numTagsInUse;
2132            break;              break;
2133     case(Points):          case Points:
2134            numTags=mesh->Points->numTagsInUse;              numTags=mesh->Points->numTagsInUse;
2135            break;              break;
2136     case(ContactElementsZero):          case ContactElementsZero:
2137     case(ReducedContactElementsZero):          case ReducedContactElementsZero:
2138     case(ContactElementsOne):          case ContactElementsOne:
2139     case(ReducedContactElementsOne):          case ReducedContactElementsOne:
2140            numTags=mesh->ContactElements->numTagsInUse;              numTags=mesh->ContactElements->numTagsInUse;
2141            break;              break;
2142     default:          default:
2143        stringstream temp;              stringstream temp;
2144        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;              temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2145        throw FinleyAdapterException(temp.str());              throw FinleyAdapterException(temp.str());
2146    }      }
2147    return numTags;      return numTags;
2148  }  }
2149    
2150  const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2151  {  {
2152    Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2153    index_t* tags=NULL;      index_t* tags=NULL;
2154    switch(functionSpaceCode) {      switch(functionSpaceCode) {
2155     case(Nodes):          case Nodes:
2156            tags=mesh->Nodes->tagsInUse;              tags=mesh->Nodes->tagsInUse;
2157            break;              break;
2158     case(ReducedNodes):          case ReducedNodes:
2159            throw FinleyAdapterException("Error - ReducedNodes does not support tags");              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2160            break;              break;
2161     case(DegreesOfFreedom):          case DegreesOfFreedom:
2162            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2163            break;              break;
2164     case(ReducedDegreesOfFreedom):          case ReducedDegreesOfFreedom:
2165            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2166            break;              break;
2167     case(Elements):          case Elements:
2168     case(ReducedElements):          case ReducedElements:
2169            tags=mesh->Elements->tagsInUse;              tags=mesh->Elements->tagsInUse;
2170            break;              break;
2171     case(FaceElements):          case FaceElements:
2172     case(ReducedFaceElements):          case ReducedFaceElements:
2173            tags=mesh->FaceElements->tagsInUse;              tags=mesh->FaceElements->tagsInUse;
2174            break;              break;
2175     case(Points):          case Points:
2176            tags=mesh->Points->tagsInUse;              tags=mesh->Points->tagsInUse;
2177            break;              break;
2178     case(ContactElementsZero):          case ContactElementsZero:
2179     case(ReducedContactElementsZero):          case ReducedContactElementsZero:
2180     case(ContactElementsOne):          case ContactElementsOne:
2181     case(ReducedContactElementsOne):          case ReducedContactElementsOne:
2182            tags=mesh->ContactElements->tagsInUse;              tags=mesh->ContactElements->tagsInUse;
2183            break;              break;
2184     default:          default:
2185        stringstream temp;              stringstream temp;
2186        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;              temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2187        throw FinleyAdapterException(temp.str());              throw FinleyAdapterException(temp.str());
2188    }      }
2189    return tags;      return tags;
2190  }  }
2191    
2192    
2193  bool MeshAdapter::canTag(int functionSpaceCode) const  bool MeshAdapter::canTag(int functionSpaceCode) const
2194  {  {
2195    switch(functionSpaceCode) {      switch(functionSpaceCode) {
2196     case(Nodes):          case Nodes:
2197     case(Elements):          case Elements:
2198     case(ReducedElements):          case ReducedElements:
2199     case(FaceElements):          case FaceElements:
2200     case(ReducedFaceElements):          case ReducedFaceElements:
2201     case(Points):          case Points:
2202     case(ContactElementsZero):          case ContactElementsZero:
2203     case(ReducedContactElementsZero):          case ReducedContactElementsZero:
2204     case(ContactElementsOne):          case ContactElementsOne:
2205     case(ReducedContactElementsOne):          case ReducedContactElementsOne:
2206            return true;              return true;
2207     case(ReducedNodes):          default:
2208     case(DegreesOfFreedom):              return false;
2209     case(ReducedDegreesOfFreedom):      }
       return false;  
    default:  
     return false;  
   }  
2210  }  }
2211    
2212  AbstractDomain::StatusType MeshAdapter::getStatus() const  escript::AbstractDomain::StatusType MeshAdapter::getStatus() const
2213  {  {
2214    Finley_Mesh* mesh=m_finleyMesh.get();      Finley_Mesh* mesh=m_finleyMesh.get();
2215    return Finley_Mesh_getStatus(mesh);      return Finley_Mesh_getStatus(mesh);
2216  }  }
2217    
2218  int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const  int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2219  {  {
2220          Finley_Mesh* mesh=m_finleyMesh.get();
2221    Finley_Mesh* mesh=m_finleyMesh.get();      int order =-1;
2222    int order =-1;      switch(functionSpaceCode) {
2223    switch(functionSpaceCode) {          case Nodes:
2224     case(Nodes):          case DegreesOfFreedom:
2225     case(DegreesOfFreedom):              order=mesh->approximationOrder;
2226            order=mesh->approximationOrder;              break;
2227            break;          case ReducedNodes:
2228     case(ReducedNodes):          case ReducedDegreesOfFreedom:
2229     case(ReducedDegreesOfFreedom):              order=mesh->reducedApproximationOrder;
2230            order=mesh->reducedApproximationOrder;              break;
2231            break;          case Elements:
2232     case(Elements):          case FaceElements:
2233     case(FaceElements):          case Points:
2234     case(Points):          case ContactElementsZero:
2235     case(ContactElementsZero):          case ContactElementsOne:
2236     case(ContactElementsOne):              order=mesh->integrationOrder;
2237            order=mesh->integrationOrder;              break;
2238            break;          case ReducedElements:
2239     case(ReducedElements):          case ReducedFaceElements:
2240     case(ReducedFaceElements):          case ReducedContactElementsZero:
2241     case(ReducedContactElementsZero):          case ReducedContactElementsOne:
2242     case(ReducedContactElementsOne):              order=mesh->reducedIntegrationOrder;
2243            order=mesh->reducedIntegrationOrder;              break;
2244            break;          default:
2245     default:              stringstream temp;
2246        stringstream temp;              temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2247        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;              throw FinleyAdapterException(temp.str());
2248        throw FinleyAdapterException(temp.str());      }
2249    }      return order;
2250    return order;  }
2251    
2252    bool MeshAdapter::supportsContactElements() const
2253    {
2254        return true;
2255  }  }
2256    
2257  ReferenceElementSetWrapper::ReferenceElementSetWrapper(ElementTypeId id, index_t order, index_t reducedOrder)  ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id,
2258            index_t order, index_t reducedOrder)
2259  {  {
2260    m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2261  }  }
2262    
2263  ReferenceElementSetWrapper::~ReferenceElementSetWrapper()  ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2264  {  {
2265    Finley_ReferenceElementSet_dealloc(m_refSet);      Finley_ReferenceElementSet_dealloc(m_refSet);
2266    }
2267    
2268    void MeshAdapter::addDiracPoints(const vector<double>& points,
2269                                     const vector<int>& tags) const
2270    {
2271        // points will be flattened
2272        const int dim = getDim();
2273        int numPoints=points.size()/dim;
2274        int numTags=tags.size();
2275        Finley_Mesh* mesh=m_finleyMesh.get();
2276    
2277        if ( points.size() % dim != 0 ) {
2278            throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2279        }
2280    
2281        if ((numTags > 0) && (numPoints != numTags))
2282            throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2283    
2284        Finley_Mesh_addPoints(mesh, numPoints, &points[0], &tags[0]);
2285        checkFinleyError();
2286  }  }
2287    
2288    // void MeshAdapter::addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2289    // {
2290    //       const int dim = getDim();
2291    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2292    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2293    //       Finley_Mesh* mesh=m_finleyMesh.get();
2294    //
2295    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2296    //       throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2297    //
2298    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2299    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2300    //
2301    //       for (int i=0;i<numPoints;++i) {
2302    //         int tag_id=-1;
2303    //         int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2304    //         if  ( numComps !=   dim ) {
2305    //                stringstream temp;
2306    //                temp << "Error - illegal number of components " << numComps << " for point " << i;
2307    //                throw FinleyAdapterException(temp.str());
2308    //         }
2309    //         points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2310    //         if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2311    //         if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2312    //
2313    //         if ( numTags > 0) {
2314    //                boost::python::extract<string> ex_str(tags[i]);
2315    //                if  ( ex_str.check() ) {
2316    //                    tag_id=getTag( ex_str());
2317    //                } else {
2318    //                     boost::python::extract<int> ex_int(tags[i]);
2319    //                     if ( ex_int.check() ) {
2320    //                         tag_id=ex_int();
2321    //                     } else {
2322    //                          stringstream temp;
2323    //                          temp << "Error - unable to extract tag for point " << i;
2324    //                          throw FinleyAdapterException(temp.str());
2325    //                    }
2326    //                }
2327    //         }
2328    //            tags_ptr[i]=tag_id;
2329    //       }
2330    //
2331    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2332    //       checkPasoError();
2333    //
2334    //       TMPMEMFREE(points_ptr);
2335    //       TMPMEMFREE(tags_ptr);
2336    // }
2337    
2338    /*
2339    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2340    {
2341        boost::python::list points =  boost::python::list();
2342        boost::python::list tags =  boost::python::list();
2343        points.append(point);
2344        tags.append(tag);
2345        addDiracPoints(points, tags);
2346    }
2347    */
2348    
2349    /*
2350    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const string& tag) const
2351    {
2352        boost::python::list points =   boost::python::list();
2353        boost::python::list tags =   boost::python::list();
2354        points.append(point);
2355        tags.append(tag);
2356        addDiracPoints(points, tags);
2357    }
2358    */
2359  }  // end of namespace  }  // end of namespace
2360    

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

  ViewVC Help
Powered by ViewVC 1.1.26