/[escript]/trunk/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Diff of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

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

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

Legend:
Removed from v.1346  
changed lines
  Added in v.3344

  ViewVC Help
Powered by ViewVC 1.1.26