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

Diff of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

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

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

Legend:
Removed from v.4346  
changed lines
  Added in v.4437

  ViewVC Help
Powered by ViewVC 1.1.26