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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26