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

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

  ViewVC Help
Powered by ViewVC 1.1.26