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

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

  ViewVC Help
Powered by ViewVC 1.1.26