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

Legend:
Removed from v.1375  
changed lines
  Added in v.3242

  ViewVC Help
Powered by ViewVC 1.1.26