/[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 2487 by jfenwick, Tue Jun 23 04:26:07 2009 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    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2010 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 29  extern "C" { Line 29  extern "C" {
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;
 const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;  
 const int MeshAdapter::ReducedContactElementsZero=FINLEY_REDUCED_CONTACT_ELEMENTS_1;  
 const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;  
 const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;  
46    
47  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)  MeshAdapter::MeshAdapter(Dudley_Mesh* dudleyMesh)
48  {  {
49     setFunctionSpaceTypeNames();     setFunctionSpaceTypeNames();
50     //     //
51     // need to use a null_deleter as Finley_Mesh_free deletes the pointer     // need to use a null_deleter as Dudley_Mesh_free deletes the pointer
52     // for us.     // for us.
53     m_finleyMesh.reset(finleyMesh,null_deleter());     m_dudleyMesh.reset(dudleyMesh,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  }  }
# Line 70  MeshAdapter::~MeshAdapter() Line 66  MeshAdapter::~MeshAdapter()
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  void MeshAdapter::MPIBarrier() const
83  {  {
84  #ifdef PASO_MPI  #ifdef PASO_MPI
85     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);     MPI_Barrier(m_dudleyMesh.get()->MPIInfo->comm);
86  #endif  #endif
87     return;     return;
88  }  }
89  bool MeshAdapter::onMasterProcessor() const  bool MeshAdapter::onMasterProcessor() const
90  {  {
91     return m_finleyMesh.get()->MPIInfo->rank == 0;     return m_dudleyMesh.get()->MPIInfo->rank == 0;
92  }  }
93    
94    
95  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  #ifdef PASO_MPI
96     return m_finleyMesh.get();    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  void MeshAdapter::write(const std::string& fileName) const  
110    Dudley_Mesh* MeshAdapter::getDudley_Mesh() const {
111       return m_dudleyMesh.get();
112    }
113    
114    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[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
132     NcVar *ids;     NcVar *ids;
133     int *int_ptr;     int *int_ptr;
134     Finley_Mesh *mesh = m_finleyMesh.get();     Dudley_Mesh *mesh = m_dudleyMesh.get();
135     Finley_TagMap* tag_map;     Dudley_TagMap* tag_map;
136     int num_Tags = 0;     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;
# Line 129  void MeshAdapter::dump(const std::string Line 140  void MeshAdapter::dump(const std::string
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;
    int num_ContactElements_numNodes = mesh->ContactElements->numNodes;  
146  #ifdef PASO_MPI  #ifdef PASO_MPI
147     MPI_Status status;     MPI_Status status;
148  #endif  #endif
# Line 149  void MeshAdapter::dump(const std::string Line 158  void MeshAdapter::dump(const std::string
158     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
159     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
160     num_Tags = 0;     num_Tags = 0;
161     if (tag_map) {     while (tag_map) {
162        while (tag_map) {        num_Tags++;
163           num_Tags++;        tag_map=tag_map->next;
          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)");
    if (num_ContactElements>0)  
       if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )  
          throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);  
198     if (num_Tags>0)     if (num_Tags>0)
199        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
200           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Tags)");
201    
202     // Attributes: 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("numDim",numDim) )     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("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
212        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
213     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
214        throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(reduced_order)");
215     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
216        throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);        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) )     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
224        throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
225     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
226        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
227     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
228        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Elements_TypeId)");
229     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
230        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
231     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
232        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Points_TypeId)");
    if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )  
       throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);  
    if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )  
       throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);  
233     if (!dataFile.add_att("num_Tags", num_Tags) )     if (!dataFile.add_att("num_Tags", num_Tags) )
234        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Tags)");
235    
236     // // // // // Nodes // // // // //     // // // // // Nodes // // // // //
237    
238     // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)     // Nodes nodeDistribution
239       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
240          throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
241       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
242       if (! (ids->put(int_ptr, mpi_size+1)) )
243          throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
244    
245       // Nodes degreesOfFreedomDistribution
246       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
247          throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
248       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
249       if (! (ids->put(int_ptr, mpi_size+1)) )
250          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) {     if (numNodes>0) {
255    
256        // Nodes Id        // Nodes Id
257        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
258           throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Id)");
259        int_ptr = &mesh->Nodes->Id[0];        int_ptr = &mesh->Nodes->Id[0];
260        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
261           throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Id)");
262    
263        // Nodes Tag        // Nodes Tag
264        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
265           throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Tag)");
266        int_ptr = &mesh->Nodes->Tag[0];        int_ptr = &mesh->Nodes->Tag[0];
267        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
268           throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Tag)");
269    
270        // Nodes gDOF        // Nodes gDOF
271        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
272           throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
273        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
274        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
275           throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gDOF)");
276    
277        // Nodes global node index        // Nodes global node index
278        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
279           throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gNI)");
280        int_ptr = &mesh->Nodes->globalNodesIndex[0];        int_ptr = &mesh->Nodes->globalNodesIndex[0];
281        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
282           throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gNI)");
283    
284        // Nodes grDof        // Nodes grDof
285        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
286           throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
287        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
288        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
289           throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grDfI)");
290    
291        // Nodes grNI        // Nodes grNI
292        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
293           throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grNI)");
294        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
295        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
296           throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grNI)");
297    
298        // Nodes Coordinates        // Nodes Coordinates
299        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
300           throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
301        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
302           throw DataException("Error - MeshAdapter::dump: copy Nodes_Coordinates to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Coordinates)");
   
       // Nodes degreesOfFreedomDistribution  
       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )  
          throw DataException("Error - MeshAdapter::dump: appending Nodes_DofDistribution to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];  
       if (! (ids->put(int_ptr, mpi_size+1)) )  
          throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);  
   
       // Nodes nodeDistribution  
       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )  
          throw DataException("Error - MeshAdapter::dump: appending Nodes_NodeDistribution to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];  
       if (! (ids->put(int_ptr, mpi_size+1)) )  
          throw DataException("Error - MeshAdapter::dump: copy Nodes_NodeDistribution to netCDF buffer failed: " + *newFileName);  
303    
304     }     }
305    
# Line 311  void MeshAdapter::dump(const std::string Line 309  void MeshAdapter::dump(const std::string
309    
310        // Elements_Id        // Elements_Id
311        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
312           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Id)");
313        int_ptr = &mesh->Elements->Id[0];        int_ptr = &mesh->Elements->Id[0];
314        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
315           throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Id)");
316    
317        // Elements_Tag        // Elements_Tag
318        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
319           throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Tag)");
320        int_ptr = &mesh->Elements->Tag[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_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Tag)");
323    
324        // Elements_Owner        // Elements_Owner
325        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
326           throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Owner)");
327        int_ptr = &mesh->Elements->Owner[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_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Owner)");
330    
331        // Elements_Color        // Elements_Color
332        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
333           throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Color)");
334        int_ptr = &mesh->Elements->Color[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_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Color)");
337    
338        // Elements_Nodes        // Elements_Nodes
339        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
340           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Nodes)");
341        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
342           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Nodes)");
343    
344     }     }
345    
# Line 351  void MeshAdapter::dump(const std::string Line 349  void MeshAdapter::dump(const std::string
349    
350        // FaceElements_Id        // FaceElements_Id
351        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
352           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Id)");
353        int_ptr = &mesh->FaceElements->Id[0];        int_ptr = &mesh->FaceElements->Id[0];
354        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
355           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Id)");
356    
357        // FaceElements_Tag        // FaceElements_Tag
358        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
359           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
360        int_ptr = &mesh->FaceElements->Tag[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_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Tag)");
363    
364        // FaceElements_Owner        // FaceElements_Owner
365        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
366           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
367        int_ptr = &mesh->FaceElements->Owner[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_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Owner)");
370    
371        // FaceElements_Color        // FaceElements_Color
372        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
373           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Color)");
374        int_ptr = &mesh->FaceElements->Color[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_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Color)");
377    
378        // FaceElements_Nodes        // FaceElements_Nodes
379        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
380           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
381        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
382           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Nodes)");
   
    }  
   
    // // // // // Contact_Elements // // // // //  
   
    if (num_ContactElements>0) {  
   
       // ContactElements_Id  
       if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )  
          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->ContactElements->Id[0];  
       if (! (ids->put(int_ptr, num_ContactElements)) )  
          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);  
   
       // ContactElements_Tag  
       if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )  
          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->ContactElements->Tag[0];  
       if (! (ids->put(int_ptr, num_ContactElements)) )  
          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);  
   
       // ContactElements_Owner  
       if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )  
          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->ContactElements->Owner[0];  
       if (! (ids->put(int_ptr, num_ContactElements)) )  
          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);  
   
       // ContactElements_Color  
       if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )  
          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->ContactElements->Color[0];  
       if (! (ids->put(int_ptr, num_ContactElements)) )  
          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);  
   
       // ContactElements_Nodes  
       if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )  
          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);  
       if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )  
          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);  
383    
384     }     }
385    
# Line 433  void MeshAdapter::dump(const std::string Line 391  void MeshAdapter::dump(const std::string
391    
392        // Points_Id        // Points_Id
393        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
394           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Id)");
395        int_ptr = &mesh->Points->Id[0];        int_ptr = &mesh->Points->Id[0];
396        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
397           throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Id)");
398    
399        // Points_Tag        // Points_Tag
400        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
401           throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Tag)");
402        int_ptr = &mesh->Points->Tag[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_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Tag)");
405    
406        // Points_Owner        // Points_Owner
407        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
408           throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Owner)");
409        int_ptr = &mesh->Points->Owner[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_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Owner)");
412    
413        // Points_Color        // Points_Color
414        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
415           throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Color)");
416        int_ptr = &mesh->Points->Color[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_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Color)");
419    
420        // Points_Nodes        // Points_Nodes
421        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
422        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
423           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Nodes)");
424        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
425           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Nodes)");
426    
427     }     }
428    
# Line 488  void MeshAdapter::dump(const std::string Line 446  void MeshAdapter::dump(const std::string
446    
447        // Tags_keys        // Tags_keys
448        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
449           throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Tags_keys)");
450        int_ptr = &Tags_keys[0];        int_ptr = &Tags_keys[0];
451        if (! (ids->put(int_ptr, num_Tags)) )        if (! (ids->put(int_ptr, num_Tags)) )
452           throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Tags_keys)");
453    
454        // Tags_names_*        // Tags_names_*
455        // This is an array of strings, it should be stored as an array but instead I have hacked in one attribute per string        // This is an array of strings, it should be stored as an array but
456        // because the NetCDF manual doesn't tell how to do an array of strings        // instead I have hacked in one attribute per string because the NetCDF
457          // manual doesn't tell how to do an array of strings
458        tag_map = mesh->TagMap;        tag_map = mesh->TagMap;
459        if (tag_map) {        if (tag_map) {
460           int i = 0;           int i = 0;
461           while (tag_map) {           while (tag_map) {
462              sprintf(name_temp, "Tags_name_%d", i);              sprintf(name_temp, "Tags_name_%d", i);
463              if (!dataFile.add_att(name_temp, tag_map->name) )              if (!dataFile.add_att(name_temp, tag_map->name) )
464                 throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);                 throw DataException(msgPrefix+"add_att(Tags_names_XX)");
465              tag_map=tag_map->next;              tag_map=tag_map->next;
466              i++;              i++;
467           }           }
468        }        }
469    
470        TMPMEMFREE(Tags_keys);        TMPMEMFREE(Tags_keys);
   
471     }     }
472    
473  /* Send token to next MPI process so he can take his turn */  /* Send token to next MPI process so he can take his turn */
# Line 520  void MeshAdapter::dump(const std::string Line 478  void MeshAdapter::dump(const std::string
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
# Line 551  bool MeshAdapter::isValidFunctionSpaceTy Line 509  bool MeshAdapter::isValidFunctionSpaceTy
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"));
    m_functionSpaceTypeNames.insert  
    (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));  
    m_functionSpaceTypeNames.insert  
    (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));  
513     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
514     (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));     (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Dudley_ReducedDegreesOfFreedom"));
515     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
516     (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));     (FunctionSpaceNamesMapType::value_type(Nodes,"Dudley_Nodes"));
517     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
518     (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));     (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Dudley_Reduced_Nodes"));
519     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
520     (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));     (FunctionSpaceNamesMapType::value_type(Elements,"Dudley_Elements"));
521     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
522     (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));     (FunctionSpaceNamesMapType::value_type(ReducedElements,"Dudley_Reduced_Elements"));
523     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
524     (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));     (FunctionSpaceNamesMapType::value_type(FaceElements,"Dudley_Face_Elements"));
525     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
526     (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));     (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Dudley_Reduced_Face_Elements"));
527     m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
528     (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));     (FunctionSpaceNamesMapType::value_type(Points,"Dudley_Points"));
    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
# Line 607  int MeshAdapter::getReducedFunctionOnBou Line 557  int MeshAdapter::getReducedFunctionOnBou
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
# Line 643  int MeshAdapter::getDiracDeltaFunctionCo Line 595  int MeshAdapter::getDiracDeltaFunctionCo
595  //  //
596  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
597  {  {
598     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
599     int numDim=Finley_Mesh_getDim(mesh);     int numDim=Dudley_Mesh_getDim(mesh);
600     checkFinleyError();     checkDudleyError();
601     return numDim;     return numDim;
602  }  }
603    
# Line 654  int MeshAdapter::getDim() const Line 606  int MeshAdapter::getDim() const
606  //  //
607  int MeshAdapter::getNumDataPointsGlobal() const  int MeshAdapter::getNumDataPointsGlobal() const
608  {  {
609     return Finley_NodeFile_getGlobalNumNodes(m_finleyMesh.get()->Nodes);     return Dudley_NodeFile_getGlobalNumNodes(m_dudleyMesh.get()->Nodes);
610  }  }
611    
612  //  //
# Line 665  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;
# Line 705  pair<int,int> MeshAdapter::getDataShape( Line 657  pair<int,int> MeshAdapter::getDataShape(
657        numSamples=mesh->Points->numElements;        numSamples=mesh->Points->numElements;
658     }     }
659     break;     break;
    case(ContactElementsZero):  
    if (mesh->ContactElements!=NULL) {  
       numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;  
       numSamples=mesh->ContactElements->numElements;  
    }  
    break;  
    case(ReducedContactElementsZero):  
    if (mesh->ContactElements!=NULL) {  
       numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;  
       numSamples=mesh->ContactElements->numElements;  
    }  
    break;  
    case(ContactElementsOne):  
    if (mesh->ContactElements!=NULL) {  
       numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;  
       numSamples=mesh->ContactElements->numElements;  
    }  
    break;  
    case(ReducedContactElementsOne):  
    if (mesh->ContactElements!=NULL) {  
       numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;  
       numSamples=mesh->ContactElements->numElements;  
    }  
    break;  
660     case(DegreesOfFreedom):     case(DegreesOfFreedom):
661     if (mesh->Nodes!=NULL) {     if (mesh->Nodes!=NULL) {
662        numDataPointsPerSample=1;        numDataPointsPerSample=1;
663        numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);        numSamples=Dudley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
664     }     }
665     break;     break;
666     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
667     if (mesh->Nodes!=NULL) {     if (mesh->Nodes!=NULL) {
668        numDataPointsPerSample=1;        numDataPointsPerSample=1;
669        numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);        numSamples=Dudley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
670     }     }
671     break;     break;
672     default:     default:
673        stringstream temp;        stringstream temp;
674        temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
675        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
676        break;        break;
677     }     }
678     return pair<int,int>(numDataPointsPerSample,numSamples);     return pair<int,int>(numDataPointsPerSample,numSamples);
# Line 771  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(
# Line 792  void  MeshAdapter::addPDEToLumpedSystem( Line 719  void  MeshAdapter::addPDEToLumpedSystem(
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 806  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 814  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  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
# Line 848  void MeshAdapter::addPDEToTransportProbl Line 776  void MeshAdapter::addPDEToTransportProbl
776     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
777     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
778    
779     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
780     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tp.getPaso_TransportProblem();
781    
782     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );     Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
783     checkFinleyError();     checkDudleyError();
784    
785     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );     Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
786     checkFinleyError();     checkDudleyError();
787    
788     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );     Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
789     checkFinleyError();     checkDudleyError();
790    
791     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );     checkDudleyError();
    checkFinleyError();  
792  }  }
793    
794  //  //
# Line 872  void MeshAdapter::interpolateOnDomain(es Line 799  void MeshAdapter::interpolateOnDomain(es
799     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
800     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
801     if (inDomain!=*this)       if (inDomain!=*this)  
802        throw FinleyAdapterException("Error - Illegal domain of interpolant.");        throw DudleyAdapterException("Error - Illegal domain of interpolant.");
803     if (targetDomain!=*this)     if (targetDomain!=*this)
804        throw FinleyAdapterException("Error - Illegal domain of interpolation target.");        throw DudleyAdapterException("Error - Illegal domain of interpolation target.");
805    
806     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
807     escriptDataC _target=target.getDataC();     escriptDataC _target=target.getDataC();
808     escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
809     switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
# Line 886  void MeshAdapter::interpolateOnDomain(es Line 813  void MeshAdapter::interpolateOnDomain(es
813        case(ReducedNodes):        case(ReducedNodes):
814        case(DegreesOfFreedom):        case(DegreesOfFreedom):
815        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
816        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
817        break;        break;
818        case(Elements):        case(Elements):
819        case(ReducedElements):        case(ReducedElements):
820        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
821        break;        break;
822        case(FaceElements):        case(FaceElements):
823        case(ReducedFaceElements):        case(ReducedFaceElements):
824        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
825        break;        break;
826        case(Points):        case(Points):
827        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
       break;  
       case(ContactElementsZero):  
       case(ReducedContactElementsZero):  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
       break;  
       case(ContactElementsOne):  
       case(ReducedContactElementsOne):  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
828        break;        break;
829        default:        default:
830           stringstream temp;           stringstream temp;
831           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
832           throw FinleyAdapterException(temp.str());           throw DudleyAdapterException(temp.str());
833           break;           break;
834        }        }
835        break;        break;
# Line 920  void MeshAdapter::interpolateOnDomain(es Line 839  void MeshAdapter::interpolateOnDomain(es
839        case(ReducedNodes):        case(ReducedNodes):
840        case(DegreesOfFreedom):        case(DegreesOfFreedom):
841        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
842        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
843        break;        break;
844        case(Elements):        case(Elements):
845        case(ReducedElements):        case(ReducedElements):
846        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
847        break;        break;
848        case(FaceElements):        case(FaceElements):
849        case(ReducedFaceElements):        case(ReducedFaceElements):
850        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
851        break;        break;
852        case(Points):        case(Points):
853        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
       break;  
       case(ContactElementsZero):  
       case(ReducedContactElementsZero):  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
       break;  
       case(ContactElementsOne):  
       case(ReducedContactElementsOne):  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
854        break;        break;
855        default:        default:
856           stringstream temp;           stringstream temp;
857           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
858           throw FinleyAdapterException(temp.str());           throw DudleyAdapterException(temp.str());
859           break;           break;
860        }        }
861        break;        break;
862     case(Elements):     case(Elements):
863        if (target.getFunctionSpace().getTypeCode()==Elements) {        if (target.getFunctionSpace().getTypeCode()==Elements) {
864           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Dudley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
865        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
866           Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);           Dudley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
867        } else {        } else {
868           throw FinleyAdapterException("Error - No interpolation with data on elements possible.");           throw DudleyAdapterException("Error - No interpolation with data on elements possible.");
869        }        }
870        break;        break;
871     case(ReducedElements):     case(ReducedElements):
872        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
873           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Dudley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
874        } else {        } else {
875           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");           throw DudleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
876        }        }
877        break;        break;
878     case(FaceElements):     case(FaceElements):
879        if (target.getFunctionSpace().getTypeCode()==FaceElements) {        if (target.getFunctionSpace().getTypeCode()==FaceElements) {
880           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Dudley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
881        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
882           Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);           Dudley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
883        } else {        } else {
884           throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");           throw DudleyAdapterException("Error - No interpolation with data on face elements possible.");
885        }        }
886        break;        break;
887     case(ReducedFaceElements):     case(ReducedFaceElements):
888        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
889           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Dudley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
890        } else {        } else {
891           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");           throw DudleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
892        }        }
893        break;        break;
894     case(Points):     case(Points):
895        if (target.getFunctionSpace().getTypeCode()==Points) {        if (target.getFunctionSpace().getTypeCode()==Points) {
896           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);           Dudley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
       } else {  
          throw FinleyAdapterException("Error - No interpolation with data on points possible.");  
       }  
       break;  
    case(ContactElementsZero):  
    case(ContactElementsOne):  
       if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {  
          Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);  
       } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {  
          Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);  
897        } else {        } else {
898           throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");           throw DudleyAdapterException("Error - No interpolation with data on points possible.");
       }  
       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.");  
899        }        }
900        break;        break;
901     case(DegreesOfFreedom):           case(DegreesOfFreedom):      
902        switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
903        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
904        case(DegreesOfFreedom):        case(DegreesOfFreedom):
905        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
906        break;        break;
907        
908        case(Nodes):        case(Nodes):
# Line 1018  void MeshAdapter::interpolateOnDomain(es Line 911  void MeshAdapter::interpolateOnDomain(es
911           escript::Data temp=escript::Data(in);           escript::Data temp=escript::Data(in);
912           temp.expand();           temp.expand();
913           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
914           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);           Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
915        } else {        } else {
916           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);           Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
917        }        }
918        break;        break;
919        case(Elements):        case(Elements):
# Line 1028  void MeshAdapter::interpolateOnDomain(es Line 921  void MeshAdapter::interpolateOnDomain(es
921        if (getMPISize()>1) {        if (getMPISize()>1) {
922           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
923           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
924           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
925        } else {        } else {
926           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
927        }        }
928        break;        break;
929        case(FaceElements):        case(FaceElements):
# Line 1038  void MeshAdapter::interpolateOnDomain(es Line 931  void MeshAdapter::interpolateOnDomain(es
931        if (getMPISize()>1) {        if (getMPISize()>1) {
932           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
933           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
934           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
935        
936        } else {        } else {
937           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
938        }        }
939        break;        break;
940        case(Points):        case(Points):
# Line 1049  void MeshAdapter::interpolateOnDomain(es Line 942  void MeshAdapter::interpolateOnDomain(es
942           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
943           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
944        } else {        } else {
945           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
       }  
       break;  
       case(ContactElementsZero):  
       case(ContactElementsOne):  
       case(ReducedContactElementsZero):  
       case(ReducedContactElementsOne):  
       if (getMPISize()>1) {  
          escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
          escriptDataC _in2 = temp.getDataC();  
          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);  
       } else {  
          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
946        }        }
947        break;        break;
948        default:        default:
949           stringstream temp;           stringstream temp;
950           temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
951           throw FinleyAdapterException(temp.str());           throw DudleyAdapterException(temp.str());
952           break;           break;
953        }        }
954        break;        break;
955     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
956        switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
957        case(Nodes):        case(Nodes):
958        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");        throw DudleyAdapterException("Error - Dudley does not support interpolation from reduced degrees of freedom to mesh nodes.");
959        break;        break;
960        case(ReducedNodes):        case(ReducedNodes):
961        if (getMPISize()>1) {        if (getMPISize()>1) {
962           escript::Data temp=escript::Data(in);           escript::Data temp=escript::Data(in);
963           temp.expand();           temp.expand();
964           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
965           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);           Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
966        } else {        } else {
967           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);           Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
968        }        }
969        break;        break;
970        case(DegreesOfFreedom):        case(DegreesOfFreedom):
971        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");        throw DudleyAdapterException("Error - Dudley does not support interpolation from reduced degrees of freedom to degrees of freedom");
972        break;        break;
973        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
974        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
975        break;        break;
976        case(Elements):        case(Elements):
977        case(ReducedElements):        case(ReducedElements):
978        if (getMPISize()>1) {        if (getMPISize()>1) {
979           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
980           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
981           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
982        } else {        } else {
983           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
984        }        }
985        break;        break;
986        case(FaceElements):        case(FaceElements):
# Line 1107  void MeshAdapter::interpolateOnDomain(es Line 988  void MeshAdapter::interpolateOnDomain(es
988        if (getMPISize()>1) {        if (getMPISize()>1) {
989           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
990           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
991           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
992        } else {        } else {
993           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
994        }        }
995        break;        break;
996        case(Points):        case(Points):
997        if (getMPISize()>1) {        if (getMPISize()>1) {
998           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
999           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1000           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1001        } else {        } else {
1002           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
       }  
       break;  
       case(ContactElementsZero):  
       case(ContactElementsOne):  
       case(ReducedContactElementsZero):  
       case(ReducedContactElementsOne):  
       if (getMPISize()>1) {  
          escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );  
          escriptDataC _in2 = temp.getDataC();  
          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);  
       } else {  
          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
1003        }        }
1004        break;        break;
1005        default:        default:
1006           stringstream temp;           stringstream temp;
1007           temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1008           throw FinleyAdapterException(temp.str());           throw DudleyAdapterException(temp.str());
1009           break;           break;
1010        }        }
1011        break;        break;
1012     default:     default:
1013        stringstream temp;        stringstream temp;
1014        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();        temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1015        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1016        break;        break;
1017     }     }
1018     checkFinleyError();     checkDudleyError();
1019  }  }
1020    
1021  //  //
# Line 1156  void MeshAdapter::setToX(escript::Data& Line 1025  void MeshAdapter::setToX(escript::Data&
1025  {  {
1026     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1027     if (argDomain!=*this)     if (argDomain!=*this)
1028        throw FinleyAdapterException("Error - Illegal domain of data point locations");        throw DudleyAdapterException("Error - Illegal domain of data point locations");
1029     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1030     // in case of values node coordinates we can do the job directly:     // in case of values node coordinates we can do the job directly:
1031     if (arg.getFunctionSpace().getTypeCode()==Nodes) {     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1032        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1033        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Dudley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1034     } else {     } else {
1035        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
1036        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1037        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Dudley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1038        // this is then interpolated onto arg:        // this is then interpolated onto arg:
1039        interpolateOnDomain(arg,tmp_data);        interpolateOnDomain(arg,tmp_data);
1040     }     }
1041     checkFinleyError();     checkDudleyError();
1042  }  }
1043    
1044  //  //
# Line 1180  void MeshAdapter::setToNormal(escript::D Line 1049  void MeshAdapter::setToNormal(escript::D
1049  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1050     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1051     if (normalDomain!=*this)     if (normalDomain!=*this)
1052        throw FinleyAdapterException("Error - Illegal domain of normal locations");        throw DudleyAdapterException("Error - Illegal domain of normal locations");
1053     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1054     escriptDataC _normal=normal.getDataC();     escriptDataC _normal=normal.getDataC();
1055     switch(normal.getFunctionSpace().getTypeCode()) {     switch(normal.getFunctionSpace().getTypeCode()) {
1056     case(Nodes):     case(Nodes):
1057     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for nodes");
1058     break;     break;
1059     case(ReducedNodes):     case(ReducedNodes):
1060     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for reduced nodes");
1061     break;     break;
1062     case(Elements):     case(Elements):
1063     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for elements");
1064     break;     break;
1065     case(ReducedElements):     case(ReducedElements):
1066     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for elements with reduced integration order");
1067     break;     break;
1068     case (FaceElements):     case (FaceElements):
1069     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     Dudley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1070     break;     break;
1071     case (ReducedFaceElements):     case (ReducedFaceElements):
1072     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     Dudley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1073     break;     break;
1074     case(Points):     case(Points):
1075     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for point elements");
    break;  
    case (ContactElementsOne):  
    case (ContactElementsZero):  
    Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);  
    break;  
    case (ReducedContactElementsOne):  
    case (ReducedContactElementsZero):  
    Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);  
1076     break;     break;
1077     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1078     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for degrees of freedom.");
1079     break;     break;
1080     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1081     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");     throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for reduced degrees of freedom.");
1082     break;     break;
1083     default:     default:
1084        stringstream temp;        stringstream temp;
1085        temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();        temp << "Error - Normal Vectors: Dudley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1086        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1087        break;        break;
1088     }     }
1089     checkFinleyError();     checkDudleyError();
1090  }  }
1091    
1092  //  //
# Line 1236  void MeshAdapter::interpolateACross(escr Line 1097  void MeshAdapter::interpolateACross(escr
1097     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1098     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1099     if (targetDomain!=this)     if (targetDomain!=this)
1100        throw FinleyAdapterException("Error - Illegal domain of interpolation target");        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();
# Line 1259  void MeshAdapter::setToIntegrals(std::ve Line 1120  void MeshAdapter::setToIntegrals(std::ve
1120     case(Nodes):     case(Nodes):
1121     temp=escript::Data( arg, escript::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, escript::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.");
    break;  
    case(ContactElementsZero):  
    Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);  
    break;  
    case(ReducedContactElementsZero):  
    Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);  
    break;  
    case(ContactElementsOne):  
    Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);  
    break;  
    case(ReducedContactElementsOne):  
    Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);  
1144     break;     break;
1145     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1146     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1147     _temp=temp.getDataC();     _temp=temp.getDataC();
1148     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1149     break;     break;
1150     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1151     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1152     _temp=temp.getDataC();     _temp=temp.getDataC();
1153     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1154     break;     break;
1155     default:     default:
1156        stringstream temp;        stringstream temp;
1157        temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();        temp << "Error - Integrals: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1158        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1159        break;        break;
1160     }     }
1161     checkFinleyError();     checkDudleyError();
1162     blocktimer_increment("integrate()", blocktimer_start);     blocktimer_increment("integrate()", blocktimer_start);
1163  }  }
1164    
# Line 1320  void MeshAdapter::setToGradient(escript: Line 1169  void MeshAdapter::setToGradient(escript:
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;
# Line 1344  void MeshAdapter::setToGradient(escript: Line 1193  void MeshAdapter::setToGradient(escript:
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.");
    break;  
    case(ContactElementsZero):  
    Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);  
    break;  
    case(ReducedContactElementsZero):  
    Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);  
    break;  
    case(ContactElementsOne):  
    Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);  
    break;  
    case(ReducedContactElementsOne):  
    Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);  
1215     break;     break;
1216     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1217     throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");     throw DudleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1218     break;     break;
1219     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1220     throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");     throw DudleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1221     break;     break;
1222     default:     default:
1223        stringstream temp;        stringstream temp;
1224        temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();        temp << "Error - Gradient: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1225        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1226        break;        break;
1227     }     }
1228     checkFinleyError();     checkDudleyError();
1229  }  }
1230    
1231  //  //
# Line 1396  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.");
    break;  
    case(ContactElementsZero):  
    case(ContactElementsOne):  
    Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);  
    break;  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);  
1259     break;     break;
1260     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1261     throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");     throw DudleyAdapterException("Error - Size of degrees of freedom is not supported.");
1262     break;     break;
1263     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1264     throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");     throw DudleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1265     break;     break;
1266     default:     default:
1267        stringstream temp;        stringstream temp;
1268        temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();        temp << "Error - Element size: Dudley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1269        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1270        break;        break;
1271     }     }
1272     checkFinleyError();     checkDudleyError();
1273  }  }
1274    
1275  //  //
# Line 1448  void MeshAdapter::setToSize(escript::Dat Line 1277  void MeshAdapter::setToSize(escript::Dat
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           tmp = new_x_inter.getDataC();
1291           Dudley_Mesh_setCoordinates(mesh,&tmp);
1292       }
1293       checkDudleyError();
1294  }  }
1295    
1296  //  //
# Line 1472  void MeshAdapter::extractArgsFromDict(co Line 1307  void MeshAdapter::extractArgsFromDict(co
1307    
1308     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1309     for (int i=0; i<numData; ++i) {     for (int i=0; i<numData; ++i) {
1310        std::string n=boost::python::extract<std::string>(keys[i]);        string n=boost::python::extract<string>(keys[i]);
1311        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1312        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1313           throw FinleyAdapterException("Error: Data must be defined on same Domain");           throw DudleyAdapterException("Error: Data must be defined on same Domain");
1314        data[i] = d.getDataC();        data[i] = d.getDataC();
1315        dataPtr[i] = &(data[i]);        dataPtr[i] = &(data[i]);
1316        names[i] = TMPMEMALLOC(n.length()+1, char);        names[i] = TMPMEMALLOC(n.length()+1, char);
# Line 1486  void MeshAdapter::extractArgsFromDict(co Line 1321  void MeshAdapter::extractArgsFromDict(co
1321  //  //
1322  // saves mesh and optionally data arrays in openDX format  // saves mesh and optionally data arrays in openDX format
1323  //  //
1324  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1325  {  {
1326     int num_data;     int num_data;
1327     char **names;     char **names;
# Line 1494  void MeshAdapter::saveDX(const std::stri Line 1329  void MeshAdapter::saveDX(const std::stri
1329     escriptDataC **ptr_data;     escriptDataC **ptr_data;
1330    
1331     extractArgsFromDict(arg, num_data, names, data, ptr_data);     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1332     Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);     Dudley_Mesh_saveDX(filename.c_str(), m_dudleyMesh.get(), num_data, names, ptr_data);
1333     checkFinleyError();     checkDudleyError();
1334    
1335     /* win32 refactor */     /* win32 refactor */
1336     TMPMEMFREE(data);     TMPMEMFREE(data);
# Line 1512  void MeshAdapter::saveDX(const std::stri Line 1347  void MeshAdapter::saveDX(const std::stri
1347  //  //
1348  // saves mesh and optionally data arrays in VTK format  // saves mesh and optionally data arrays in VTK format
1349  //  //
1350  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg,  const std::string& metadata, const std::string& metadata_schema) const  void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1351  {  {
1352     int num_data;     int num_data;
1353     char **names;     char **names;
# Line 1520  void MeshAdapter::saveVTK(const std::str Line 1355  void MeshAdapter::saveVTK(const std::str
1355     escriptDataC **ptr_data;     escriptDataC **ptr_data;
1356    
1357     extractArgsFromDict(arg, num_data, names, data, ptr_data);     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1358     Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());     Dudley_Mesh_saveVTK(filename.c_str(), m_dudleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());
1359     checkFinleyError();     checkDudleyError();
1360    
1361     /* win32 refactor */     /* win32 refactor */
1362     TMPMEMFREE(data);     TMPMEMFREE(data);
# Line 1533  void MeshAdapter::saveVTK(const std::str Line 1368  void MeshAdapter::saveVTK(const std::str
1368     TMPMEMFREE(names);     TMPMEMFREE(names);
1369  }  }
1370    
1371    bool MeshAdapter::ownSample(int fs_code, index_t id) const
1372    {
1373    #ifdef PASO_MPI
1374        index_t myFirstNode=0, myLastNode=0, k=0;
1375        index_t* globalNodeIndex=0;
1376        Dudley_Mesh* mesh_p=m_dudleyMesh.get();
1377        if (fs_code == DUDLEY_REDUCED_NODES)
1378        {
1379        myFirstNode = Dudley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1380        myLastNode = Dudley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1381        globalNodeIndex = Dudley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1382        }
1383        else
1384        {
1385        myFirstNode = Dudley_NodeFile_getFirstNode(mesh_p->Nodes);
1386        myLastNode = Dudley_NodeFile_getLastNode(mesh_p->Nodes);
1387        globalNodeIndex = Dudley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1388        }
1389        k=globalNodeIndex[id];
1390        return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1391    #endif
1392        return true;
1393    }
1394    
1395    
1396    
1397  //  //
1398  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1399  //  //
# Line 1548  SystemMatrixAdapter MeshAdapter::newSyst Line 1409  SystemMatrixAdapter MeshAdapter::newSyst
1409     // is the domain right?     // is the domain right?
1410     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1411     if (row_domain!=*this)     if (row_domain!=*this)
1412        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");        throw DudleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1413     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1414     if (col_domain!=*this)     if (col_domain!=*this)
1415        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        throw DudleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1416     // is the function space type right     // is the function space type right
1417     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1418        reduceRowOrder=0;        reduceRowOrder=0;
1419     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1420        reduceRowOrder=1;        reduceRowOrder=1;
1421     } else {     } else {
1422        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");        throw DudleyAdapterException("Error - illegal function space type for system matrix rows.");
1423     }     }
1424     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1425        reduceColOrder=0;        reduceColOrder=0;
1426     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1427        reduceColOrder=1;        reduceColOrder=1;
1428     } else {     } else {
1429        throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");        throw DudleyAdapterException("Error - illegal function space type for system matrix columns.");
1430     }     }
1431     // generate matrix:     // generate matrix:
1432    
1433     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Dudley_getPattern(getDudley_Mesh(),reduceRowOrder,reduceColOrder);
1434     checkFinleyError();     checkDudleyError();
1435     Paso_SystemMatrix* fsystemMatrix;     Paso_SystemMatrix* fsystemMatrix;
1436     int trilinos = 0;     int trilinos = 0;
1437     if (trilinos) {     if (trilinos) {
# Line 1579  SystemMatrixAdapter MeshAdapter::newSyst Line 1440  SystemMatrixAdapter MeshAdapter::newSyst
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);
# Line 1590  SystemMatrixAdapter MeshAdapter::newSyst Line 1451  SystemMatrixAdapter MeshAdapter::newSyst
1451  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1452  //  //
1453  TransportProblemAdapter MeshAdapter::newTransportProblem(  TransportProblemAdapter MeshAdapter::newTransportProblem(
1454                                                           const double theta,                                                           const bool useBackwardEuler,
1455                                                           const int blocksize,                                                           const int blocksize,
1456                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1457                                                           const int type) const                                                           const int type) const
# Line 1599  TransportProblemAdapter MeshAdapter::new Line 1460  TransportProblemAdapter MeshAdapter::new
1460     // is the domain right?     // is the domain right?
1461     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1462     if (domain!=*this)     if (domain!=*this)
1463        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");        throw DudleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1464     // is the function space type right     // is the function space type right
1465     if (functionspace.getTypeCode()==DegreesOfFreedom) {     if (functionspace.getTypeCode()==DegreesOfFreedom) {
1466        reduceOrder=0;        reduceOrder=0;
1467     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1468        reduceOrder=1;        reduceOrder=1;
1469     } else {     } else {
1470        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");        throw DudleyAdapterException("Error - illegal function space type for system matrix rows.");
1471     }     }
1472     // generate matrix:     // generate matrix:
1473    
1474     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Dudley_getPattern(getDudley_Mesh(),reduceOrder,reduceOrder);
1475     checkFinleyError();     checkDudleyError();
1476     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1477     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1478     checkPasoError();     checkPasoError();
1479     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1480     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1481  }  }
1482    
1483  //  //
# Line 1637  bool MeshAdapter::isCellOriented(int fun Line 1498  bool MeshAdapter::isCellOriented(int fun
1498     case(Elements):     case(Elements):
1499     case(FaceElements):     case(FaceElements):
1500     case(Points):     case(Points):
    case(ContactElementsZero):  
    case(ContactElementsOne):  
1501     case(ReducedElements):     case(ReducedElements):
1502     case(ReducedFaceElements):     case(ReducedFaceElements):
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
1503     return true;     return true;
1504     break;     break;
1505     default:     default:
1506        stringstream temp;        stringstream temp;
1507        temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;        temp << "Error - Cell: Dudley does not know anything about function space type " << functionSpaceCode;
1508        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1509        break;        break;
1510     }     }
1511     checkFinleyError();     checkDudleyError();
1512     return false;     return false;
1513  }  }
1514    
1515    bool
1516    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());
    default:  
       stringstream temp;  
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;  
       throw FinleyAdapterException(temp.str());  
1659     }     }
1660     break;     break;
1661     case(ReducedNodes):     case(ReducedNodes):
1662     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1663     case(ReducedNodes):      case(ReducedNodes):
1664     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1665     case(Elements):      case(Elements):
1666     case(ReducedElements):      case(ReducedElements):
1667     case(FaceElements):      case(FaceElements):
1668     case(ReducedFaceElements):      case(ReducedFaceElements):
1669     case(Points):      case(Points):
1670     case(ContactElementsZero):      return true;
1671     case(ReducedContactElementsZero):      case(Nodes):
1672     case(ContactElementsOne):      case(DegreesOfFreedom):
1673     case(ReducedContactElementsOne):      return false;
1674     return true;      default:
1675     case(Nodes):          stringstream temp;
1676     case(DegreesOfFreedom):          temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1677     return false;          throw DudleyAdapterException(temp.str());
    default:  
       stringstream temp;  
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;  
       throw FinleyAdapterException(temp.str());  
1678     }     }
1679     break;     break;
1680     case(Elements):     case(Elements):
1681     if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1682        return true;        return true;
1683     } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1684        return true;        return true;
1685     } else {          } else {
1686        return false;            return false;
1687     }          }
1688     case(ReducedElements):     case(ReducedElements):
1689     if (functionSpaceType_target==ReducedElements) {      if (functionSpaceType_target==ReducedElements) {
1690        return true;        return true;
1691     } else {      } else {
1692        return false;            return false;
1693     }      }
1694     case(FaceElements):     case(FaceElements):
1695     if (functionSpaceType_target==FaceElements) {      if (functionSpaceType_target==FaceElements) {
1696        return true;              return true;
1697     } else if (functionSpaceType_target==ReducedFaceElements) {      } else if (functionSpaceType_target==ReducedFaceElements) {
1698        return true;              return true;
1699     } else {      } else {
1700        return false;              return false;
1701     }      }
1702     case(ReducedFaceElements):     case(ReducedFaceElements):
1703     if (functionSpaceType_target==ReducedFaceElements) {      if (functionSpaceType_target==ReducedFaceElements) {
1704        return true;              return true;
1705     } else {      } else {
1706        return false;          return false;
1707     }      }
1708     case(Points):     case(Points):
1709     if (functionSpaceType_target==Points) {      if (functionSpaceType_target==Points) {
1710        return true;              return true;
1711     } else {      } else {
1712        return false;              return false;
1713     }      }
    case(ContactElementsZero):  
    case(ContactElementsOne):  
    if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {  
       return true;  
    } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {  
       return true;  
    } else {  
       return false;  
    }  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {  
       return true;  
    } else {  
       return false;  
    }  
1714     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1715     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1716     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1717     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1718     case(Nodes):      case(Nodes):
1719     case(ReducedNodes):      case(ReducedNodes):
1720     case(Elements):      case(Elements):
1721     case(ReducedElements):      case(ReducedElements):
1722     case(Points):      case(Points):
1723     case(FaceElements):      case(FaceElements):
1724     case(ReducedFaceElements):      case(ReducedFaceElements):
1725     case(ContactElementsZero):      return true;
1726     case(ReducedContactElementsZero):      default:
1727     case(ContactElementsOne):          stringstream temp;
1728     case(ReducedContactElementsOne):          temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1729     return true;          throw DudleyAdapterException(temp.str());
1730     default:      }
1731        stringstream temp;      break;
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;  
       throw FinleyAdapterException(temp.str());  
    }  
    break;  
1732     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1733     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
1734     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1735     case(ReducedNodes):      case(ReducedNodes):
1736     case(Elements):      case(Elements):
1737     case(ReducedElements):      case(ReducedElements):
1738     case(FaceElements):      case(FaceElements):
1739     case(ReducedFaceElements):      case(ReducedFaceElements):
1740     case(Points):      case(Points):
1741     case(ContactElementsZero):      return true;
1742     case(ReducedContactElementsZero):      case(Nodes):
1743     case(ContactElementsOne):      case(DegreesOfFreedom):
1744     case(ReducedContactElementsOne):      return false;
1745     return true;      default:
1746     case(Nodes):          stringstream temp;
1747     case(DegreesOfFreedom):          temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1748     return false;          throw DudleyAdapterException(temp.str());
1749     default:      }
1750        stringstream temp;      break;
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;  
       throw FinleyAdapterException(temp.str());  
    }  
    break;  
1751     default:     default:
1752        stringstream temp;        stringstream temp;
1753        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;        temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_source;
1754        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1755        break;        break;
1756     }     }
1757     checkFinleyError();     checkDudleyError();
1758     return false;     return false;
1759  }  }
1760    
# Line 1817  bool MeshAdapter::operator==(const Abstr Line 1767  bool MeshAdapter::operator==(const Abstr
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     }     }
# Line 1830  bool MeshAdapter::operator!=(const Abstr Line 1780  bool MeshAdapter::operator!=(const Abstr
1780    
1781  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1782  {  {
1783     Finley_Mesh* mesh=m_finleyMesh.get();     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);     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1785     checkPasoError();     checkPasoError();
1786     return out;     return out;
1787  }  }
1788  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1789  {  {
1790     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1791     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1792     checkPasoError();     checkPasoError();
1793     return out;     return out;
1794  }  }
# Line 1861  escript::Data MeshAdapter::getSize() con Line 1811  escript::Data MeshAdapter::getSize() con
1811  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1812  {  {
1813     int *out = NULL;     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;
# Line 1884  const int* MeshAdapter::borrowSampleRefe Line 1834  const int* MeshAdapter::borrowSampleRefe
1834     case(Points):     case(Points):
1835     out=mesh->Points->Id;     out=mesh->Points->Id;
1836     break;     break;
    case(ContactElementsZero):  
    out=mesh->ContactElements->Id;  
    break;  
    case(ReducedContactElementsZero):  
    out=mesh->ContactElements->Id;  
    break;  
    case(ContactElementsOne):  
    out=mesh->ContactElements->Id;  
    break;  
    case(ReducedContactElementsOne):  
    out=mesh->ContactElements->Id;  
    break;  
1837     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1838     out=mesh->Nodes->degreesOfFreedomId;     out=mesh->Nodes->degreesOfFreedomId;
1839     break;     break;
# Line 1905  const int* MeshAdapter::borrowSampleRefe Line 1843  const int* MeshAdapter::borrowSampleRefe
1843     default:     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;
# Line 1913  const int* MeshAdapter::borrowSampleRefe Line 1851  const int* MeshAdapter::borrowSampleRefe
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];
# Line 1936  int MeshAdapter::getTagFromSampleNo(int Line 1874  int MeshAdapter::getTagFromSampleNo(int
1874     case(Points):     case(Points):
1875     out=mesh->Points->Tag[sampleNo];     out=mesh->Points->Tag[sampleNo];
1876     break;     break;
    case(ContactElementsZero):  
    out=mesh->ContactElements->Tag[sampleNo];  
    break;  
    case(ReducedContactElementsZero):  
    out=mesh->ContactElements->Tag[sampleNo];  
    break;  
    case(ContactElementsOne):  
    out=mesh->ContactElements->Tag[sampleNo];  
    break;  
    case(ReducedContactElementsOne):  
    out=mesh->ContactElements->Tag[sampleNo];  
    break;  
1877     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1878     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");     throw DudleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1879     break;     break;
1880     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1881     throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");     throw DudleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1882     break;     break;
1883     default:     default:
1884        stringstream temp;        stringstream temp;
1885        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1886        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1887        break;        break;
1888     }     }
1889     return out;     return out;
# Line 1966  int MeshAdapter::getTagFromSampleNo(int Line 1892  int MeshAdapter::getTagFromSampleNo(int
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;     break;
1901     case(ReducedNodes):     case(ReducedNodes):
1902     throw FinleyAdapterException("Error - ReducedNodes does not support tags");     throw DudleyAdapterException("Error - ReducedNodes does not support tags");
1903     break;     break;
1904     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1905     throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");     throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
1906     break;     break;
1907     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1908     throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");     throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1909     break;     break;
1910     case(Elements):     case(Elements):
1911     Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     Dudley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1912     break;     break;
1913     case(ReducedElements):     case(ReducedElements):
1914     Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     Dudley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1915     break;     break;
1916     case(FaceElements):     case(FaceElements):
1917     Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);     Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1918     break;     break;
1919     case(ReducedFaceElements):     case(ReducedFaceElements):
1920     Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);     Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1921     break;     break;
1922     case(Points):     case(Points):
1923     Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);     Dudley_ElementFile_setTags(mesh->Points,newTag,&tmp);
    break;  
    case(ContactElementsZero):  
    Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);  
    break;  
    case(ReducedContactElementsZero):  
    Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);  
    break;  
    case(ContactElementsOne):  
    Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);  
    break;  
    case(ReducedContactElementsOne):  
    Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);  
1924     break;     break;
1925     default:     default:
1926        stringstream temp;        stringstream temp;
1927        temp << "Error - Finley does not know anything about function space type " << functionSpaceType;        temp << "Error - Dudley does not know anything about function space type " << functionSpaceType;
1928        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
1929     }     }
1930     checkFinleyError();     checkDudleyError();
1931     return;     return;
1932  }  }
1933    
1934  void MeshAdapter::setTagMap(const std::string& name,  int tag)  void MeshAdapter::setTagMap(const string& name,  int tag)
1935  {  {
1936     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1937     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Dudley_Mesh_addTagMap(mesh, name.c_str(),tag);
1938     checkPasoError();     checkPasoError();
1939     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1940  }  }
1941    
1942  int MeshAdapter::getTag(const std::string& name) const  int MeshAdapter::getTag(const string& name) const
1943  {  {
1944     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1945     int tag=0;     int tag=0;
1946     tag=Finley_Mesh_getTag(mesh, name.c_str());     tag=Dudley_Mesh_getTag(mesh, name.c_str());
1947     checkPasoError();     checkPasoError();
1948     // throwStandardException("MeshAdapter::getTag is not implemented.");     // throwStandardException("MeshAdapter::getTag is not implemented.");
1949     return tag;     return tag;
1950  }  }
1951    
1952  bool MeshAdapter::isValidTagName(const std::string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
1953  {  {
1954     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1955     return Finley_Mesh_isValidTagName(mesh,name.c_str());     return Dudley_Mesh_isValidTagName(mesh,name.c_str());
1956  }  }
1957    
1958  std::string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
1959  {  {
1960     stringstream temp;     stringstream temp;
1961     Finley_Mesh* mesh=m_finleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1962     Finley_TagMap* tag_map=mesh->TagMap;     Dudley_TagMap* tag_map=mesh->TagMap;
1963     while (tag_map) {     while (tag_map) {
1964        temp << tag_map->name;        temp << tag_map->name;
1965        tag_map=tag_map->next;        tag_map=tag_map->next;
# Line 2056  std::string MeshAdapter::showTagNames() Line 1970  std::string MeshAdapter::showTagNames()
1970    
1971  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
1972  {  {
1973    Finley_Mesh* mesh=m_finleyMesh.get();    Dudley_Mesh* mesh=m_dudleyMesh.get();
1974    dim_t numTags=0;    dim_t numTags=0;
1975    switch(functionSpaceCode) {    switch(functionSpaceCode) {
1976     case(Nodes):     case(Nodes):
1977            numTags=mesh->Nodes->numTagsInUse;            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     case(ReducedElements):     case(ReducedElements):
# Line 2082  int MeshAdapter::getNumberOfTagsInUse(in Line 1996  int MeshAdapter::getNumberOfTagsInUse(in
1996     case(Points):     case(Points):
1997            numTags=mesh->Points->numTagsInUse;            numTags=mesh->Points->numTagsInUse;
1998            break;            break;
    case(ContactElementsZero):  
    case(ReducedContactElementsZero):  
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
           numTags=mesh->ContactElements->numTagsInUse;  
           break;  
1999     default:     default:
2000        stringstream temp;        stringstream temp;
2001        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;        temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2002        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
2003    }    }
2004    return numTags;    return numTags;
2005  }  }
2006    
2007  const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2008  {  {
2009    Finley_Mesh* mesh=m_finleyMesh.get();    Dudley_Mesh* mesh=m_dudleyMesh.get();
2010    index_t* tags=NULL;    index_t* tags=NULL;
2011    switch(functionSpaceCode) {    switch(functionSpaceCode) {
2012     case(Nodes):     case(Nodes):
2013            tags=mesh->Nodes->tagsInUse;            tags=mesh->Nodes->tagsInUse;
2014            break;            break;
2015     case(ReducedNodes):     case(ReducedNodes):
2016            throw FinleyAdapterException("Error - ReducedNodes does not support tags");            throw DudleyAdapterException("Error - ReducedNodes does not support tags");
2017            break;            break;
2018     case(DegreesOfFreedom):     case(DegreesOfFreedom):
2019            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");            throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
2020            break;            break;
2021     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2022            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");            throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2023            break;            break;
2024     case(Elements):     case(Elements):
2025     case(ReducedElements):     case(ReducedElements):
# Line 2124  const int* MeshAdapter::borrowListOfTags Line 2032  const int* MeshAdapter::borrowListOfTags
2032     case(Points):     case(Points):
2033            tags=mesh->Points->tagsInUse;            tags=mesh->Points->tagsInUse;
2034            break;            break;
    case(ContactElementsZero):  
    case(ReducedContactElementsZero):  
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
           tags=mesh->ContactElements->tagsInUse;  
           break;  
2035     default:     default:
2036        stringstream temp;        stringstream temp;
2037        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;        temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2038        throw FinleyAdapterException(temp.str());        throw DudleyAdapterException(temp.str());
2039    }    }
2040    return tags;    return tags;
2041  }  }
# Line 2148  bool MeshAdapter::canTag(int functionSpa Line 2050  bool MeshAdapter::canTag(int functionSpa
2050     case(FaceElements):     case(FaceElements):
2051     case(ReducedFaceElements):     case(ReducedFaceElements):
2052     case(Points):     case(Points):
    case(ContactElementsZero):  
    case(ReducedContactElementsZero):  
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
2053            return true;            return true;
2054     case(ReducedNodes):     case(ReducedNodes):
2055     case(DegreesOfFreedom):     case(DegreesOfFreedom):
# Line 2162  bool MeshAdapter::canTag(int functionSpa Line 2060  bool MeshAdapter::canTag(int functionSpa
2060    }    }
2061  }  }
2062    
2063    AbstractDomain::StatusType MeshAdapter::getStatus() const
2064    {
2065      Dudley_Mesh* mesh=m_dudleyMesh.get();
2066      return Dudley_Mesh_getStatus(mesh);
2067    }
2068    
2069    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2070    {
2071      
2072      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    
2101    bool MeshAdapter::supportsContactElements() const
2102    {
2103        return false;
2104    }
2105    
2106    ReferenceElementSetWrapper::ReferenceElementSetWrapper(ElementTypeId id, index_t order, index_t reducedOrder)
2107    {
2108      m_refSet = Dudley_ReferenceElementSet_alloc(id, order, reducedOrder);
2109    }
2110    
2111    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2112    {
2113      Dudley_ReferenceElementSet_dealloc(m_refSet);
2114    }
2115    
2116  }  // end of namespace  }  // end of namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26