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

Legend:
Removed from v.2940  
changed lines
  Added in v.4492

  ViewVC Help
Powered by ViewVC 1.1.26