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

Legend:
Removed from v.1343  
changed lines
  Added in v.3090

  ViewVC Help
Powered by ViewVC 1.1.26