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

Legend:
Removed from v.3317  
changed lines
  Added in v.4800

  ViewVC Help
Powered by ViewVC 1.1.26