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

Legend:
Removed from v.3911  
changed lines
  Added in v.4612

  ViewVC Help
Powered by ViewVC 1.1.26