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

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

  ViewVC Help
Powered by ViewVC 1.1.26