/[escript]/branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Diff of /branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp

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

trunk/finley/src/CPPAdapter/MeshAdapter.cpp revision 2385 by gross, Wed Apr 15 03:52:38 2009 UTC branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp revision 3239 by jfenwick, Tue Oct 5 05:31:20 2010 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2010 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 18  Line 18 
18  #ifdef USE_NETCDF  #ifdef USE_NETCDF
19  #include <netcdfcpp.h>  #include <netcdfcpp.h>
20  #endif  #endif
21  #ifdef PASO_MPI  #ifdef ESYS_MPI
22  #include <mpi.h>  #include "esysUtils/Esys_MPI.h"
 #include "paso/Paso_MPI.h"  
23  #endif  #endif
24  extern "C" {  extern "C" {
25  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
# Line 29  extern "C" { Line 28  extern "C" {
28  using namespace std;  using namespace std;
29  using namespace escript;  using namespace escript;
30    
31  namespace finley {  namespace dudley {
32    
33  //  //
34  // define the static constants  // define the static constants
35  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
36  const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;  const int MeshAdapter::DegreesOfFreedom=DUDLEY_DEGREES_OF_FREEDOM;
37  const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;  const int MeshAdapter::ReducedDegreesOfFreedom=DUDLEY_REDUCED_DEGREES_OF_FREEDOM;
38  const int MeshAdapter::Nodes=FINLEY_NODES;  const int MeshAdapter::Nodes=DUDLEY_NODES;
39  const int MeshAdapter::ReducedNodes=FINLEY_REDUCED_NODES;  const int MeshAdapter::ReducedNodes=DUDLEY_REDUCED_NODES;
40  const int MeshAdapter::Elements=FINLEY_ELEMENTS;  const int MeshAdapter::Elements=DUDLEY_ELEMENTS;
41  const int MeshAdapter::ReducedElements=FINLEY_REDUCED_ELEMENTS;  const int MeshAdapter::ReducedElements=DUDLEY_REDUCED_ELEMENTS;
42  const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;  const int MeshAdapter::FaceElements=DUDLEY_FACE_ELEMENTS;
43  const int MeshAdapter::ReducedFaceElements=FINLEY_REDUCED_FACE_ELEMENTS;  const int MeshAdapter::ReducedFaceElements=DUDLEY_REDUCED_FACE_ELEMENTS;
44  const int MeshAdapter::Points=FINLEY_POINTS;  const int MeshAdapter::Points=DUDLEY_POINTS;
 const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;  
 const int MeshAdapter::ReducedContactElementsZero=FINLEY_REDUCED_CONTACT_ELEMENTS_1;  
 const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;  
 const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;  
45    
46  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)  MeshAdapter::MeshAdapter(Dudley_Mesh* dudleyMesh)
47  {  {
48     setFunctionSpaceTypeNames();     setFunctionSpaceTypeNames();
49     //     //
50     // need to use a null_deleter as Finley_Mesh_free deletes the pointer     // need to use a null_deleter as Dudley_Mesh_free deletes the pointer
51     // for us.     // for us.
52     m_finleyMesh.reset(finleyMesh,null_deleter());     m_dudleyMesh.reset(dudleyMesh,null_deleter());
53  }  }
54    
55  //  //
56  // The copy constructor should just increment the use count  // The copy constructor should just increment the use count
57  MeshAdapter::MeshAdapter(const MeshAdapter& in):  MeshAdapter::MeshAdapter(const MeshAdapter& in):
58  m_finleyMesh(in.m_finleyMesh)  m_dudleyMesh(in.m_dudleyMesh)
59  {  {
60     setFunctionSpaceTypeNames();     setFunctionSpaceTypeNames();
61  }  }
# Line 70  MeshAdapter::~MeshAdapter() Line 65  MeshAdapter::~MeshAdapter()
65     //     //
66     // I hope the case for the pointer being zero has been taken care of.     // I hope the case for the pointer being zero has been taken care of.
67     //  cout << "In MeshAdapter destructor." << endl;     //  cout << "In MeshAdapter destructor." << endl;
68     if (m_finleyMesh.unique()) {     if (m_dudleyMesh.unique()) {
69        Finley_Mesh_free(m_finleyMesh.get());        Dudley_Mesh_free(m_dudleyMesh.get());
70     }     }
71  }  }
72    
73  int MeshAdapter::getMPISize() const  int MeshAdapter::getMPISize() const
74  {  {
75     return m_finleyMesh.get()->MPIInfo->size;     return m_dudleyMesh.get()->MPIInfo->size;
76  }  }
77  int MeshAdapter::getMPIRank() const  int MeshAdapter::getMPIRank() const
78  {  {
79     return m_finleyMesh.get()->MPIInfo->rank;     return m_dudleyMesh.get()->MPIInfo->rank;
80  }  }
81  void MeshAdapter::MPIBarrier() const  void MeshAdapter::MPIBarrier() const
82  {  {
83  #ifdef PASO_MPI  #ifdef ESYS_MPI
84     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);     MPI_Barrier(m_dudleyMesh.get()->MPIInfo->comm);
85  #endif  #endif
86     return;     return;
87  }  }
88  bool MeshAdapter::onMasterProcessor() const  bool MeshAdapter::onMasterProcessor() const
89  {  {
90     return m_finleyMesh.get()->MPIInfo->rank == 0;     return m_dudleyMesh.get()->MPIInfo->rank == 0;
91  }  }
92    
93    
94  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  #ifdef ESYS_MPI
95     return m_finleyMesh.get();    MPI_Comm
96    #else
97      unsigned int
98    #endif
99    MeshAdapter::getMPIComm() const
100    {
101    #ifdef ESYS_MPI
102        return m_dudleyMesh->MPIInfo->comm;
103    #else
104        return 0;
105    #endif
106  }  }
107    
108  void MeshAdapter::write(const std::string& fileName) const  
109    Dudley_Mesh* MeshAdapter::getDudley_Mesh() const {
110       return m_dudleyMesh.get();
111    }
112    
113    void MeshAdapter::write(const string& fileName) const
114  {  {
115     char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;     char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
116     strcpy(fName,fileName.c_str());     strcpy(fName,fileName.c_str());
117     Finley_Mesh_write(m_finleyMesh.get(),fName);     Dudley_Mesh_write(m_dudleyMesh.get(),fName);
118     checkFinleyError();     checkDudleyError();
119     TMPMEMFREE(fName);     TMPMEMFREE(fName);
120  }  }
121    
122  void MeshAdapter::Print_Mesh_Info(const bool full=false) const  void MeshAdapter::Print_Mesh_Info(const bool full) const
123  {  {
124     Finley_PrintMesh_Info(m_finleyMesh.get(), full);     Dudley_PrintMesh_Info(m_dudleyMesh.get(), full);
125  }  }
126    
127  void MeshAdapter::dump(const std::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();     Dudley_Mesh *mesh = m_dudleyMesh.get();
134     Finley_TagMap* tag_map;     Dudley_TagMap* tag_map;
135     int num_Tags = 0;     int num_Tags = 0;
136     int mpi_size             = mesh->MPIInfo->size;     int mpi_size             = mesh->MPIInfo->size;
137     int mpi_rank             = mesh->MPIInfo->rank;     int mpi_rank             = mesh->MPIInfo->rank;
# Line 129  void MeshAdapter::dump(const std::string Line 139  void MeshAdapter::dump(const std::string
139     int numNodes             = mesh->Nodes->numNodes;     int numNodes             = mesh->Nodes->numNodes;
140     int num_Elements         = mesh->Elements->numElements;     int num_Elements         = mesh->Elements->numElements;
141     int num_FaceElements         = mesh->FaceElements->numElements;     int num_FaceElements         = mesh->FaceElements->numElements;
    int num_ContactElements      = mesh->ContactElements->numElements;  
142     int num_Points           = mesh->Points->numElements;     int num_Points           = mesh->Points->numElements;
143     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
144     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
145     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;  #ifdef ESYS_MPI
 #ifdef PASO_MPI  
146     MPI_Status status;     MPI_Status status;
147  #endif  #endif
148    
149  /* Incoming token indicates it's my turn to write */  /* Incoming token indicates it's my turn to write */
150  #ifdef PASO_MPI  #ifdef ESYS_MPI
151     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
152  #endif  #endif
153    
154     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
155                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
156    
157     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
158     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
159     num_Tags = 0;     num_Tags = 0;
160     if (tag_map) {     while (tag_map) {
161        while (tag_map) {        num_Tags++;
162           num_Tags++;        tag_map=tag_map->next;
          tag_map=tag_map->next;  
       }  
163     }     }
164    
165     // NetCDF error handler     // NetCDF error handler
166     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
167     // Create the file.     // Create the file.
168     NcFile dataFile(newFileName, NcFile::Replace);     NcFile dataFile(newFileName, NcFile::Replace);
169       string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
170     // check if writing was successful     // check if writing was successful
171     if (!dataFile.is_valid())     if (!dataFile.is_valid())
172        throw DataException("Error - MeshAdapter::dump: opening of NetCDF file for output failed: " + *newFileName);        throw DataException(msgPrefix+"Open file for output");
173    
174     // Define dimensions (num_Elements and dim_Elements are identical, dim_Elements only appears if > 0)     // Define dimensions (num_Elements and dim_Elements are identical,
175       // dim_Elements only appears if > 0)
176     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
177        throw DataException("Error - MeshAdapter::dump: appending dimension numNodes to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numNodes)");
178     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
179        throw DataException("Error - MeshAdapter::dump: appending dimension numDim to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numDim)");
180     if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )     if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )
181        throw DataException("Error - MeshAdapter::dump: appending dimension mpi_size to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(mpi_size)");
182     if (num_Elements>0)     if (num_Elements>0)
183        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
184           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements)");
185     if (num_FaceElements>0)     if (num_FaceElements>0)
186        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
187           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
    if (num_ContactElements>0)  
       if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )  
          throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements to netCDF file failed: " + *newFileName);  
188     if (num_Points>0)     if (num_Points>0)
189        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
190           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Points to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Points)");
191     if (num_Elements>0)     if (num_Elements>0)
192        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
193           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
194     if (num_FaceElements>0)     if (num_FaceElements>0)
195        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
196           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements_numNodes to netCDF file failed: " + *newFileName);           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("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);  
197     if (num_Tags>0)     if (num_Tags>0)
198        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
199           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Tags)");
200    
201     // Attributes: MPI size, MPI rank, Name, order, reduced_order     // Attributes: MPI size, MPI rank, Name, order, reduced_order
202     if (!dataFile.add_att("mpi_size", mpi_size) )     if (!dataFile.add_att("mpi_size", mpi_size) )
203        throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_size)");
204     if (!dataFile.add_att("mpi_rank", mpi_rank) )     if (!dataFile.add_att("mpi_rank", mpi_rank) )
205        throw DataException("Error - MeshAdapter::dump: appending mpi_rank to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_rank)");
206     if (!dataFile.add_att("Name",mesh->Name) )     if (!dataFile.add_att("Name",mesh->Name) )
207        throw DataException("Error - MeshAdapter::dump: appending Name to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Name)");
208     if (!dataFile.add_att("numDim",numDim) )     if (!dataFile.add_att("numDim",numDim) )
209        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
210     if (!dataFile.add_att("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
211        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
212     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
213        throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(reduced_order)");
214     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
215        throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(numNodes)");
216     if (!dataFile.add_att("num_Elements",num_Elements) )     if (!dataFile.add_att("num_Elements",num_Elements) )
217        throw DataException("Error - MeshAdapter::dump: appending num_Elements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements)");
218     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
219        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements)");
    if (!dataFile.add_att("num_ContactElements",num_ContactElements) )  
       throw DataException("Error - MeshAdapter::dump: appending num_ContactElements to NetCDF file failed: " + *newFileName);  
220     if (!dataFile.add_att("num_Points",num_Points) )     if (!dataFile.add_att("num_Points",num_Points) )
221        throw DataException("Error - MeshAdapter::dump: appending num_Points to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Points)");
222     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
223        throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
224     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
225        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
226     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->etype) )
227        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Elements_TypeId)");
228     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->etype) )
229        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
230     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Points_TypeId", mesh->Points->etype) )
231        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Points_TypeId)");
    if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )  
       throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);  
    if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )  
       throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);  
232     if (!dataFile.add_att("num_Tags", num_Tags) )     if (!dataFile.add_att("num_Tags", num_Tags) )
233        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Tags)");
234    
235     // // // // // Nodes // // // // //     // // // // // Nodes // // // // //
236    
237     // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)     // Nodes nodeDistribution
238       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
239          throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
240       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
241       if (! (ids->put(int_ptr, mpi_size+1)) )
242          throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
243    
244       // Nodes degreesOfFreedomDistribution
245       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
246          throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
247       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
248       if (! (ids->put(int_ptr, mpi_size+1)) )
249          throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
250    
251       // Only write nodes if non-empty because NetCDF doesn't like empty arrays
252       // (it treats them as NC_UNLIMITED)
253     if (numNodes>0) {     if (numNodes>0) {
254    
255        // Nodes Id        // Nodes Id
256        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
257           throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Id)");
258        int_ptr = &mesh->Nodes->Id[0];        int_ptr = &mesh->Nodes->Id[0];
259        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
260           throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Id)");
261    
262        // Nodes Tag        // Nodes Tag
263        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
264           throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Tag)");
265        int_ptr = &mesh->Nodes->Tag[0];        int_ptr = &mesh->Nodes->Tag[0];
266        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
267           throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Tag)");
268    
269        // Nodes gDOF        // Nodes gDOF
270        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
271           throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
272        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
273        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
274           throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gDOF)");
275    
276        // Nodes global node index        // Nodes global node index
277        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
278           throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gNI)");
279        int_ptr = &mesh->Nodes->globalNodesIndex[0];        int_ptr = &mesh->Nodes->globalNodesIndex[0];
280        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
281           throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gNI)");
282    
283        // Nodes grDof        // Nodes grDof
284        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
285           throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
286        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
287        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
288           throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grDfI)");
289    
290        // Nodes grNI        // Nodes grNI
291        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
292           throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grNI)");
293        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
294        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
295           throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grNI)");
296    
297        // Nodes Coordinates        // Nodes Coordinates
298        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
299           throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
300        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
301           throw DataException("Error - MeshAdapter::dump: copy Nodes_Coordinates to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Coordinates)");
   
       // Nodes degreesOfFreedomDistribution  
       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )  
          throw DataException("Error - MeshAdapter::dump: appending Nodes_DofDistribution to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];  
       if (! (ids->put(int_ptr, mpi_size+1)) )  
          throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);  
   
       // Nodes nodeDistribution  
       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )  
          throw DataException("Error - MeshAdapter::dump: appending Nodes_NodeDistribution to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];  
       if (! (ids->put(int_ptr, mpi_size+1)) )  
          throw DataException("Error - MeshAdapter::dump: copy Nodes_NodeDistribution to netCDF buffer failed: " + *newFileName);  
302    
303     }     }
304    
# Line 311  void MeshAdapter::dump(const std::string Line 308  void MeshAdapter::dump(const std::string
308    
309        // Elements_Id        // Elements_Id
310        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
311           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Id)");
312        int_ptr = &mesh->Elements->Id[0];        int_ptr = &mesh->Elements->Id[0];
313        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
314           throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Id)");
315    
316        // Elements_Tag        // Elements_Tag
317        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
318           throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Tag)");
319        int_ptr = &mesh->Elements->Tag[0];        int_ptr = &mesh->Elements->Tag[0];
320        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
321           throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Tag)");
322    
323        // Elements_Owner        // Elements_Owner
324        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
325           throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Owner)");
326        int_ptr = &mesh->Elements->Owner[0];        int_ptr = &mesh->Elements->Owner[0];
327        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
328           throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Owner)");
329    
330        // Elements_Color        // Elements_Color
331        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
332           throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Color)");
333        int_ptr = &mesh->Elements->Color[0];        int_ptr = &mesh->Elements->Color[0];
334        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
335           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Color)");
336    
337        // Elements_Nodes        // Elements_Nodes
338        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
339           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Nodes)");
340        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
341           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Nodes)");
342    
343     }     }
344    
# Line 351  void MeshAdapter::dump(const std::string Line 348  void MeshAdapter::dump(const std::string
348    
349        // FaceElements_Id        // FaceElements_Id
350        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
351           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Id)");
352        int_ptr = &mesh->FaceElements->Id[0];        int_ptr = &mesh->FaceElements->Id[0];
353        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
354           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Id)");
355    
356        // FaceElements_Tag        // FaceElements_Tag
357        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
358           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
359        int_ptr = &mesh->FaceElements->Tag[0];        int_ptr = &mesh->FaceElements->Tag[0];
360        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
361           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Tag)");
362    
363        // FaceElements_Owner        // FaceElements_Owner
364        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
365           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
366        int_ptr = &mesh->FaceElements->Owner[0];        int_ptr = &mesh->FaceElements->Owner[0];
367        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
368           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Owner)");
369    
370        // FaceElements_Color        // FaceElements_Color
371        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
372           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Color)");
373        int_ptr = &mesh->FaceElements->Color[0];        int_ptr = &mesh->FaceElements->Color[0];
374        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
375           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Color)");
376    
377        // FaceElements_Nodes        // FaceElements_Nodes
378        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
379           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
380        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
381           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);           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("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->ContactElements->Id[0];  
       if (! (ids->put(int_ptr, num_ContactElements)) )  
          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);  
   
       // ContactElements_Tag  
       if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )  
          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->ContactElements->Tag[0];  
       if (! (ids->put(int_ptr, num_ContactElements)) )  
          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);  
   
       // ContactElements_Owner  
       if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )  
          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->ContactElements->Owner[0];  
       if (! (ids->put(int_ptr, num_ContactElements)) )  
          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);  
   
       // ContactElements_Color  
       if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )  
          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->ContactElements->Color[0];  
       if (! (ids->put(int_ptr, num_ContactElements)) )  
          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);  
   
       // ContactElements_Nodes  
       if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )  
          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);  
       if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )  
          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);  
382    
383     }     }
384    
# Line 433  void MeshAdapter::dump(const std::string Line 390  void MeshAdapter::dump(const std::string
390    
391        // Points_Id        // Points_Id
392        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
393           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Id)");
394        int_ptr = &mesh->Points->Id[0];        int_ptr = &mesh->Points->Id[0];
395        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
396           throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Id)");
397    
398        // Points_Tag        // Points_Tag
399        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
400           throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Tag)");
401        int_ptr = &mesh->Points->Tag[0];        int_ptr = &mesh->Points->Tag[0];
402        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
403           throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Tag)");
404    
405        // Points_Owner        // Points_Owner
406        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
407           throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Owner)");
408        int_ptr = &mesh->Points->Owner[0];        int_ptr = &mesh->Points->Owner[0];
409        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
410           throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Owner)");
411    
412        // Points_Color        // Points_Color
413        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
414           throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Color)");
415        int_ptr = &mesh->Points->Color[0];        int_ptr = &mesh->Points->Color[0];
416        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
417           throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Color)");
418    
419        // Points_Nodes        // Points_Nodes
420        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
421        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
422           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Nodes)");
423        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
424           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Nodes)");
425    
426     }     }
427    
# Line 488  void MeshAdapter::dump(const std::string Line 445  void MeshAdapter::dump(const std::string
445    
446        // Tags_keys        // Tags_keys
447        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
448           throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Tags_keys)");
449        int_ptr = &Tags_keys[0];        int_ptr = &Tags_keys[0];
450        if (! (ids->put(int_ptr, num_Tags)) )        if (! (ids->put(int_ptr, num_Tags)) )
451           throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Tags_keys)");
452    
453        // Tags_names_*        // Tags_names_*
454        // This is an array of strings, it should be stored as an array but instead I have hacked in one attribute per string        // This is an array of strings, it should be stored as an array but
455        // because the NetCDF manual doesn't tell how to do an array of strings        // instead I have hacked in one attribute per string because the NetCDF
456          // manual doesn't tell how to do an array of strings
457        tag_map = mesh->TagMap;        tag_map = mesh->TagMap;
458        if (tag_map) {        if (tag_map) {
459           int i = 0;           int i = 0;
460           while (tag_map) {           while (tag_map) {
461              sprintf(name_temp, "Tags_name_%d", i);              sprintf(name_temp, "Tags_name_%d", i);
462              if (!dataFile.add_att(name_temp, tag_map->name) )              if (!dataFile.add_att(name_temp, tag_map->name) )
463                 throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);                 throw DataException(msgPrefix+"add_att(Tags_names_XX)");
464              tag_map=tag_map->next;              tag_map=tag_map->next;
465              i++;              i++;
466           }           }
467        }        }
468    
469        TMPMEMFREE(Tags_keys);        TMPMEMFREE(Tags_keys);
   
470     }     }
471    
472  /* Send token to next MPI process so he can take his turn */  /* Send token to next MPI process so he can take his turn */
473  #ifdef PASO_MPI  #ifdef ESYS_MPI
474     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) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
475  #endif  #endif
476    
477     // NetCDF file is closed by destructor of NcFile object     // NetCDF file is closed by destructor of NcFile object
478    
479  #else  #else
480     Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");     Dudley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
481  #endif  /* USE_NETCDF */  #endif  /* USE_NETCDF */
482     checkFinleyError();     checkDudleyError();
483  }  }
484    
485  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
486  {  {
487     return "FinleyMesh";     return "DudleyMesh";
488  }  }
489    
490  string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const  string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
# Line 551  bool MeshAdapter::isValidFunctionSpaceTy Line 508  bool MeshAdapter::isValidFunctionSpaceTy
508  void MeshAdapter::setFunctionSpaceTypeNames()  void MeshAdapter::setFunctionSpaceTypeNames()
509  {  {
510     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
511     (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));     (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Dudley_DegreesOfFreedom"));
    m_functionSpaceTypeNames.insert  
    (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));  
512     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
513     (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));     (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Dudley_ReducedDegreesOfFreedom"));
514     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
515     (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));     (FunctionSpaceNamesMapType::value_type(Nodes,"Dudley_Nodes"));
516     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
517     (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));     (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Dudley_Reduced_Nodes"));
518     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
519     (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));     (FunctionSpaceNamesMapType::value_type(Elements,"Dudley_Elements"));
520     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
521     (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));     (FunctionSpaceNamesMapType::value_type(ReducedElements,"Dudley_Reduced_Elements"));
522     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
523     (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));     (FunctionSpaceNamesMapType::value_type(FaceElements,"Dudley_Face_Elements"));
524     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
525     (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));     (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Dudley_Reduced_Face_Elements"));
526     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
527     (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));     (FunctionSpaceNamesMapType::value_type(Points,"Dudley_Points"));
    m_functionSpaceTypeNames.insert  
    (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));  
    m_functionSpaceTypeNames.insert  
    (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));  
    m_functionSpaceTypeNames.insert  
    (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));  
528  }  }
529    
530  int MeshAdapter::getContinuousFunctionCode() const  int MeshAdapter::getContinuousFunctionCode() const
# Line 607  int MeshAdapter::getReducedFunctionOnBou Line 556  int MeshAdapter::getReducedFunctionOnBou
556    
557  int MeshAdapter::getFunctionOnContactZeroCode() const  int MeshAdapter::getFunctionOnContactZeroCode() const
558  {  {
559     return ContactElementsZero;     throw DudleyAdapterException("Dudley does not support contact elements.");
560  }  }
561    
562  int MeshAdapter::getReducedFunctionOnContactZeroCode() const  int MeshAdapter::getReducedFunctionOnContactZeroCode() const
563  {  {
564     return ReducedContactElementsZero;     throw DudleyAdapterException("Dudley does not support contact elements.");
565  }  }
566    
567  int MeshAdapter::getFunctionOnContactOneCode() const  int MeshAdapter::getFunctionOnContactOneCode() const
568  {  {
569     return ContactElementsOne;     throw DudleyAdapterException("Dudley does not support contact elements.");
570  }  }
571    
572  int MeshAdapter::getReducedFunctionOnContactOneCode() const  int MeshAdapter::getReducedFunctionOnContactOneCode() const
573  {  {
574     return ReducedContactElementsOne;     throw DudleyAdapterException("Dudley does not support contact elements.");
575  }  }
576    
577  int MeshAdapter::getSolutionCode() const  int MeshAdapter::getSolutionCode() const
# Line 643  int MeshAdapter::getDiracDeltaFunctionCo Line 594  int MeshAdapter::getDiracDeltaFunctionCo
594  //  //
595  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
596  {  {
597     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
598     int numDim=Finley_Mesh_getDim(mesh);     int numDim=Dudley_Mesh_getDim(mesh);
599     checkFinleyError();     checkDudleyError();
600     return numDim;     return numDim;
601  }  }
602    
# Line 654  int MeshAdapter::getDim() const Line 605  int MeshAdapter::getDim() const
605  //  //
606  int MeshAdapter::getNumDataPointsGlobal() const  int MeshAdapter::getNumDataPointsGlobal() const
607  {  {
608     return Finley_NodeFile_getGlobalNumNodes(m_finleyMesh.get()->Nodes);     return Dudley_NodeFile_getGlobalNumNodes(m_dudleyMesh.get()->Nodes);
609  }  }
610    
611  //  //
# Line 665  pair<int,int> MeshAdapter::getDataShape( Line 616  pair<int,int> MeshAdapter::getDataShape(
616  {  {
617     int numDataPointsPerSample=0;     int numDataPointsPerSample=0;
618     int numSamples=0;     int numSamples=0;
619     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
620     switch (functionSpaceCode) {     switch (functionSpaceCode) {
621     case(Nodes):     case(Nodes):
622     numDataPointsPerSample=1;     numDataPointsPerSample=1;
623     numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);     numSamples=Dudley_NodeFile_getNumNodes(mesh->Nodes);
624     break;     break;
625     case(ReducedNodes):     case(ReducedNodes):
626     numDataPointsPerSample=1;     numDataPointsPerSample=1;
627     numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);     numSamples=Dudley_NodeFile_getNumReducedNodes(mesh->Nodes);
628     break;     break;
629     case(Elements):     case(Elements):
630     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
631        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
632        numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->Elements->numLocalDim+1/*referenceElementSet->referenceElement->BasisFunctions->numQuadNodes*/;
633     }     }
634     break;     break;
635     case(ReducedElements):     case(ReducedElements):
636     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
637        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
638        numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=(mesh->Elements->numLocalDim==0)?0:1;
639     }     }
640     break;     break;
641     case(FaceElements):     case(FaceElements):
642     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
643        numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->numLocalDim+1/*referenceElementSet->referenceElement->BasisFunctions->numQuadNodes*/;
644        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
645     }     }
646     break;     break;
647     case(ReducedFaceElements):     case(ReducedFaceElements):
648     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
649        numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=(mesh->FaceElements->numLocalDim==0)?0:1/*referenceElementSet->referenceElementReducedQuadrature->BasisFunctions->numQuadNodes*/;
650        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
651     }     }
652     break;     break;
# Line 705  pair<int,int> MeshAdapter::getDataShape( Line 656  pair<int,int> MeshAdapter::getDataShape(
656        numSamples=mesh->Points->numElements;        numSamples=mesh->Points->numElements;
657     }     }
658     break;     break;
    case(ContactElementsZero):  
    if (mesh->ContactElements!=NULL) {  
       numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;  
       numSamples=mesh->ContactElements->numElements;  
    }  
    break;  
    case(ReducedContactElementsZero):  
    if (mesh->ContactElements!=NULL) {  
       numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;  
       numSamples=mesh->ContactElements->numElements;  
    }  
    break;  
    case(ContactElementsOne):  
    if (mesh->ContactElements!=NULL) {  
       numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;  
       numSamples=mesh->ContactElements->numElements;  
    }  
    break;  
    case(ReducedContactElementsOne):  
    if (mesh->ContactElements!=NULL) {  
       numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;  
       numSamples=mesh->ContactElements->numElements;  
    }  
    break;  
659     case(DegreesOfFreedom):     case(DegreesOfFreedom):
660     if (mesh->Nodes!=NULL) {     if (mesh->Nodes!=NULL) {
661        numDataPointsPerSample=1;        numDataPointsPerSample=1;
662        numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);        numSamples=Dudley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
663     }     }
664     break;     break;
665     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
666     if (mesh->Nodes!=NULL) {     if (mesh->Nodes!=NULL) {
667        numDataPointsPerSample=1;        numDataPointsPerSample=1;
668        numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);        numSamples=Dudley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
669     }     }
670     break;     break;
671     default:     default:
672        stringstream temp;        stringstream temp;
673        temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
674        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
675        break;        break;
676     }     }
677     return pair<int,int>(numDataPointsPerSample,numSamples);     return pair<int,int>(numDataPointsPerSample,numSamples);
# Line 756  pair<int,int> MeshAdapter::getDataShape( Line 683  pair<int,int> MeshAdapter::getDataShape(
683  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
684                                   SystemMatrixAdapter& mat, escript::Data& rhs,                                   SystemMatrixAdapter& mat, escript::Data& rhs,
685                                   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,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
686                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y) const
                                  const escript::Data& d_contact,const escript::Data& y_contact) const  
687  {  {
688     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
689     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
# Line 768  void MeshAdapter::addPDEToSystem( Line 694  void MeshAdapter::addPDEToSystem(
694     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
695     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
696     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
    escriptDataC _d_contact=d_contact.getDataC();  
    escriptDataC _y_contact=y_contact.getDataC();  
697    
698     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
699    
700       Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
701       checkDudleyError();
702    
    Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );  
    checkFinleyError();  
703    
704     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );     Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
705     checkFinleyError();     checkDudleyError();
706    
707     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );  
708     checkFinleyError();     checkDudleyError();
709  }  }
710    
711  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
# Line 792  void  MeshAdapter::addPDEToLumpedSystem( Line 717  void  MeshAdapter::addPDEToLumpedSystem(
717     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
718     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
719    
720     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
721    
722     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
723     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     checkDudleyError();
724      
725       Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
726       checkDudleyError();
727    
    checkFinleyError();  
728  }  }
729    
730    
731  //  //
732  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
733  //  //
734  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,const  escript::Data& Y, const escript::Data& y) const
735  {  {
736     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
737    
738     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
739     escriptDataC _X=X.getDataC();     escriptDataC _X=X.getDataC();
740     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
741     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
    escriptDataC _y_contact=y_contact.getDataC();  
742    
743     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );     Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
744     checkFinleyError();     checkDudleyError();
745    
746     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );     Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
747     checkFinleyError();     checkDudleyError();
748    
749     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );     checkDudleyError();
    checkFinleyError();  
750  }  }
751  //  //
752  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
# Line 848  void MeshAdapter::addPDEToTransportProbl Line 773  void MeshAdapter::addPDEToTransportProbl
773     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
774     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
775    
776     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
777     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tp.getPaso_TransportProblem();
778    
779     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );     Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
780     checkFinleyError();     checkDudleyError();
781    
782     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );     Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
783     checkFinleyError();     checkDudleyError();
784    
785     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );     Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
786     checkFinleyError();     checkDudleyError();
787    
788     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );     checkDudleyError();
    checkFinleyError();  
789  }  }
790    
791  //  //
# Line 872  void MeshAdapter::interpolateOnDomain(es Line 796  void MeshAdapter::interpolateOnDomain(es
796     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
797     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
798     if (inDomain!=*this)       if (inDomain!=*this)  
799        throw FinleyAdapterException("Error - Illegal domain of interpolant.");        throw DudleyAdapterException("Error - Illegal domain of interpolant.");
800     if (targetDomain!=*this)     if (targetDomain!=*this)
801        throw FinleyAdapterException("Error - Illegal domain of interpolation target.");        throw DudleyAdapterException("Error - Illegal domain of interpolation target.");
802    
803     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
804     escriptDataC _target=target.getDataC();     escriptDataC _target=target.getDataC();
805     escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
806     switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
# Line 886  void MeshAdapter::interpolateOnDomain(es Line 810  void MeshAdapter::interpolateOnDomain(es
810        case(ReducedNodes):        case(ReducedNodes):
811        case(DegreesOfFreedom):        case(DegreesOfFreedom):
812        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
813        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
814        break;        break;
815        case(Elements):        case(Elements):
816        case(ReducedElements):        case(ReducedElements):
817        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
818        break;        break;
819        case(FaceElements):        case(FaceElements):
820        case(ReducedFaceElements):        case(ReducedFaceElements):
821        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
822        break;        break;
823        case(Points):        case(Points):
824        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
       break;  
       case(ContactElementsZero):  
       case(ReducedContactElementsZero):  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
       break;  
       case(ContactElementsOne):  
       case(ReducedContactElementsOne):  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
825        break;        break;
826        default:        default:
827           stringstream temp;           stringstream temp;
828           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
829           throw FinleyAdapterException(temp.str());           throw DudleyAdapterException(temp.str());
830           break;           break;
831        }        }
832        break;        break;
# Line 920  void MeshAdapter::interpolateOnDomain(es Line 836  void MeshAdapter::interpolateOnDomain(es
836        case(ReducedNodes):        case(ReducedNodes):
837        case(DegreesOfFreedom):        case(DegreesOfFreedom):
838        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
839        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
840        break;        break;
841        case(Elements):        case(Elements):
842        case(ReducedElements):        case(ReducedElements):
843        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
844        break;        break;
845        case(FaceElements):        case(FaceElements):
846        case(ReducedFaceElements):        case(ReducedFaceElements):
847        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
848        break;        break;
849        case(Points):        case(Points):
850        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
       break;  
       case(ContactElementsZero):  
       case(ReducedContactElementsZero):  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
       break;  
       case(ContactElementsOne):  
       case(ReducedContactElementsOne):  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
851        break;        break;
852        default:        default:
853           stringstream temp;           stringstream temp;
854           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
855           throw FinleyAdapterException(temp.str());           throw DudleyAdapterException(temp.str());
856           break;           break;
857        }        }
858        break;        break;
859     case(Elements):     case(Elements):
860        if (target.getFunctionSpace().getTypeCode()==Elements) {        if (target.getFunctionSpace().getTypeCode()==Elements) {
861           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Dudley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
862        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
863           Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);           Dudley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
864        } else {        } else {
865           throw FinleyAdapterException("Error - No interpolation with data on elements possible.");           throw DudleyAdapterException("Error - No interpolation with data on elements possible.");
866        }        }
867        break;        break;
868     case(ReducedElements):     case(ReducedElements):
869        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
870           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Dudley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
871        } else {        } else {
872           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");           throw DudleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
873        }        }
874        break;        break;
875     case(FaceElements):     case(FaceElements):
876        if (target.getFunctionSpace().getTypeCode()==FaceElements) {        if (target.getFunctionSpace().getTypeCode()==FaceElements) {
877           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Dudley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
878        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
879           Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);           Dudley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
880        } else {        } else {
881           throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");           throw DudleyAdapterException("Error - No interpolation with data on face elements possible.");
882        }        }
883        break;        break;
884     case(ReducedFaceElements):     case(ReducedFaceElements):
885        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
886           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Dudley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
887        } else {        } else {
888           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");           throw DudleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
889        }        }
890        break;        break;
891     case(Points):     case(Points):
892        if (target.getFunctionSpace().getTypeCode()==Points) {        if (target.getFunctionSpace().getTypeCode()==Points) {
893           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);           Dudley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
       } else {  
          throw FinleyAdapterException("Error - No interpolation with data on points possible.");  
       }  
       break;  
    case(ContactElementsZero):  
    case(ContactElementsOne):  
       if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {  
          Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);  
       } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {  
          Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);  
       } else {  
          throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");  
       }  
       break;  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
       if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {  
          Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);  
894        } else {        } else {
895           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");           throw DudleyAdapterException("Error - No interpolation with data on points possible.");
896        }        }
897        break;        break;
898     case(DegreesOfFreedom):           case(DegreesOfFreedom):      
899        switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
900        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
901        case(DegreesOfFreedom):        case(DegreesOfFreedom):
902        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
903        break;        break;
904        
905        case(Nodes):        case(Nodes):
# Line 1018  void MeshAdapter::interpolateOnDomain(es Line 908  void MeshAdapter::interpolateOnDomain(es
908           escript::Data temp=escript::Data(in);           escript::Data temp=escript::Data(in);
909           temp.expand();           temp.expand();
910           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
911           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);           Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
912        } else {        } else {
913           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);           Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
914        }        }
915        break;        break;
916        case(Elements):        case(Elements):
# Line 1028  void MeshAdapter::interpolateOnDomain(es Line 918  void MeshAdapter::interpolateOnDomain(es
918        if (getMPISize()>1) {        if (getMPISize()>1) {
919           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
920           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
921           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
922        } else {        } else {
923           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
924        }        }
925        break;        break;
926        case(FaceElements):        case(FaceElements):
# Line 1038  void MeshAdapter::interpolateOnDomain(es Line 928  void MeshAdapter::interpolateOnDomain(es
928        if (getMPISize()>1) {        if (getMPISize()>1) {
929           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
930           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
931           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
932        
933        } else {        } else {
934           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
935        }        }
936        break;        break;
937        case(Points):        case(Points):
# Line 1049  void MeshAdapter::interpolateOnDomain(es Line 939  void MeshAdapter::interpolateOnDomain(es
939           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
940           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
941        } else {        } else {
942           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
       }  
       break;  
       case(ContactElementsZero):  
       case(ContactElementsOne):  
       case(ReducedContactElementsZero):  
       case(ReducedContactElementsOne):  
       if (getMPISize()>1) {  
          escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
          escriptDataC _in2 = temp.getDataC();  
          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);  
       } else {  
          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
943        }        }
944        break;        break;
945        default:        default:
946           stringstream temp;           stringstream temp;
947           temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
948           throw FinleyAdapterException(temp.str());           throw DudleyAdapterException(temp.str());
949           break;           break;
950        }        }
951        break;        break;
952     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
953        switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
954        case(Nodes):        case(Nodes):
955        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");        throw DudleyAdapterException("Error - Dudley does not support interpolation from reduced degrees of freedom to mesh nodes.");
956        break;        break;
957        case(ReducedNodes):        case(ReducedNodes):
958        if (getMPISize()>1) {        if (getMPISize()>1) {
959           escript::Data temp=escript::Data(in);           escript::Data temp=escript::Data(in);
960           temp.expand();           temp.expand();
961           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
962           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);           Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
963        } else {        } else {
964           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);           Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
965        }        }
966        break;        break;
967        case(DegreesOfFreedom):        case(DegreesOfFreedom):
968        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");        throw DudleyAdapterException("Error - Dudley does not support interpolation from reduced degrees of freedom to degrees of freedom");
969        break;        break;
970        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
971        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
972        break;        break;
973        case(Elements):        case(Elements):
974        case(ReducedElements):        case(ReducedElements):
975        if (getMPISize()>1) {        if (getMPISize()>1) {
976           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
977           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
978           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
979        } else {        } else {
980           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
981        }        }
982        break;        break;
983        case(FaceElements):        case(FaceElements):
# Line 1107  void MeshAdapter::interpolateOnDomain(es Line 985  void MeshAdapter::interpolateOnDomain(es
985        if (getMPISize()>1) {        if (getMPISize()>1) {
986           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
987           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
988           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
989        } else {        } else {
990           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
991        }        }
992        break;        break;
993        case(Points):        case(Points):
994        if (getMPISize()>1) {        if (getMPISize()>1) {
995           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
996           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
997           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
998        } else {        } else {
999           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
       }  
       break;  
       case(ContactElementsZero):  
       case(ContactElementsOne):  
       case(ReducedContactElementsZero):  
       case(ReducedContactElementsOne):  
       if (getMPISize()>1) {  
          escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );  
          escriptDataC _in2 = temp.getDataC();  
          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);  
       } else {  
          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
1000        }        }
1001        break;        break;
1002        default:        default:
1003           stringstream temp;           stringstream temp;
1004           temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1005           throw FinleyAdapterException(temp.str());           throw DudleyAdapterException(temp.str());
1006           break;           break;
1007        }        }
1008        break;        break;
1009     default:     default:
1010        stringstream temp;        stringstream temp;
1011        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();        temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1012        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1013        break;        break;
1014     }     }
1015     checkFinleyError();     checkDudleyError();
1016  }  }
1017    
1018  //  //
# Line 1156  void MeshAdapter::setToX(escript::Data& Line 1022  void MeshAdapter::setToX(escript::Data&
1022  {  {
1023     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1024     if (argDomain!=*this)     if (argDomain!=*this)
1025        throw FinleyAdapterException("Error - Illegal domain of data point locations");        throw DudleyAdapterException("Error - Illegal domain of data point locations");
1026     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1027     // in case of values node coordinates we can do the job directly:     // in case of values node coordinates we can do the job directly:
1028     if (arg.getFunctionSpace().getTypeCode()==Nodes) {     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1029        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1030        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Dudley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1031     } else {     } else {
1032        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
1033        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1034        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Dudley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1035        // this is then interpolated onto arg:        // this is then interpolated onto arg:
1036        interpolateOnDomain(arg,tmp_data);        interpolateOnDomain(arg,tmp_data);
1037     }     }
1038     checkFinleyError();     checkDudleyError();
1039  }  }
1040    
1041  //  //
# Line 1180  void MeshAdapter::setToNormal(escript::D Line 1046  void MeshAdapter::setToNormal(escript::D
1046  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1047     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1048     if (normalDomain!=*this)     if (normalDomain!=*this)
1049        throw FinleyAdapterException("Error - Illegal domain of normal locations");        throw DudleyAdapterException("Error - Illegal domain of normal locations");
1050     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1051     escriptDataC _normal=normal.getDataC();     escriptDataC _normal=normal.getDataC();
1052     switch(normal.getFunctionSpace().getTypeCode()) {     switch(normal.getFunctionSpace().getTypeCode()) {
1053     case(Nodes):     case(Nodes):
1054     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for nodes");
1055     break;     break;
1056     case(ReducedNodes):     case(ReducedNodes):
1057     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for reduced nodes");
1058     break;     break;
1059     case(Elements):     case(Elements):
1060     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for elements");
1061     break;     break;
1062     case(ReducedElements):     case(ReducedElements):
1063     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for elements with reduced integration order");
1064     break;     break;
1065     case (FaceElements):     case (FaceElements):
1066     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     Dudley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1067     break;     break;
1068     case (ReducedFaceElements):     case (ReducedFaceElements):
1069     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     Dudley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1070     break;     break;
1071     case(Points):     case(Points):
1072     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for point elements");
    break;  
    case (ContactElementsOne):  
    case (ContactElementsZero):  
    Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);  
    break;  
    case (ReducedContactElementsOne):  
    case (ReducedContactElementsZero):  
    Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);  
1073     break;     break;
1074     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1075     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for degrees of freedom.");
1076     break;     break;
1077     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1078     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for reduced degrees of freedom.");
1079     break;     break;
1080     default:     default:
1081        stringstream temp;        stringstream temp;
1082        temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();        temp << "Error - Normal Vectors: Dudley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1083        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1084        break;        break;
1085     }     }
1086     checkFinleyError();     checkDudleyError();
1087  }  }
1088    
1089  //  //
# Line 1236  void MeshAdapter::interpolateACross(escr Line 1094  void MeshAdapter::interpolateACross(escr
1094     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1095     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1096     if (targetDomain!=this)     if (targetDomain!=this)
1097        throw FinleyAdapterException("Error - Illegal domain of interpolation target");        throw DudleyAdapterException("Error - Illegal domain of interpolation target");
1098    
1099     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");     throw DudleyAdapterException("Error - Dudley does not allow interpolation across domains yet.");
1100  }  }
1101    
1102  //  //
1103  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
1104  //  //
1105  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1106  {  {
1107     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1108     if (argDomain!=*this)     if (argDomain!=*this)
1109        throw FinleyAdapterException("Error - Illegal domain of integration kernel");        throw DudleyAdapterException("Error - Illegal domain of integration kernel");
1110    
1111     double blocktimer_start = blocktimer_time();     double blocktimer_start = blocktimer_time();
1112     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1113     escriptDataC _temp;     escriptDataC _temp;
1114     escript::Data temp;     escript::Data temp;
1115     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
# Line 1259  void MeshAdapter::setToIntegrals(std::ve Line 1117  void MeshAdapter::setToIntegrals(std::ve
1117     case(Nodes):     case(Nodes):
1118     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1119     _temp=temp.getDataC();     _temp=temp.getDataC();
1120     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1121     break;     break;
1122     case(ReducedNodes):     case(ReducedNodes):
1123     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1124     _temp=temp.getDataC();     _temp=temp.getDataC();
1125     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1126     break;     break;
1127     case(Elements):     case(Elements):
1128     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);     Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1129     break;     break;
1130     case(ReducedElements):     case(ReducedElements):
1131     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);     Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1132     break;     break;
1133     case(FaceElements):     case(FaceElements):
1134     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);     Dudley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1135     break;     break;
1136     case(ReducedFaceElements):     case(ReducedFaceElements):
1137     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);     Dudley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1138     break;     break;
1139     case(Points):     case(Points):
1140     throw FinleyAdapterException("Error - Integral of data on points is not supported.");     throw DudleyAdapterException("Error - Integral of data on points is not supported.");
    break;  
    case(ContactElementsZero):  
    Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);  
    break;  
    case(ReducedContactElementsZero):  
    Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);  
    break;  
    case(ContactElementsOne):  
    Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);  
    break;  
    case(ReducedContactElementsOne):  
    Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);  
1141     break;     break;
1142     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1143     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1144     _temp=temp.getDataC();     _temp=temp.getDataC();
1145     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1146     break;     break;
1147     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1148     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1149     _temp=temp.getDataC();     _temp=temp.getDataC();
1150     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1151     break;     break;
1152     default:     default:
1153        stringstream temp;        stringstream temp;
1154        temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();        temp << "Error - Integrals: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1155        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1156        break;        break;
1157     }     }
1158     checkFinleyError();     checkDudleyError();
1159     blocktimer_increment("integrate()", blocktimer_start);     blocktimer_increment("integrate()", blocktimer_start);
1160  }  }
1161    
# Line 1320  void MeshAdapter::setToGradient(escript: Line 1166  void MeshAdapter::setToGradient(escript:
1166  {  {
1167     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1168     if (argDomain!=*this)     if (argDomain!=*this)
1169        throw FinleyAdapterException("Error - Illegal domain of gradient argument");        throw DudleyAdapterException("Error - Illegal domain of gradient argument");
1170     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1171     if (gradDomain!=*this)     if (gradDomain!=*this)
1172        throw FinleyAdapterException("Error - Illegal domain of gradient");        throw DudleyAdapterException("Error - Illegal domain of gradient");
1173    
1174     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1175     escriptDataC _grad=grad.getDataC();     escriptDataC _grad=grad.getDataC();
1176     escriptDataC nodeDataC;     escriptDataC nodeDataC;
1177     escript::Data temp;     escript::Data temp;
# Line 1344  void MeshAdapter::setToGradient(escript: Line 1190  void MeshAdapter::setToGradient(escript:
1190     }     }
1191     switch(grad.getFunctionSpace().getTypeCode()) {     switch(grad.getFunctionSpace().getTypeCode()) {
1192     case(Nodes):     case(Nodes):
1193     throw FinleyAdapterException("Error - Gradient at nodes is not supported.");     throw DudleyAdapterException("Error - Gradient at nodes is not supported.");
1194     break;     break;
1195     case(ReducedNodes):     case(ReducedNodes):
1196     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");     throw DudleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1197     break;     break;
1198     case(Elements):     case(Elements):
1199     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);     Dudley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1200     break;     break;
1201     case(ReducedElements):     case(ReducedElements):
1202     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);     Dudley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1203     break;     break;
1204     case(FaceElements):     case(FaceElements):
1205     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);     Dudley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1206     break;     break;
1207     case(ReducedFaceElements):     case(ReducedFaceElements):
1208     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);     Dudley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1209     break;     break;
1210     case(Points):     case(Points):
1211     throw FinleyAdapterException("Error - Gradient at points is not supported.");     throw DudleyAdapterException("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);  
1212     break;     break;
1213     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1214     throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");     throw DudleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1215     break;     break;
1216     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1217     throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");     throw DudleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1218     break;     break;
1219     default:     default:
1220        stringstream temp;        stringstream temp;
1221        temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();        temp << "Error - Gradient: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1222        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1223        break;        break;
1224     }     }
1225     checkFinleyError();     checkDudleyError();
1226  }  }
1227    
1228  //  //
# Line 1396  void MeshAdapter::setToGradient(escript: Line 1230  void MeshAdapter::setToGradient(escript:
1230  //  //
1231  void MeshAdapter::setToSize(escript::Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
1232  {  {
1233     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1234     escriptDataC tmp=size.getDataC();     escriptDataC tmp=size.getDataC();
1235     switch(size.getFunctionSpace().getTypeCode()) {     switch(size.getFunctionSpace().getTypeCode()) {
1236     case(Nodes):     case(Nodes):
1237     throw FinleyAdapterException("Error - Size of nodes is not supported.");     throw DudleyAdapterException("Error - Size of nodes is not supported.");
1238     break;     break;
1239     case(ReducedNodes):     case(ReducedNodes):
1240     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");     throw DudleyAdapterException("Error - Size of reduced nodes is not supported.");
1241     break;     break;
1242     case(Elements):     case(Elements):
1243     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     Dudley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1244     break;     break;
1245     case(ReducedElements):     case(ReducedElements):
1246     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     Dudley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1247     break;     break;
1248     case(FaceElements):     case(FaceElements):
1249     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Dudley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1250     break;     break;
1251     case(ReducedFaceElements):     case(ReducedFaceElements):
1252     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Dudley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1253     break;     break;
1254     case(Points):     case(Points):
1255     throw FinleyAdapterException("Error - Size of point elements is not supported.");     throw DudleyAdapterException("Error - Size of point elements is not supported.");
    break;  
    case(ContactElementsZero):  
    case(ContactElementsOne):  
    Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);  
    break;  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);  
1256     break;     break;
1257     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1258     throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");     throw DudleyAdapterException("Error - Size of degrees of freedom is not supported.");
1259     break;     break;
1260     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1261     throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");     throw DudleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1262     break;     break;
1263     default:     default:
1264        stringstream temp;        stringstream temp;
1265        temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();        temp << "Error - Element size: Dudley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1266        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1267        break;        break;
1268     }     }
1269     checkFinleyError();     checkDudleyError();
1270  }  }
1271    
1272  //  //
# Line 1448  void MeshAdapter::setToSize(escript::Dat Line 1274  void MeshAdapter::setToSize(escript::Dat
1274  //  //
1275  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1276  {  {
1277     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1278     escriptDataC tmp;     escriptDataC tmp;
1279     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1280     if (newDomain!=*this)     if (newDomain!=*this)
1281        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw DudleyAdapterException("Error - Illegal domain of new point locations");
1282     tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {
1283     Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1284     checkFinleyError();         Dudley_Mesh_setCoordinates(mesh,&tmp);
1285       } else {
1286           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );
1287           tmp = new_x_inter.getDataC();
1288           Dudley_Mesh_setCoordinates(mesh,&tmp);
1289       }
1290       checkDudleyError();
1291  }  }
1292    
1293  //  //
# Line 1472  void MeshAdapter::extractArgsFromDict(co Line 1304  void MeshAdapter::extractArgsFromDict(co
1304    
1305     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1306     for (int i=0; i<numData; ++i) {     for (int i=0; i<numData; ++i) {
1307        std::string n=boost::python::extract<std::string>(keys[i]);        string n=boost::python::extract<string>(keys[i]);
1308        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1309        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1310           throw FinleyAdapterException("Error: Data must be defined on same Domain");           throw DudleyAdapterException("Error: Data must be defined on same Domain");
1311        data[i] = d.getDataC();        data[i] = d.getDataC();
1312        dataPtr[i] = &(data[i]);        dataPtr[i] = &(data[i]);
1313        names[i] = TMPMEMALLOC(n.length()+1, char);        names[i] = TMPMEMALLOC(n.length()+1, char);
# Line 1486  void MeshAdapter::extractArgsFromDict(co Line 1318  void MeshAdapter::extractArgsFromDict(co
1318  //  //
1319  // saves mesh and optionally data arrays in openDX format  // saves mesh and optionally data arrays in openDX format
1320  //  //
1321  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1322  {  {
1323     int num_data;     int num_data;
1324     char **names;     char **names;
# Line 1494  void MeshAdapter::saveDX(const std::stri Line 1326  void MeshAdapter::saveDX(const std::stri
1326     escriptDataC **ptr_data;     escriptDataC **ptr_data;
1327    
1328     extractArgsFromDict(arg, num_data, names, data, ptr_data);     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1329     Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);     Dudley_Mesh_saveDX(filename.c_str(), m_dudleyMesh.get(), num_data, names, ptr_data);
1330     checkFinleyError();     checkDudleyError();
1331    
1332     /* win32 refactor */     /* win32 refactor */
1333     TMPMEMFREE(data);     TMPMEMFREE(data);
# Line 1512  void MeshAdapter::saveDX(const std::stri Line 1344  void MeshAdapter::saveDX(const std::stri
1344  //  //
1345  // saves mesh and optionally data arrays in VTK format  // saves mesh and optionally data arrays in VTK format
1346  //  //
1347  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1348  {  {
1349     int num_data;     int num_data;
1350     char **names;     char **names;
# Line 1520  void MeshAdapter::saveVTK(const std::str Line 1352  void MeshAdapter::saveVTK(const std::str
1352     escriptDataC **ptr_data;     escriptDataC **ptr_data;
1353    
1354     extractArgsFromDict(arg, num_data, names, data, ptr_data);     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1355     Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);     Dudley_Mesh_saveVTK(filename.c_str(), m_dudleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());
1356     checkFinleyError();     checkDudleyError();
1357    
1358     /* win32 refactor */     /* win32 refactor */
1359     TMPMEMFREE(data);     TMPMEMFREE(data);
# Line 1533  void MeshAdapter::saveVTK(const std::str Line 1365  void MeshAdapter::saveVTK(const std::str
1365     TMPMEMFREE(names);     TMPMEMFREE(names);
1366  }  }
1367    
1368    bool MeshAdapter::ownSample(int fs_code, index_t id) const
1369    {
1370    #ifdef ESYS_MPI
1371        index_t myFirstNode=0, myLastNode=0, k=0;
1372        index_t* globalNodeIndex=0;
1373        Dudley_Mesh* mesh_p=m_dudleyMesh.get();
1374        if (fs_code == DUDLEY_REDUCED_NODES)
1375        {
1376        myFirstNode = Dudley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1377        myLastNode = Dudley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1378        globalNodeIndex = Dudley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1379        }
1380        else
1381        {
1382        myFirstNode = Dudley_NodeFile_getFirstNode(mesh_p->Nodes);
1383        myLastNode = Dudley_NodeFile_getLastNode(mesh_p->Nodes);
1384        globalNodeIndex = Dudley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1385        }
1386        k=globalNodeIndex[id];
1387        return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1388    #endif
1389        return true;
1390    }
1391    
1392    
1393    
1394  //  //
1395  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1396  //  //
1397  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  ASM_ptr MeshAdapter::newSystemMatrix(
1398                                                   const int row_blocksize,                                                   const int row_blocksize,
1399                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1400                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1548  SystemMatrixAdapter MeshAdapter::newSyst Line 1406  SystemMatrixAdapter MeshAdapter::newSyst
1406     // is the domain right?     // is the domain right?
1407     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1408     if (row_domain!=*this)     if (row_domain!=*this)
1409        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");        throw DudleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1410     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1411     if (col_domain!=*this)     if (col_domain!=*this)
1412        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        throw DudleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1413     // is the function space type right     // is the function space type right
1414     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1415        reduceRowOrder=0;        reduceRowOrder=0;
1416     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1417        reduceRowOrder=1;        reduceRowOrder=1;
1418     } else {     } else {
1419        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");        throw DudleyAdapterException("Error - illegal function space type for system matrix rows.");
1420     }     }
1421     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1422        reduceColOrder=0;        reduceColOrder=0;
1423     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1424        reduceColOrder=1;        reduceColOrder=1;
1425     } else {     } else {
1426        throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");        throw DudleyAdapterException("Error - illegal function space type for system matrix columns.");
1427     }     }
1428     // generate matrix:     // generate matrix:
1429    
1430     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Dudley_getPattern(getDudley_Mesh(),reduceRowOrder,reduceColOrder);
1431     checkFinleyError();     checkDudleyError();
1432     Paso_SystemMatrix* fsystemMatrix;     Paso_SystemMatrix* fsystemMatrix;
1433     int trilinos = 0;     int trilinos = 0;
1434     if (trilinos) {     if (trilinos) {
# Line 1579  SystemMatrixAdapter MeshAdapter::newSyst Line 1437  SystemMatrixAdapter MeshAdapter::newSyst
1437  #endif  #endif
1438     }     }
1439     else {     else {
1440        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1441     }     }
1442     checkPasoError();     checkPasoError();
1443     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1444     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace, column_blocksize,column_functionspace);
1445       return ASM_ptr(sma);
1446  }  }
1447    
1448  //  //
1449  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1450  //  //
1451  TransportProblemAdapter MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
1452                                                           const double theta,                                                           const bool useBackwardEuler,
1453                                                           const int blocksize,                                                           const int blocksize,
1454                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1455                                                           const int type) const                                                           const int type) const
# Line 1599  TransportProblemAdapter MeshAdapter::new Line 1458  TransportProblemAdapter MeshAdapter::new
1458     // is the domain right?     // is the domain right?
1459     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1460     if (domain!=*this)     if (domain!=*this)
1461        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");        throw DudleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1462     // is the function space type right     // is the function space type right
1463     if (functionspace.getTypeCode()==DegreesOfFreedom) {     if (functionspace.getTypeCode()==DegreesOfFreedom) {
1464        reduceOrder=0;        reduceOrder=0;
1465     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1466        reduceOrder=1;        reduceOrder=1;
1467     } else {     } else {
1468        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");        throw DudleyAdapterException("Error - illegal function space type for system matrix rows.");
1469     }     }
1470     // generate matrix:     // generate matrix:
1471    
1472     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Dudley_getPattern(getDudley_Mesh(),reduceOrder,reduceOrder);
1473     checkFinleyError();     checkDudleyError();
1474     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1475     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1476     checkPasoError();     checkPasoError();
1477     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1478     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     AbstractTransportProblem* atp=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1479       return ATP_ptr(atp);
1480  }  }
1481    
1482  //  //
# Line 1637  bool MeshAdapter::isCellOriented(int fun Line 1497  bool MeshAdapter::isCellOriented(int fun
1497     case(Elements):     case(Elements):
1498     case(FaceElements):     case(FaceElements):
1499     case(Points):     case(Points):
    case(ContactElementsZero):  
    case(ContactElementsOne):  
1500     case(ReducedElements):     case(ReducedElements):
1501     case(ReducedFaceElements):     case(ReducedFaceElements):
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
1502     return true;     return true;
1503     break;     break;
1504     default:     default:
1505        stringstream temp;        stringstream temp;
1506        temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;        temp << "Error - Cell: Dudley does not know anything about function space type " << functionSpaceCode;
1507        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1508        break;        break;
1509     }     }
1510     checkFinleyError();     checkDudleyError();
1511     return false;     return false;
1512  }  }
1513    
1514    bool
1515    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1516    {
1517       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1518        class 1: DOF <-> Nodes
1519        class 2: ReducedDOF <-> ReducedNodes
1520        class 3: Points
1521        class 4: Elements
1522        class 5: ReducedElements
1523        class 6: FaceElements
1524        class 7: ReducedFaceElements
1525        class 8: ContactElementZero <-> ContactElementOne
1526        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1527    
1528       There is also a set of lines. Interpolation is possible down a line but not between lines.
1529       class 1 and 2 belong to all lines so aren't considered.
1530        line 0: class 3
1531        line 1: class 4,5
1532        line 2: class 6,7
1533        line 3: class 8,9
1534    
1535       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1536       eg hasnodes is true if we have at least one instance of Nodes.
1537       */
1538        if (fs.empty())
1539        {
1540            return false;
1541        }
1542        vector<int> hasclass(10);
1543        vector<int> hasline(4);
1544        bool hasnodes=false;
1545        bool hasrednodes=false;
1546        for (int i=0;i<fs.size();++i)
1547        {
1548        switch(fs[i])
1549        {
1550        case(Nodes):   hasnodes=true;   // no break is deliberate
1551        case(DegreesOfFreedom):
1552            hasclass[1]=1;
1553            break;
1554        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1555        case(ReducedDegreesOfFreedom):
1556            hasclass[2]=1;
1557            break;
1558        case(Points):
1559            hasline[0]=1;
1560            hasclass[3]=1;
1561            break;
1562        case(Elements):
1563            hasclass[4]=1;
1564            hasline[1]=1;
1565            break;
1566        case(ReducedElements):
1567            hasclass[5]=1;
1568            hasline[1]=1;
1569            break;
1570        case(FaceElements):
1571            hasclass[6]=1;
1572            hasline[2]=1;
1573            break;
1574        case(ReducedFaceElements):
1575            hasclass[7]=1;
1576            hasline[2]=1;
1577            break;
1578        default:
1579            return false;
1580        }
1581        }
1582        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1583        // fail if we have more than one leaf group
1584    
1585        if (totlines>1)
1586        {
1587        return false;   // there are at least two branches we can't interpolate between
1588        }
1589        else if (totlines==1)
1590        {
1591        if (hasline[0]==1)      // we have points
1592        {
1593            resultcode=Points;
1594        }
1595        else if (hasline[1]==1)
1596        {
1597            if (hasclass[5]==1)
1598            {
1599            resultcode=ReducedElements;
1600            }
1601            else
1602            {
1603            resultcode=Elements;
1604            }
1605        }
1606        else if (hasline[2]==1)
1607        {
1608            if (hasclass[7]==1)
1609            {
1610            resultcode=ReducedFaceElements;
1611            }
1612            else
1613            {
1614            resultcode=FaceElements;
1615            }
1616        }
1617        else    // so we must be in line3
1618        {
1619    
1620            throw DudleyAdapterException("Programmer Error - choosing between contact elements - we should never get here.");
1621    
1622        }
1623        }
1624        else    // totlines==0
1625        {
1626        if (hasclass[2]==1)
1627        {
1628            // something from class 2
1629            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1630        }
1631        else
1632        {   // something from class 1
1633            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1634        }
1635        }
1636        return true;
1637    }
1638    
1639  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1640  {  {
1641     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1642     case(Nodes):     case(Nodes):
1643     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1644     case(Nodes):      case(Nodes):
1645     case(ReducedNodes):      case(ReducedNodes):
1646     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1647     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1648     case(Elements):      case(Elements):
1649     case(ReducedElements):      case(ReducedElements):
1650     case(FaceElements):      case(FaceElements):
1651     case(ReducedFaceElements):      case(ReducedFaceElements):
1652     case(Points):      case(Points):
1653     case(ContactElementsZero):      return true;
1654     case(ReducedContactElementsZero):      default:
1655     case(ContactElementsOne):            stringstream temp;
1656     case(ReducedContactElementsOne):            temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1657     return true;            throw DudleyAdapterException(temp.str());
    default:  
       stringstream temp;  
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;  
       throw FinleyAdapterException(temp.str());  
1658     }     }
1659     break;     break;
1660     case(ReducedNodes):     case(ReducedNodes):
1661     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1662     case(ReducedNodes):      case(ReducedNodes):
1663     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1664     case(Elements):      case(Elements):
1665     case(ReducedElements):      case(ReducedElements):
1666     case(FaceElements):      case(FaceElements):
1667     case(ReducedFaceElements):      case(ReducedFaceElements):
1668     case(Points):      case(Points):
1669     case(ContactElementsZero):      return true;
1670     case(ReducedContactElementsZero):      case(Nodes):
1671     case(ContactElementsOne):      case(DegreesOfFreedom):
1672     case(ReducedContactElementsOne):      return false;
1673     return true;      default:
1674     case(Nodes):          stringstream temp;
1675     case(DegreesOfFreedom):          temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1676     return false;          throw DudleyAdapterException(temp.str());
    default:  
       stringstream temp;  
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;  
       throw FinleyAdapterException(temp.str());  
1677     }     }
1678     break;     break;
1679     case(Elements):     case(Elements):
1680     if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1681        return true;        return true;
1682     } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1683        return true;        return true;
1684     } else {          } else {
1685        return false;            return false;
1686     }          }
1687     case(ReducedElements):     case(ReducedElements):
1688     if (functionSpaceType_target==ReducedElements) {      if (functionSpaceType_target==ReducedElements) {
1689        return true;        return true;
1690     } else {      } else {
1691        return false;            return false;
1692     }      }
1693     case(FaceElements):     case(FaceElements):
1694     if (functionSpaceType_target==FaceElements) {      if (functionSpaceType_target==FaceElements) {
1695        return true;              return true;
1696     } else if (functionSpaceType_target==ReducedFaceElements) {      } else if (functionSpaceType_target==ReducedFaceElements) {
1697        return true;              return true;
1698     } else {      } else {
1699        return false;              return false;
1700     }      }
1701     case(ReducedFaceElements):     case(ReducedFaceElements):
1702     if (functionSpaceType_target==ReducedFaceElements) {      if (functionSpaceType_target==ReducedFaceElements) {
1703        return true;              return true;
1704     } else {      } else {
1705        return false;          return false;
1706     }      }
1707     case(Points):     case(Points):
1708     if (functionSpaceType_target==Points) {      if (functionSpaceType_target==Points) {
1709        return true;              return true;
1710     } else {      } else {
1711        return false;              return false;
1712     }      }
    case(ContactElementsZero):  
    case(ContactElementsOne):  
    if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {  
       return true;  
    } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {  
       return true;  
    } else {  
       return false;  
    }  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {  
       return true;  
    } else {  
       return false;  
    }  
1713     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1714     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1715     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1716     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1717     case(Nodes):      case(Nodes):
1718     case(ReducedNodes):      case(ReducedNodes):
1719     case(Elements):      case(Elements):
1720     case(ReducedElements):      case(ReducedElements):
1721     case(Points):      case(Points):
1722     case(FaceElements):      case(FaceElements):
1723     case(ReducedFaceElements):      case(ReducedFaceElements):
1724     case(ContactElementsZero):      return true;
1725     case(ReducedContactElementsZero):      default:
1726     case(ContactElementsOne):          stringstream temp;
1727     case(ReducedContactElementsOne):          temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1728     return true;          throw DudleyAdapterException(temp.str());
1729     default:      }
1730        stringstream temp;      break;
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;  
       throw FinleyAdapterException(temp.str());  
    }  
    break;  
1731     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1732     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
1733     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1734     case(ReducedNodes):      case(ReducedNodes):
1735     case(Elements):      case(Elements):
1736     case(ReducedElements):      case(ReducedElements):
1737     case(FaceElements):      case(FaceElements):
1738     case(ReducedFaceElements):      case(ReducedFaceElements):
1739     case(Points):      case(Points):
1740     case(ContactElementsZero):      return true;
1741     case(ReducedContactElementsZero):      case(Nodes):
1742     case(ContactElementsOne):      case(DegreesOfFreedom):
1743     case(ReducedContactElementsOne):      return false;
1744     return true;      default:
1745     case(Nodes):          stringstream temp;
1746     case(DegreesOfFreedom):          temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1747     return false;          throw DudleyAdapterException(temp.str());
1748     default:      }
1749        stringstream temp;      break;
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;  
       throw FinleyAdapterException(temp.str());  
    }  
    break;  
1750     default:     default:
1751        stringstream temp;        stringstream temp;
1752        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;        temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_source;
1753        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1754        break;        break;
1755     }     }
1756     checkFinleyError();     checkDudleyError();
1757     return false;     return false;
1758  }  }
1759    
# Line 1817  bool MeshAdapter::operator==(const Abstr Line 1766  bool MeshAdapter::operator==(const Abstr
1766  {  {
1767     const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);     const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1768     if (temp!=0) {     if (temp!=0) {
1769        return (m_finleyMesh==temp->m_finleyMesh);        return (m_dudleyMesh==temp->m_dudleyMesh);
1770     } else {     } else {
1771        return false;        return false;
1772     }     }
# Line 1830  bool MeshAdapter::operator!=(const Abstr Line 1779  bool MeshAdapter::operator!=(const Abstr
1779    
1780  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
1781  {  {
1782     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1783     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1784     checkPasoError();     checkPasoError();
1785     return out;     return out;
1786  }  }
1787  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
1788  {  {
1789     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1790     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1791     checkPasoError();     checkPasoError();
1792     return out;     return out;
1793  }  }
# Line 1858  escript::Data MeshAdapter::getSize() con Line 1807  escript::Data MeshAdapter::getSize() con
1807     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(asAbstractContinuousDomain()).getSize();
1808  }  }
1809    
1810  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1811  {  {
1812     int *out = NULL;     int *out = NULL;
1813     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1814     switch (functionSpaceType) {     switch (functionSpaceType) {
1815     case(Nodes):     case(Nodes):
1816     out=mesh->Nodes->Id;     out=mesh->Nodes->Id;
# Line 1884  int* MeshAdapter::borrowSampleReferenceI Line 1833  int* MeshAdapter::borrowSampleReferenceI
1833     case(Points):     case(Points):
1834     out=mesh->Points->Id;     out=mesh->Points->Id;
1835     break;     break;
    case(ContactElementsZero):  
    out=mesh->ContactElements->Id;  
    break;  
    case(ReducedContactElementsZero):  
    out=mesh->ContactElements->Id;  
    break;  
    case(ContactElementsOne):  
    out=mesh->ContactElements->Id;  
    break;  
    case(ReducedContactElementsOne):  
    out=mesh->ContactElements->Id;  
    break;  
1836     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1837     out=mesh->Nodes->degreesOfFreedomId;     out=mesh->Nodes->degreesOfFreedomId;
1838     break;     break;
# Line 1905  int* MeshAdapter::borrowSampleReferenceI Line 1842  int* MeshAdapter::borrowSampleReferenceI
1842     default:     default:
1843        stringstream temp;        stringstream temp;
1844        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1845        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1846        break;        break;
1847     }     }
1848     return out;     return out;
# Line 1913  int* MeshAdapter::borrowSampleReferenceI Line 1850  int* MeshAdapter::borrowSampleReferenceI
1850  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1851  {  {
1852     int out=0;     int out=0;
1853     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1854     switch (functionSpaceType) {     switch (functionSpaceType) {
1855     case(Nodes):     case(Nodes):
1856     out=mesh->Nodes->Tag[sampleNo];     out=mesh->Nodes->Tag[sampleNo];
1857     break;     break;
1858     case(ReducedNodes):     case(ReducedNodes):
1859     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");     throw DudleyAdapterException(" Error - ReducedNodes does not support tags.");
1860     break;     break;
1861     case(Elements):     case(Elements):
1862     out=mesh->Elements->Tag[sampleNo];     out=mesh->Elements->Tag[sampleNo];
# Line 1936  int MeshAdapter::getTagFromSampleNo(int Line 1873  int MeshAdapter::getTagFromSampleNo(int
1873     case(Points):     case(Points):
1874     out=mesh->Points->Tag[sampleNo];     out=mesh->Points->Tag[sampleNo];
1875     break;     break;
    case(ContactElementsZero):  
    out=mesh->ContactElements->Tag[sampleNo];  
    break;  
    case(ReducedContactElementsZero):  
    out=mesh->ContactElements->Tag[sampleNo];  
    break;  
    case(ContactElementsOne):  
    out=mesh->ContactElements->Tag[sampleNo];  
    break;  
    case(ReducedContactElementsOne):  
    out=mesh->ContactElements->Tag[sampleNo];  
    break;  
1876     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1877     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");     throw DudleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1878     break;     break;
1879     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1880     throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");     throw DudleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1881     break;     break;
1882     default:     default:
1883        stringstream temp;        stringstream temp;
1884        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1885        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1886        break;        break;
1887     }     }
1888     return out;     return out;
# Line 1966  int MeshAdapter::getTagFromSampleNo(int Line 1891  int MeshAdapter::getTagFromSampleNo(int
1891    
1892  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
1893  {  {
1894     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1895     escriptDataC tmp=mask.getDataC();     escriptDataC tmp=mask.getDataC();
1896     switch(functionSpaceType) {     switch(functionSpaceType) {
1897     case(Nodes):     case(Nodes):
1898     Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);     Dudley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1899     break;     break;
1900     case(ReducedNodes):     case(ReducedNodes):
1901     throw FinleyAdapterException("Error - ReducedNodes does not support tags");     throw DudleyAdapterException("Error - ReducedNodes does not support tags");
1902     break;     break;
1903     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1904     throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");     throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
1905     break;     break;
1906     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1907     throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");     throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1908     break;     break;
1909     case(Elements):     case(Elements):
1910     Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     Dudley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1911     break;     break;
1912     case(ReducedElements):     case(ReducedElements):
1913     Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     Dudley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1914     break;     break;
1915     case(FaceElements):     case(FaceElements):
1916     Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);     Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1917     break;     break;
1918     case(ReducedFaceElements):     case(ReducedFaceElements):
1919     Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);     Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1920     break;     break;
1921     case(Points):     case(Points):
1922     Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);     Dudley_ElementFile_setTags(mesh->Points,newTag,&tmp);
    break;  
    case(ContactElementsZero):  
    Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);  
    break;  
    case(ReducedContactElementsZero):  
    Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);  
    break;  
    case(ContactElementsOne):  
    Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);  
    break;  
    case(ReducedContactElementsOne):  
    Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);  
1923     break;     break;
1924     default:     default:
1925        stringstream temp;        stringstream temp;
1926        temp << "Error - Finley does not know anything about function space type " << functionSpaceType;        temp << "Error - Dudley does not know anything about function space type " << functionSpaceType;
1927        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1928     }     }
1929     checkFinleyError();     checkDudleyError();
1930     return;     return;
1931  }  }
1932    
1933  void MeshAdapter::setTagMap(const std::string& name,  int tag)  void MeshAdapter::setTagMap(const string& name,  int tag)
1934  {  {
1935     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1936     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Dudley_Mesh_addTagMap(mesh, name.c_str(),tag);
1937     checkPasoError();     checkPasoError();
1938     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1939  }  }
1940    
1941  int MeshAdapter::getTag(const std::string& name) const  int MeshAdapter::getTag(const string& name) const
1942  {  {
1943     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1944     int tag=0;     int tag=0;
1945     tag=Finley_Mesh_getTag(mesh, name.c_str());     tag=Dudley_Mesh_getTag(mesh, name.c_str());
1946     checkPasoError();     checkPasoError();
1947     // throwStandardException("MeshAdapter::getTag is not implemented.");     // throwStandardException("MeshAdapter::getTag is not implemented.");
1948     return tag;     return tag;
1949  }  }
1950    
1951  bool MeshAdapter::isValidTagName(const std::string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
1952  {  {
1953     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1954     return Finley_Mesh_isValidTagName(mesh,name.c_str());     return Dudley_Mesh_isValidTagName(mesh,name.c_str());
1955  }  }
1956    
1957  std::string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
1958  {  {
1959     stringstream temp;     stringstream temp;
1960     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1961     Finley_TagMap* tag_map=mesh->TagMap;     Dudley_TagMap* tag_map=mesh->TagMap;
1962     while (tag_map) {     while (tag_map) {
1963        temp << tag_map->name;        temp << tag_map->name;
1964        tag_map=tag_map->next;        tag_map=tag_map->next;
# Line 2056  std::string MeshAdapter::showTagNames() Line 1969  std::string MeshAdapter::showTagNames()
1969    
1970  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
1971  {  {
1972    Finley_Mesh* mesh=m_finleyMesh.get();    Dudley_Mesh* mesh=m_dudleyMesh.get();
1973    dim_t numTags=0;    dim_t numTags=0;
1974    switch(functionSpaceCode) {    switch(functionSpaceCode) {
1975     case(Nodes):     case(Nodes):
1976            numTags=mesh->Nodes->numTagsInUse;            numTags=mesh->Nodes->numTagsInUse;
1977            break;            break;
1978     case(ReducedNodes):     case(ReducedNodes):
1979            throw FinleyAdapterException("Error - ReducedNodes does not support tags");            throw DudleyAdapterException("Error - ReducedNodes does not support tags");
1980            break;            break;
1981     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1982            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");            throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
1983            break;            break;
1984     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1985            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");            throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1986            break;            break;
1987     case(Elements):     case(Elements):
1988     case(ReducedElements):     case(ReducedElements):
# Line 2082  int MeshAdapter::getNumberOfTagsInUse(in Line 1995  int MeshAdapter::getNumberOfTagsInUse(in
1995     case(Points):     case(Points):
1996            numTags=mesh->Points->numTagsInUse;            numTags=mesh->Points->numTagsInUse;
1997            break;            break;
    case(ContactElementsZero):  
    case(ReducedContactElementsZero):  
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
           numTags=mesh->ContactElements->numTagsInUse;  
           break;  
1998     default:     default:
1999        stringstream temp;        stringstream temp;
2000        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;        temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2001        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
2002    }    }
2003    return numTags;    return numTags;
2004  }  }
2005  int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  
2006    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2007  {  {
2008    Finley_Mesh* mesh=m_finleyMesh.get();    Dudley_Mesh* mesh=m_dudleyMesh.get();
2009    index_t* tags=NULL;    index_t* tags=NULL;
2010    switch(functionSpaceCode) {    switch(functionSpaceCode) {
2011     case(Nodes):     case(Nodes):
2012            tags=mesh->Nodes->tagsInUse;            tags=mesh->Nodes->tagsInUse;
2013            break;            break;
2014     case(ReducedNodes):     case(ReducedNodes):
2015            throw FinleyAdapterException("Error - ReducedNodes does not support tags");            throw DudleyAdapterException("Error - ReducedNodes does not support tags");
2016            break;            break;
2017     case(DegreesOfFreedom):     case(DegreesOfFreedom):
2018            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");            throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
2019            break;            break;
2020     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2021            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");            throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2022            break;            break;
2023     case(Elements):     case(Elements):
2024     case(ReducedElements):     case(ReducedElements):
# Line 2123  int* MeshAdapter::borrowListOfTagsInUse( Line 2031  int* MeshAdapter::borrowListOfTagsInUse(
2031     case(Points):     case(Points):
2032            tags=mesh->Points->tagsInUse;            tags=mesh->Points->tagsInUse;
2033            break;            break;
    case(ContactElementsZero):  
    case(ReducedContactElementsZero):  
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
           tags=mesh->ContactElements->tagsInUse;  
           break;  
2034     default:     default:
2035        stringstream temp;        stringstream temp;
2036        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;        temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2037        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
2038    }    }
2039    return tags;    return tags;
2040  }  }
# Line 2147  bool MeshAdapter::canTag(int functionSpa Line 2049  bool MeshAdapter::canTag(int functionSpa
2049     case(FaceElements):     case(FaceElements):
2050     case(ReducedFaceElements):     case(ReducedFaceElements):
2051     case(Points):     case(Points):
    case(ContactElementsZero):  
    case(ReducedContactElementsZero):  
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
2052            return true;            return true;
2053     case(ReducedNodes):     case(ReducedNodes):
2054     case(DegreesOfFreedom):     case(DegreesOfFreedom):
# Line 2161  bool MeshAdapter::canTag(int functionSpa Line 2059  bool MeshAdapter::canTag(int functionSpa
2059    }    }
2060  }  }
2061    
2062    AbstractDomain::StatusType MeshAdapter::getStatus() const
2063    {
2064      Dudley_Mesh* mesh=m_dudleyMesh.get();
2065      return Dudley_Mesh_getStatus(mesh);
2066    }
2067    
2068    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2069    {
2070      
2071      Dudley_Mesh* mesh=m_dudleyMesh.get();
2072      int order =-1;
2073      switch(functionSpaceCode) {
2074       case(Nodes):
2075       case(DegreesOfFreedom):
2076              order=mesh->approximationOrder;
2077              break;
2078       case(ReducedNodes):
2079       case(ReducedDegreesOfFreedom):
2080              order=mesh->reducedApproximationOrder;
2081              break;
2082       case(Elements):
2083       case(FaceElements):
2084       case(Points):
2085              order=mesh->integrationOrder;
2086              break;
2087       case(ReducedElements):
2088       case(ReducedFaceElements):
2089              order=mesh->reducedIntegrationOrder;
2090              break;
2091       default:
2092          stringstream temp;
2093          temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2094          throw DudleyAdapterException(temp.str());
2095      }
2096      return order;
2097    }
2098    
2099    
2100    bool MeshAdapter::supportsContactElements() const
2101    {
2102        return false;
2103    }
2104    
2105  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.2385  
changed lines
  Added in v.3239

  ViewVC Help
Powered by ViewVC 1.1.26