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

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

  ViewVC Help
Powered by ViewVC 1.1.26