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

Diff of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

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

revision 1407 by gross, Mon Feb 4 06:45:48 2008 UTC revision 4346 by jfenwick, Tue Apr 2 04:46:45 2013 UTC
# Line 1  Line 1 
1    
2  /* $Id$ */  /*****************************************************************************
3    *
4  /*******************************************************  * Copyright (c) 2003-2013 by University of Queensland
5   *  * http://www.uq.edu.au
6   *           Copyright 2003-2007 by ACceSS MNRF  *
7   *       Copyright 2007 by University of Queensland  * Primary Business: Queensland, Australia
8   *  * Licensed under the Open Software License version 3.0
9   *                http://esscc.uq.edu.au  * http://www.opensource.org/licenses/osl-3.0.php
10   *        Primary Business: Queensland, Australia  *
11   *  Licensed under the Open Software License version 3.0  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12   *     http://www.opensource.org/licenses/osl-3.0.php  * Development since 2012 by School of Earth Sciences
13   *  *
14   *******************************************************/  *****************************************************************************/
15    
16    #include <pasowrap/PasoException.h>
17    #include <pasowrap/TransportProblemAdapter.h>
18  #include "MeshAdapter.h"  #include "MeshAdapter.h"
19  #include "escript/Data.h"  #include "escript/Data.h"
20  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
21  #ifdef USE_NETCDF  #ifdef USE_NETCDF
22  #include <netcdfcpp.h>  #include <netcdfcpp.h>
23  #endif  #endif
24  extern "C" {  #include "esysUtils/blocktimer.h"
25  #include "escript/blocktimer.h"  
26  }  #include <boost/python/import.hpp>
 #include <vector>  
27    
28  using namespace std;  using namespace std;
29  using namespace escript;  using namespace escript;
30    using namespace paso;
31    
32  namespace finley {  namespace finley {
33    
# 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);  
449    
450       // Points_Id        // Points_Id
451       if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
452          throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Id)");
453       int_ptr = &mesh->Points->Id[0];        int_ptr = &mesh->Points->Id[0];
454       if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
455          throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Id)");
456    
457       // Points_Tag        // Points_Tag
458       if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
459          throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Tag)");
460       int_ptr = &mesh->Points->Tag[0];        int_ptr = &mesh->Points->Tag[0];
461       if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
462          throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Tag)");
463    
464       // Points_Owner        // Points_Owner
465       if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
466          throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Owner)");
467       int_ptr = &mesh->Points->Owner[0];        int_ptr = &mesh->Points->Owner[0];
468       if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
469          throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Owner)");
470    
471       // Points_Color        // Points_Color
472       if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
473          throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Color)");
474       int_ptr = &mesh->Points->Color[0];        int_ptr = &mesh->Points->Color[0];
475       if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
476          throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Color)");
477    
478       // Points_Nodes        // Points_Nodes
479       // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
480       for (int i=0; i<num_Points; i++)        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
481         Points_Nodes[i] = mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]];           throw DataException(msgPrefix+"add_var(Points_Nodes)");
482       if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
483          throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Nodes)");
      if (! (ids->put(&(Points_Nodes[0]), num_Points)) )  
         throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);  
   
      TMPMEMFREE(Points_Nodes);  
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::getDiracDeltaFunctionsCode() 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,                                   AbstractSystemMatrix& 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,
777  {                                   const escript::Data& d_dirac,const escript::Data& y_dirac) const
778    {
779       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
780       if (smat==0)
781       {
782        throw FinleyAdapterException("finley only supports Paso system matrices.");
783       }
784     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
785     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
786     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 756  void MeshAdapter::addPDEToSystem( Line 792  void MeshAdapter::addPDEToSystem(
792     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
793     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
794     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
795       escriptDataC _d_dirac=d_dirac.getDataC();
796       escriptDataC _y_dirac=y_dirac.getDataC();
797    
798     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
799    
800     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
801       checkFinleyError();
802    
803       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
804     checkFinleyError();     checkFinleyError();
805    
806     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
807     checkFinleyError();     checkFinleyError();
808    
809     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );      Finley_Assemble_PDE(mesh->Nodes,mesh->Points, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
810     checkFinleyError();     checkFinleyError();
811  }  }
812    
813  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
814                       escript::Data& mat,                                          escript::Data& mat,
815                       const escript::Data& D,                                          const escript::Data& D,
816                       const escript::Data& d) const                                          const escript::Data& d,
817                                            const escript::Data& d_dirac,
818                                            const bool useHRZ) const
819  {  {
820     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
821     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
822     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
823       escriptDataC _d_dirac=d_dirac.getDataC();
824    
825     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
826    
827     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
828     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     checkFinleyError();
829      
830       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
831       checkFinleyError();
832    
833       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Points,&_mat, &_d_dirac, useHRZ);
834     checkFinleyError();     checkFinleyError();
835    
836  }  }
837    
838    
839  //  //
840  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
841  //  //
842  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact,  const escript::Data& y_dirac) const
843  {  {
844     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
845    
# Line 799  void MeshAdapter::addPDEToRHS( escript:: Line 848  void MeshAdapter::addPDEToRHS( escript::
848     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
849     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
850     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
851       escriptDataC _y_dirac=y_dirac.getDataC();
852    
853     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
854     checkFinleyError();     checkFinleyError();
# Line 808  void MeshAdapter::addPDEToRHS( escript:: Line 858  void MeshAdapter::addPDEToRHS( escript::
858    
859     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
860     checkFinleyError();     checkFinleyError();
861    
862       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs , 0, 0, 0, 0, 0, &_y_dirac );
863       checkFinleyError();
864    
865  }  }
866  //  //
867  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
868  //  //
869  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
870                       TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
871                       const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
872                       const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
873                       const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
874                       const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact,
875  {                                             const escript::Data& d_dirac, const escript::Data& y_dirac) const
876     DataArrayView::ShapeType shape;  {
877       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
878       if (tpa==0)
879       {
880        throw FinleyAdapterException("finley only supports Paso transport problems.");
881       }
882    
883    
884       DataTypes::ShapeType shape;
885     source.expand();     source.expand();
886     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
887     escriptDataC _M=M.getDataC();     escriptDataC _M=M.getDataC();
# Line 833  void MeshAdapter::addPDEToTransportProbl Line 895  void MeshAdapter::addPDEToTransportProbl
895     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
896     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
897     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
898       escriptDataC _d_dirac=d_dirac.getDataC();
899       escriptDataC _y_dirac=y_dirac.getDataC();
900    
901    
902     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
903     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
904      
905     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
906     checkFinleyError();     checkFinleyError();
907    
# Line 848  void MeshAdapter::addPDEToTransportProbl Line 913  void MeshAdapter::addPDEToTransportProbl
913    
914     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
915     checkFinleyError();     checkFinleyError();
916    
917       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
918       checkFinleyError();
919    
920  }  }
921    
922  //  //
# Line 855  void MeshAdapter::addPDEToTransportProbl Line 924  void MeshAdapter::addPDEToTransportProbl
924  //  //
925  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
926  {  {
927    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
928    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
929    if (inDomain!=*this)       if (inDomain!=*this)  
930      throw FinleyAdapterException("Error - Illegal domain of interpolant.");        throw FinleyAdapterException("Error - Illegal domain of interpolant.");
931    if (targetDomain!=*this)     if (targetDomain!=*this)
932      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");        throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
933    
934    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
935    escriptDataC _target=target.getDataC();     escriptDataC _target=target.getDataC();
936    escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
937    switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
938       case(Nodes):     case(Nodes):
939          switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
940             case(Nodes):        case(Nodes):
941             case(ReducedNodes):        case(ReducedNodes):
942             case(DegreesOfFreedom):        case(DegreesOfFreedom):
943             case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
944                 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
945                 break;        break;
946             case(Elements):        case(Elements):
947             case(ReducedElements):        case(ReducedElements):
948                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
949                 break;        break;
950             case(FaceElements):        case(FaceElements):
951             case(ReducedFaceElements):        case(ReducedFaceElements):
952                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
953                 break;        break;
954             case(Points):        case(Points):
955                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
956                 break;        break;
957             case(ContactElementsZero):        case(ContactElementsZero):
958             case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
959                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
960                 break;        break;
961             case(ContactElementsOne):        case(ContactElementsOne):
962             case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
963                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
964                 break;        break;
965             default:        default:
966                 stringstream temp;           stringstream temp;
967                 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();
968                 throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
969                 break;           break;
970          }        }
971          break;        break;
972       case(ReducedNodes):     case(ReducedNodes):
973          switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
974             case(Nodes):        case(Nodes):
975             case(ReducedNodes):        case(ReducedNodes):
976             case(DegreesOfFreedom):        case(DegreesOfFreedom):
977             case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
978                 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
979                 break;        break;
980             case(Elements):        case(Elements):
981             case(ReducedElements):        case(ReducedElements):
982                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
983                 break;        break;
984             case(FaceElements):        case(FaceElements):
985             case(ReducedFaceElements):        case(ReducedFaceElements):
986                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
987                 break;        break;
988             case(Points):        case(Points):
989                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
990                 break;        break;
991             case(ContactElementsZero):        case(ContactElementsZero):
992             case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
993                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
994                 break;        break;
995             case(ContactElementsOne):        case(ContactElementsOne):
996             case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
997                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
998                 break;        break;
999             default:        default:
1000                 stringstream temp;           stringstream temp;
1001                 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();
1002                 throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
1003                 break;           break;
1004          }        }
1005          break;        break;
1006       case(Elements):     case(Elements):
1007          if (target.getFunctionSpace().getTypeCode()==Elements) {        if (target.getFunctionSpace().getTypeCode()==Elements) {
1008             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1009          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1010             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
1011          } else {        } else {
1012             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
1013          }        }
1014          break;        break;
1015       case(ReducedElements):     case(ReducedElements):
1016          if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1017             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1018          } else {        } else {
1019             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.");
1020          }        }
1021          break;        break;
1022       case(FaceElements):     case(FaceElements):
1023          if (target.getFunctionSpace().getTypeCode()==FaceElements) {        if (target.getFunctionSpace().getTypeCode()==FaceElements) {
1024             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1025          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1026             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1027          } else {        } else {
1028             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1029          }        }
1030          break;        break;
1031       case(ReducedFaceElements):     case(ReducedFaceElements):
1032          if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1033             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1034          } else {        } else {
1035             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.");
1036         }        }
1037         break;        break;
1038       case(Points):     case(Points):
1039          if (target.getFunctionSpace().getTypeCode()==Points) {        if (target.getFunctionSpace().getTypeCode()==Points) {
1040             Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1041          } else {        } else {
1042             throw FinleyAdapterException("Error - No interpolation with data on points possible.");           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1043          }        }
1044          break;        break;
1045       case(ContactElementsZero):     case(ContactElementsZero):
1046       case(ContactElementsOne):     case(ContactElementsOne):
1047          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {        if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1048             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1049          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1050             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1051          } else {        } else {
1052             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1053          }        }
1054          break;        break;
1055       case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1056       case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1057          if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1058             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1059          } else {        } else {
1060             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.");
1061          }        }
1062          break;        break;
1063       case(DegreesOfFreedom):           case(DegreesOfFreedom):      
1064          switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
1065             case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
1066             case(DegreesOfFreedom):        case(DegreesOfFreedom):
1067                Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1068                break;        break;
1069      
1070             case(Nodes):        case(Nodes):
1071             case(ReducedNodes):        case(ReducedNodes):
1072                if (getMPISize()>1) {        if (getMPISize()>1) {
1073                    escript::Data temp=escript::Data(in);           escript::Data temp=escript::Data(in);
1074                    temp.expand();           temp.expand();
1075                    escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1076                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1077                } else {        } else {
1078                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1079                }        }
1080                break;        break;
1081             case(Elements):        case(Elements):
1082             case(ReducedElements):        case(ReducedElements):
1083                if (getMPISize()>1) {        if (getMPISize()>1) {
1084                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1085                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1086                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1087                } else {        } else {
1088                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1089                }        }
1090                break;        break;
1091             case(FaceElements):        case(FaceElements):
1092             case(ReducedFaceElements):        case(ReducedFaceElements):
1093                if (getMPISize()>1) {        if (getMPISize()>1) {
1094                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1095                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1096                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1097      
1098                } else {        } else {
1099                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1100                }        }
1101                break;        break;
1102             case(Points):        case(Points):
1103                if (getMPISize()>1) {        if (getMPISize()>1) {
1104                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1105                   escriptDataC _in2 = temp.getDataC();           //escriptDataC _in2 = temp.getDataC();
1106                } else {        } else {
1107                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1108                }        }
1109                break;        break;
1110             case(ContactElementsZero):        case(ContactElementsZero):
1111             case(ContactElementsOne):        case(ContactElementsOne):
1112             case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1113             case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1114                if (getMPISize()>1) {        if (getMPISize()>1) {
1115                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1116                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1117                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1118                } else {        } else {
1119                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1120                }        }
1121               break;        break;
1122             default:        default:
1123               stringstream temp;           stringstream temp;
1124               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();
1125               throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
1126               break;           break;
1127          }        }
1128          break;        break;
1129       case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1130         switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
1131            case(Nodes):        case(Nodes):
1132               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.");
1133               break;        break;
1134            case(ReducedNodes):        case(ReducedNodes):
1135                if (getMPISize()>1) {        if (getMPISize()>1) {
1136                    escript::Data temp=escript::Data(in);           escript::Data temp=escript::Data(in);
1137                    temp.expand();           temp.expand();
1138                    escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1139                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1140                } else {        } else {
1141                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1142                }        }
1143                break;        break;
1144            case(DegreesOfFreedom):        case(DegreesOfFreedom):
1145               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");
1146               break;        break;
1147            case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
1148               Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1149               break;        break;
1150            case(Elements):        case(Elements):
1151            case(ReducedElements):        case(ReducedElements):
1152                if (getMPISize()>1) {        if (getMPISize()>1) {
1153                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1154                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1155                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1156                } else {        } else {
1157                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1158               }        }
1159               break;        break;
1160            case(FaceElements):        case(FaceElements):
1161            case(ReducedFaceElements):        case(ReducedFaceElements):
1162                if (getMPISize()>1) {        if (getMPISize()>1) {
1163                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1164                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1165                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1166                } else {        } else {
1167                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1168                }        }
1169               break;        break;
1170            case(Points):        case(Points):
1171                if (getMPISize()>1) {        if (getMPISize()>1) {
1172                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1173                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1174                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1175                } else {        } else {
1176                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1177               }        }
1178               break;        break;
1179            case(ContactElementsZero):        case(ContactElementsZero):
1180            case(ContactElementsOne):        case(ContactElementsOne):
1181            case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1182            case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1183                if (getMPISize()>1) {        if (getMPISize()>1) {
1184                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1185                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1186                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1187                } else {        } else {
1188                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1189               }        }
1190               break;        break;
1191            default:        default:
1192               stringstream temp;           stringstream temp;
1193               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();
1194               throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
1195               break;           break;
1196         }        }
1197         break;        break;
1198       default:     default:
1199          stringstream temp;        stringstream temp;
1200          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();
1201          throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1202          break;        break;
1203    }     }
1204    checkFinleyError();     checkFinleyError();
1205  }  }
1206    
1207  //  //
# Line 1140  void MeshAdapter::interpolateOnDomain(es Line 1209  void MeshAdapter::interpolateOnDomain(es
1209  //  //
1210  void MeshAdapter::setToX(escript::Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
1211  {  {
1212    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1213    if (argDomain!=*this)     if (argDomain!=*this)
1214       throw FinleyAdapterException("Error - Illegal domain of data point locations");        throw FinleyAdapterException("Error - Illegal domain of data point locations");
1215    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1216    // in case of values node coordinates we can do the job directly:     // in case of values node coordinates we can do the job directly:
1217    if (arg.getFunctionSpace().getTypeCode()==Nodes) {     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1218       escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1219       Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1220    } else {     } else {
1221       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1222       escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1223       Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1224       // this is then interpolated onto arg:        // this is then interpolated onto arg:
1225       interpolateOnDomain(arg,tmp_data);        interpolateOnDomain(arg,tmp_data);
1226    }     }
1227    checkFinleyError();     checkFinleyError();
1228  }  }
1229    
1230  //  //
# Line 1163  void MeshAdapter::setToX(escript::Data& Line 1232  void MeshAdapter::setToX(escript::Data&
1232  //  //
1233  void MeshAdapter::setToNormal(escript::Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
1234  {  {
1235    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1236    if (normalDomain!=*this)     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1237       throw FinleyAdapterException("Error - Illegal domain of normal locations");     if (normalDomain!=*this)
1238    Finley_Mesh* mesh=m_finleyMesh.get();        throw FinleyAdapterException("Error - Illegal domain of normal locations");
1239    escriptDataC _normal=normal.getDataC();     Finley_Mesh* mesh=m_finleyMesh.get();
1240    switch(normal.getFunctionSpace().getTypeCode()) {     escriptDataC _normal=normal.getDataC();
1241      case(Nodes):     switch(normal.getFunctionSpace().getTypeCode()) {
1242        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");     case(Nodes):
1243        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1244      case(ReducedNodes):     break;
1245        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");     case(ReducedNodes):
1246        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1247      case(Elements):     break;
1248        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");     case(Elements):
1249        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1250      case(ReducedElements):     break;
1251        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");     case(ReducedElements):
1252        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1253      case (FaceElements):     break;
1254        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     case (FaceElements):
1255        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1256      case (ReducedFaceElements):     break;
1257        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     case (ReducedFaceElements):
1258        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1259      case(Points):     break;
1260        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");     case(Points):
1261        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1262      case (ContactElementsOne):     break;
1263      case (ContactElementsZero):     case (ContactElementsOne):
1264        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);     case (ContactElementsZero):
1265        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1266      case (ReducedContactElementsOne):     break;
1267      case (ReducedContactElementsZero):     case (ReducedContactElementsOne):
1268        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);     case (ReducedContactElementsZero):
1269        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1270      case(DegreesOfFreedom):     break;
1271        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");     case(DegreesOfFreedom):
1272        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1273      case(ReducedDegreesOfFreedom):     break;
1274        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");     case(ReducedDegreesOfFreedom):
1275        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1276      default:     break;
1277       default:
1278        stringstream temp;        stringstream temp;
1279        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();
1280        throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1281        break;        break;
1282    }     }
1283    checkFinleyError();     checkFinleyError();
1284  }  }
1285    
1286  //  //
# Line 1218  void MeshAdapter::setToNormal(escript::D Line 1288  void MeshAdapter::setToNormal(escript::D
1288  //  //
1289  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1290  {  {
1291    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1292    if (targetDomain!=*this)     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1293       throw FinleyAdapterException("Error - Illegal domain of interpolation target");     if (targetDomain!=this)
1294          throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1295    
1296    throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
1297  }  }
1298    
1299  //  //
1300  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
1301  //  //
1302  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1303  {  {
1304    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1305    if (argDomain!=*this)     if (argDomain!=*this)
1306       throw FinleyAdapterException("Error - Illegal domain of integration kernel");        throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1307    
1308    double blocktimer_start = blocktimer_time();     double blocktimer_start = blocktimer_time();
1309    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1310    escriptDataC _temp;     escriptDataC _temp;
1311    escript::Data temp;     escript::Data temp;
1312    escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1313    switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1314       case(Nodes):     case(Nodes):
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(ReducedNodes):     case(ReducedNodes):
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       case(Elements):     case(Elements):
1325          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1326          break;     break;
1327       case(ReducedElements):     case(ReducedElements):
1328          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1329          break;     break;
1330       case(FaceElements):     case(FaceElements):
1331          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1332          break;     break;
1333       case(ReducedFaceElements):     case(ReducedFaceElements):
1334          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1335          break;     break;
1336       case(Points):     case(Points):
1337          throw FinleyAdapterException("Error - Integral of data on points is not supported.");     throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1338          break;     break;
1339       case(ContactElementsZero):     case(ContactElementsZero):
1340          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1341          break;     break;
1342       case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1343          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1344          break;     break;
1345       case(ContactElementsOne):     case(ContactElementsOne):
1346          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1347          break;     break;
1348       case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1349          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1350          break;     break;
1351       case(DegreesOfFreedom):     case(DegreesOfFreedom):
1352          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1353          _temp=temp.getDataC();     _temp=temp.getDataC();
1354          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1355          break;     break;
1356       case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1357          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1358          _temp=temp.getDataC();     _temp=temp.getDataC();
1359          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1360          break;     break;
1361       default:     default:
1362          stringstream temp;        stringstream temp;
1363          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();
1364          throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1365          break;        break;
1366    }     }
1367    checkFinleyError();     checkFinleyError();
1368    blocktimer_increment("integrate()", blocktimer_start);     blocktimer_increment("integrate()", blocktimer_start);
1369  }  }
1370    
1371  //  //
# Line 1302  void MeshAdapter::setToIntegrals(std::ve Line 1373  void MeshAdapter::setToIntegrals(std::ve
1373  //  //
1374  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1375  {  {
1376    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1377    if (argDomain!=*this)     if (argDomain!=*this)
1378      throw FinleyAdapterException("Error - Illegal domain of gradient argument");        throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1379    const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1380    if (gradDomain!=*this)     if (gradDomain!=*this)
1381       throw FinleyAdapterException("Error - Illegal domain of gradient");        throw FinleyAdapterException("Error - Illegal domain of gradient");
1382    
1383    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1384    escriptDataC _grad=grad.getDataC();     escriptDataC _grad=grad.getDataC();
1385    escriptDataC nodeDataC;     escriptDataC nodeDataC;
1386    escript::Data temp;     escript::Data temp;
1387    if (getMPISize()>1) {     if (getMPISize()>1) {
1388        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1389          temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1390          nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1391        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1392          temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1393          nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1394        } else {        } else {
1395          nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
1396        }        }
1397    } else {     } else {
1398       nodeDataC = arg.getDataC();        nodeDataC = arg.getDataC();
1399    }     }
1400    switch(grad.getFunctionSpace().getTypeCode()) {     switch(grad.getFunctionSpace().getTypeCode()) {
1401         case(Nodes):     case(Nodes):
1402            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");     throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1403            break;     break;
1404         case(ReducedNodes):     case(ReducedNodes):
1405            throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1406            break;     break;
1407         case(Elements):     case(Elements):
1408            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1409            break;     break;
1410         case(ReducedElements):     case(ReducedElements):
1411            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1412            break;     break;
1413         case(FaceElements):     case(FaceElements):
1414            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1415            break;     break;
1416         case(ReducedFaceElements):     case(ReducedFaceElements):
1417            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1418            break;     break;
1419         case(Points):     case(Points):
1420            throw FinleyAdapterException("Error - Gradient at points is not supported.");     throw FinleyAdapterException("Error - Gradient at points is not supported.");
1421            break;     break;
1422         case(ContactElementsZero):     case(ContactElementsZero):
1423            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1424            break;     break;
1425         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1426            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1427            break;     break;
1428         case(ContactElementsOne):     case(ContactElementsOne):
1429            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1430            break;     break;
1431         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1432            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1433            break;     break;
1434         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1435            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1436            break;     break;
1437         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1438            throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1439            break;     break;
1440         default:     default:
1441            stringstream temp;        stringstream temp;
1442            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();
1443            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1444            break;        break;
1445    }     }
1446    checkFinleyError();     checkFinleyError();
1447  }  }
1448    
1449  //  //
# Line 1380  void MeshAdapter::setToGradient(escript: Line 1451  void MeshAdapter::setToGradient(escript:
1451  //  //
1452  void MeshAdapter::setToSize(escript::Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
1453  {  {
1454    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1455    escriptDataC tmp=size.getDataC();     escriptDataC tmp=size.getDataC();
1456    switch(size.getFunctionSpace().getTypeCode()) {     switch(size.getFunctionSpace().getTypeCode()) {
1457         case(Nodes):     case(Nodes):
1458            throw FinleyAdapterException("Error - Size of nodes is not supported.");     throw FinleyAdapterException("Error - Size of nodes is not supported.");
1459            break;     break;
1460         case(ReducedNodes):     case(ReducedNodes):
1461            throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1462            break;     break;
1463         case(Elements):     case(Elements):
1464            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1465            break;     break;
1466         case(ReducedElements):     case(ReducedElements):
1467            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1468            break;     break;
1469         case(FaceElements):     case(FaceElements):
1470            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1471            break;     break;
1472         case(ReducedFaceElements):     case(ReducedFaceElements):
1473            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1474            break;     break;
1475         case(Points):     case(Points):
1476            throw FinleyAdapterException("Error - Size of point elements is not supported.");     throw FinleyAdapterException("Error - Size of point elements is not supported.");
1477            break;     break;
1478         case(ContactElementsZero):     case(ContactElementsZero):
1479         case(ContactElementsOne):     case(ContactElementsOne):
1480            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1481            break;     break;
1482         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1483         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1484            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1485            break;     break;
1486         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1487            throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1488            break;     break;
1489         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1490            throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1491            break;     break;
1492         default:     default:
1493            stringstream temp;        stringstream temp;
1494            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();
1495            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1496            break;        break;
1497    }     }
1498    checkFinleyError();     checkFinleyError();
1499  }  }
1500    
1501  // sets the location of nodes:  //
1502    // sets the location of nodes
1503    //
1504  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1505  {  {
1506    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1507    escriptDataC tmp;     escriptDataC tmp;
1508    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1509    if (newDomain!=*this)     if (newDomain!=*this)
1510       throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1511    tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1512    Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1513    checkFinleyError();         Finley_Mesh_setCoordinates(mesh,&tmp);
1514  }     } else {
1515           throw FinleyAdapterException("As of version escript3.3 - SetNewX only accepts ContinuousFunction arguments please interpolate.");    
1516  // saves a data array in openDX format:  /*       escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1517  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const         tmp = new_x_inter.getDataC();
1518  {         Finley_Mesh_setCoordinates(mesh,&tmp);*/
1519      int MAX_namelength=256;     }
1520      const int num_data=boost::python::extract<int>(arg.attr("__len__")());     checkFinleyError();
   /* 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);  
   
     return;  
1521  }  }
1522    
1523  // saves a data array in openVTK format:  bool MeshAdapter::ownSample(int fs_code, index_t id) const
 void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  
1524  {  {
1525      int MAX_namelength=256;      if (getMPISize() > 1) {
1526      const int num_data=boost::python::extract<int>(arg.attr("__len__")());  #ifdef ESYS_MPI
1527    /* win32 refactor */          index_t myFirstNode=0, myLastNode=0, k=0;
1528    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;          index_t* globalNodeIndex=0;
1529    for(int i=0;i<num_data;i++)          Finley_Mesh* mesh_p=m_finleyMesh.get();
1530    {          /*
1531      names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;           * this method is only used by saveDataCSV which would use the returned
1532    }           * values for reduced nodes wrongly so this case is disabled for now
1533            if (fs_code == FINLEY_REDUCED_NODES)
1534            {
1535                myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1536                myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1537                globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1538            }
1539            else
1540            */
1541            if (fs_code == FINLEY_NODES)
1542            {
1543                myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1544                myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1545                globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1546            }
1547            else
1548            {
1549                throw FinleyAdapterException("Unsupported function space type for ownSample()");
1550            }
1551    
1552    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;          k=globalNodeIndex[id];
1553    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;          return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1554    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  #endif
   
     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());  
          }  
1555      }      }
1556      Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);      return true;
1557    }
1558    
   checkFinleyError();  
   /* win32 refactor */  
   TMPMEMFREE(c_names);  
   TMPMEMFREE(data);  
   TMPMEMFREE(ptr_data);  
   for(int i=0;i<num_data;i++)  
   {  
     TMPMEMFREE(names[i]);  
   }  
   TMPMEMFREE(names);  
1559    
1560      return;  //
1561  }  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1562                                                                                                                                                                            //
1563                                                                                                                                                                            ASM_ptr MeshAdapter::newSystemMatrix(
1564  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:                                                   const int row_blocksize,
1565  SystemMatrixAdapter MeshAdapter::newSystemMatrix(                                                   const escript::FunctionSpace& row_functionspace,
1566                        const int row_blocksize,                                                   const int column_blocksize,
1567                        const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& column_functionspace,
1568                        const int column_blocksize,                                                   const int type) const
1569                        const escript::FunctionSpace& column_functionspace,  {
1570                        const int type) const     int reduceRowOrder=0;
1571  {     int reduceColOrder=0;
1572      int reduceRowOrder=0;     // is the domain right?
1573      int reduceColOrder=0;     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1574      // is the domain right?     if (row_domain!=*this)
1575      const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1576      if (row_domain!=*this)     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1577            throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");     if (col_domain!=*this)
1578      const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());        throw FinleyAdapterException("Error - domain of column function space does not match the domain of matrix generator.");
1579      if (col_domain!=*this)     // is the function space type right
1580            throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1581      // is the function space type right        reduceRowOrder=0;
1582      if (row_functionspace.getTypeCode()==DegreesOfFreedom) {     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1583          reduceRowOrder=0;        reduceRowOrder=1;
1584      } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     } else {
1585          reduceRowOrder=1;        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1586      } else {     }
1587          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1588      }        reduceColOrder=0;
1589      if (column_functionspace.getTypeCode()==DegreesOfFreedom) {     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1590          reduceColOrder=0;        reduceColOrder=1;
1591      } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     } else {
1592          reduceColOrder=1;        throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1593      } else {     }
1594          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");     // generate matrix:
1595      }  
1596      // generate matrix:     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1597           checkFinleyError();
1598      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);     Paso_SystemMatrix* fsystemMatrix;
1599      checkFinleyError();     int trilinos = 0;
1600      Paso_SystemMatrix* fsystemMatrix;     if (trilinos) {
     int trilinos = 0;  
     if (trilinos) {  
1601  #ifdef TRILINOS  #ifdef TRILINOS
1602        /* Allocation Epetra_VrbMatrix here */        /* Allocation Epetra_VrbMatrix here */
1603  #endif  #endif
1604      }     }
1605      else {     else {
1606        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1607      }     }
1608      checkPasoError();     checkPasoError();
1609      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1610      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1611       return ASM_ptr(sma);
1612    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1613  }  }
1614    
1615    //
1616  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1617  TransportProblemAdapter MeshAdapter::newTransportProblem(  //
1618                        const double theta,  ATP_ptr MeshAdapter::newTransportProblem(
1619                        const int blocksize,                                                           const int blocksize,
1620                        const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1621                        const int type) const                                                           const int type) const
1622  {  {
1623      int reduceOrder=0;     int reduceOrder=0;
1624      // is the domain right?     // is the domain right?
1625      const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(functionspace.getDomain());     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1626      if (domain!=*this)     if (domain!=*this)
1627            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.");
1628      // is the function space type right     // is the function space type right
1629      if (functionspace.getTypeCode()==DegreesOfFreedom) {     if (functionspace.getTypeCode()==DegreesOfFreedom) {
1630          reduceOrder=0;        reduceOrder=0;
1631      } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1632          reduceOrder=1;        reduceOrder=1;
1633      } else {     } else {
1634          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1635      }     }
1636      // generate matrix:     // generate matrix:
1637        
1638      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1639      checkFinleyError();     checkFinleyError();
1640      Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1641      transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1642      checkPasoError();     checkPasoError();
1643      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1644      return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1645       return ATP_ptr(tpa);
1646    //   return TransportProblemAdapter(transportProblem, blocksize,functionspace);
1647  }  }
1648    
1649  //  //
# Line 1627  TransportProblemAdapter MeshAdapter::new Line 1655  TransportProblemAdapter MeshAdapter::new
1655  // 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:
1656  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1657  {  {
1658    switch(functionSpaceCode) {     switch(functionSpaceCode) {
1659         case(Nodes):     case(Nodes):
1660         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1661         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1662            return false;     return false;
1663            break;     break;
1664         case(Elements):     case(Elements):
1665         case(FaceElements):     case(FaceElements):
1666         case(Points):     case(Points):
1667         case(ContactElementsZero):     case(ContactElementsZero):
1668         case(ContactElementsOne):     case(ContactElementsOne):
1669         case(ReducedElements):     case(ReducedElements):
1670         case(ReducedFaceElements):     case(ReducedFaceElements):
1671         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1672         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1673            return true;     return true;
1674            break;     break;
1675         default:     default:
1676            stringstream temp;        stringstream temp;
1677            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;
1678            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1679            break;        break;
1680    }     }
1681    checkFinleyError();     checkFinleyError();
1682    return false;     return false;
1683    }
1684    
1685    bool
1686    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1687    {
1688       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1689        class 1: DOF <-> Nodes
1690        class 2: ReducedDOF <-> ReducedNodes
1691        class 3: Points
1692        class 4: Elements
1693        class 5: ReducedElements
1694        class 6: FaceElements
1695        class 7: ReducedFaceElements
1696        class 8: ContactElementZero <-> ContactElementOne
1697        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1698    
1699       There is also a set of lines. Interpolation is possible down a line but not between lines.
1700       class 1 and 2 belong to all lines so aren't considered.
1701        line 0: class 3
1702        line 1: class 4,5
1703        line 2: class 6,7
1704        line 3: class 8,9
1705    
1706       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1707       eg hasnodes is true if we have at least one instance of Nodes.
1708       */
1709        if (fs.empty())
1710        {
1711            return false;
1712        }
1713        vector<int> hasclass(10);
1714        vector<int> hasline(4);
1715        bool hasnodes=false;
1716        bool hasrednodes=false;
1717        bool hascez=false;
1718        bool hasrcez=false;
1719        for (int i=0;i<fs.size();++i)
1720        {
1721        switch(fs[i])
1722        {
1723        case(Nodes):   hasnodes=true;   // no break is deliberate
1724        case(DegreesOfFreedom):
1725            hasclass[1]=1;
1726            break;
1727        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1728        case(ReducedDegreesOfFreedom):
1729            hasclass[2]=1;
1730            break;
1731        case(Points):
1732            hasline[0]=1;
1733            hasclass[3]=1;
1734            break;
1735        case(Elements):
1736            hasclass[4]=1;
1737            hasline[1]=1;
1738            break;
1739        case(ReducedElements):
1740            hasclass[5]=1;
1741            hasline[1]=1;
1742            break;
1743        case(FaceElements):
1744            hasclass[6]=1;
1745            hasline[2]=1;
1746            break;
1747        case(ReducedFaceElements):
1748            hasclass[7]=1;
1749            hasline[2]=1;
1750            break;
1751        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1752        case(ContactElementsOne):
1753            hasclass[8]=1;
1754            hasline[3]=1;
1755            break;
1756        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1757        case(ReducedContactElementsOne):
1758            hasclass[9]=1;
1759            hasline[3]=1;
1760            break;
1761        default:
1762            return false;
1763        }
1764        }
1765        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1766        // fail if we have more than one leaf group
1767    
1768        if (totlines>1)
1769        {
1770        return false;   // there are at least two branches we can't interpolate between
1771        }
1772        else if (totlines==1)
1773        {
1774        if (hasline[0]==1)      // we have points
1775        {
1776            resultcode=Points;
1777        }
1778        else if (hasline[1]==1)
1779        {
1780            if (hasclass[5]==1)
1781            {
1782            resultcode=ReducedElements;
1783            }
1784            else
1785            {
1786            resultcode=Elements;
1787            }
1788        }
1789        else if (hasline[2]==1)
1790        {
1791            if (hasclass[7]==1)
1792            {
1793            resultcode=ReducedFaceElements;
1794            }
1795            else
1796            {
1797            resultcode=FaceElements;
1798            }
1799        }
1800        else    // so we must be in line3
1801        {
1802            if (hasclass[9]==1)
1803            {
1804            // need something from class 9
1805            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1806            }
1807            else
1808            {
1809            // something from class 8
1810            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1811            }
1812        }
1813        }
1814        else    // totlines==0
1815        {
1816        if (hasclass[2]==1)
1817        {
1818            // something from class 2
1819            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1820        }
1821        else
1822        {   // something from class 1
1823            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1824        }
1825        }
1826        return true;
1827  }  }
1828    
1829  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1830  {  {
1831    switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1832       case(Nodes):     case(Nodes):
1833          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1834             case(Nodes):      case(Nodes):
1835             case(ReducedNodes):      case(ReducedNodes):
1836             case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1837             case(DegreesOfFreedom):      case(DegreesOfFreedom):
1838             case(Elements):      case(Elements):
1839             case(ReducedElements):      case(ReducedElements):
1840             case(FaceElements):      case(FaceElements):
1841             case(ReducedFaceElements):      case(ReducedFaceElements):
1842             case(Points):      case(Points):
1843             case(ContactElementsZero):      case(ContactElementsZero):
1844             case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1845             case(ContactElementsOne):      case(ContactElementsOne):
1846             case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1847                 return true;      return true;
1848             default:      default:
1849                 stringstream temp;            stringstream temp;
1850                 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;
1851                 throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1852          }     }
1853          break;     break;
1854       case(ReducedNodes):     case(ReducedNodes):
1855          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1856             case(ReducedNodes):      case(ReducedNodes):
1857             case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1858             case(Elements):      case(Elements):
1859             case(ReducedElements):      case(ReducedElements):
1860             case(FaceElements):      case(FaceElements):
1861             case(ReducedFaceElements):      case(ReducedFaceElements):
1862             case(Points):      case(Points):
1863             case(ContactElementsZero):      case(ContactElementsZero):
1864             case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1865             case(ContactElementsOne):      case(ContactElementsOne):
1866             case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1867                 return true;      return true;
1868            case(Nodes):      case(Nodes):
1869            case(DegreesOfFreedom):      case(DegreesOfFreedom):
1870               return false;      return false;
1871             default:      default:
1872                 stringstream temp;          stringstream temp;
1873                 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;
1874                 throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1875          }     }
1876          break;     break;
1877       case(Elements):     case(Elements):
1878          if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1879             return true;        return true;
1880          } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1881             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;  
1882          } else {          } else {
1883             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());  
1884          }          }
1885          break;     case(ReducedElements):
1886       case(ReducedDegreesOfFreedom):      if (functionSpaceType_target==ReducedElements) {
1887         switch(functionSpaceType_target) {        return true;
1888            case(ReducedDegreesOfFreedom):      } else {
1889            case(ReducedNodes):            return false;
1890            case(Elements):      }
1891            case(ReducedElements):     case(FaceElements):
1892            case(FaceElements):      if (functionSpaceType_target==FaceElements) {
1893            case(ReducedFaceElements):              return true;
1894            case(Points):      } else if (functionSpaceType_target==ReducedFaceElements) {
1895            case(ContactElementsZero):              return true;
1896            case(ReducedContactElementsZero):      } else {
1897            case(ContactElementsOne):              return false;
1898            case(ReducedContactElementsOne):      }
1899                return true;     case(ReducedFaceElements):
1900            case(Nodes):      if (functionSpaceType_target==ReducedFaceElements) {
1901            case(DegreesOfFreedom):              return true;
1902               return false;      } else {
1903            default:          return false;
1904               stringstream temp;      }
1905               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;     case(Points):
1906               throw FinleyAdapterException(temp.str());      if (functionSpaceType_target==Points) {
1907         }              return true;
1908         break;      } else {
1909       default:              return false;
1910          stringstream temp;      }
1911          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;     case(ContactElementsZero):
1912          throw FinleyAdapterException(temp.str());     case(ContactElementsOne):
1913          break;      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1914    }              return true;
1915    checkFinleyError();      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1916    return false;              return true;
1917        } else {
1918                return false;
1919        }
1920       case(ReducedContactElementsZero):
1921       case(ReducedContactElementsOne):
1922        if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1923                return true;
1924        } else {
1925                return false;
1926        }
1927       case(DegreesOfFreedom):
1928        switch(functionSpaceType_target) {
1929        case(ReducedDegreesOfFreedom):
1930        case(DegreesOfFreedom):
1931        case(Nodes):
1932        case(ReducedNodes):
1933        case(Elements):
1934        case(ReducedElements):
1935        case(Points):
1936        case(FaceElements):
1937        case(ReducedFaceElements):
1938        case(ContactElementsZero):
1939        case(ReducedContactElementsZero):
1940        case(ContactElementsOne):
1941        case(ReducedContactElementsOne):
1942        return true;
1943        default:
1944            stringstream temp;
1945            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1946            throw FinleyAdapterException(temp.str());
1947        }
1948        break;
1949       case(ReducedDegreesOfFreedom):
1950       switch(functionSpaceType_target) {
1951        case(ReducedDegreesOfFreedom):
1952        case(ReducedNodes):
1953        case(Elements):
1954        case(ReducedElements):
1955        case(FaceElements):
1956        case(ReducedFaceElements):
1957        case(Points):
1958        case(ContactElementsZero):
1959        case(ReducedContactElementsZero):
1960        case(ContactElementsOne):
1961        case(ReducedContactElementsOne):
1962        return true;
1963        case(Nodes):
1964        case(DegreesOfFreedom):
1965        return false;
1966        default:
1967            stringstream temp;
1968            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1969            throw FinleyAdapterException(temp.str());
1970        }
1971        break;
1972       default:
1973          stringstream temp;
1974          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1975          throw FinleyAdapterException(temp.str());
1976          break;
1977       }
1978       checkFinleyError();
1979       return false;
1980  }  }
1981    
1982    signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1983    {
1984        if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1985        {
1986            return 1;
1987        }
1988        else if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1989        {
1990            return -1;
1991        }
1992        return 0;
1993    }  
1994      
1995  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
1996  {  {
1997      return false;     return false;
1998  }  }
1999    
2000  bool MeshAdapter::operator==(const AbstractDomain& other) const  bool MeshAdapter::operator==(const AbstractDomain& other) const
2001  {  {
2002    const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);     const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2003    if (temp!=0) {     if (temp!=0) {
2004      return (m_finleyMesh==temp->m_finleyMesh);        return (m_finleyMesh==temp->m_finleyMesh);
2005    } else {     } else {
2006      return false;        return false;
2007    }     }
2008  }  }
2009    
2010  bool MeshAdapter::operator!=(const AbstractDomain& other) const  bool MeshAdapter::operator!=(const AbstractDomain& other) const
2011  {  {
2012    return !(operator==(other));     return !(operator==(other));
2013  }  }
2014    
2015  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
2016  {  {
2017     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2018     checkPasoError();     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2019     return out;             package, symmetry, mesh->MPIInfo);
2020    }
2021    
2022    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2023    {
2024       Finley_Mesh* mesh=m_finleyMesh.get();
2025       return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2026               package, symmetry, mesh->MPIInfo);
2027  }  }
2028    
2029  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2030  {  {
2031    return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2032  }  }
2033    
2034  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2035  {  {
2036    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2037  }  }
2038    
2039  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2040  {  {
2041    return function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2042  }  }
2043    
2044  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2045  {  {
2046    int *out=0,i;     int *out = NULL;
2047    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2048    switch (functionSpaceType) {     switch (functionSpaceType) {
2049      case(Nodes):     case(Nodes):
2050        out=mesh->Nodes->Id;     out=mesh->Nodes->Id;
2051        break;     break;
2052      case(ReducedNodes):     case(ReducedNodes):
2053        out=mesh->Nodes->reducedNodesId;     out=mesh->Nodes->reducedNodesId;
2054        break;     break;
2055      case(Elements):     case(Elements):
2056        out=mesh->Elements->Id;     out=mesh->Elements->Id;
2057        break;     break;
2058      case(ReducedElements):     case(ReducedElements):
2059        out=mesh->Elements->Id;     out=mesh->Elements->Id;
2060        break;     break;
2061      case(FaceElements):     case(FaceElements):
2062        out=mesh->FaceElements->Id;     out=mesh->FaceElements->Id;
2063        break;     break;
2064      case(ReducedFaceElements):     case(ReducedFaceElements):
2065        out=mesh->FaceElements->Id;     out=mesh->FaceElements->Id;
2066        break;     break;
2067      case(Points):     case(Points):
2068        out=mesh->Points->Id;     out=mesh->Points->Id;
2069        break;     break;
2070      case(ContactElementsZero):     case(ContactElementsZero):
2071        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2072        break;     break;
2073      case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
2074        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2075        break;     break;
2076      case(ContactElementsOne):     case(ContactElementsOne):
2077        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2078        break;     break;
2079      case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
2080        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2081        break;     break;
2082      case(DegreesOfFreedom):     case(DegreesOfFreedom):
2083        out=mesh->Nodes->degreesOfFreedomId;     out=mesh->Nodes->degreesOfFreedomId;
2084        break;     break;
2085      case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2086        out=mesh->Nodes->reducedDegreesOfFreedomId;     out=mesh->Nodes->reducedDegreesOfFreedomId;
2087        break;     break;
2088      default:     default:
2089        stringstream temp;        stringstream temp;
2090        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2091        throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2092        break;        break;
2093    }     }
2094    return out;     return out;
2095  }  }
2096  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2097  {  {
2098    int out=0;     int out=0;
2099    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2100    switch (functionSpaceType) {     switch (functionSpaceType) {
2101    case(Nodes):     case(Nodes):
2102      out=mesh->Nodes->Tag[sampleNo];     out=mesh->Nodes->Tag[sampleNo];
2103      break;     break;
2104    case(ReducedNodes):     case(ReducedNodes):
2105      throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2106      break;     break;
2107    case(Elements):     case(Elements):
2108      out=mesh->Elements->Tag[sampleNo];     out=mesh->Elements->Tag[sampleNo];
2109      break;     break;
2110    case(ReducedElements):     case(ReducedElements):
2111      out=mesh->Elements->Tag[sampleNo];     out=mesh->Elements->Tag[sampleNo];
2112      break;     break;
2113    case(FaceElements):     case(FaceElements):
2114      out=mesh->FaceElements->Tag[sampleNo];     out=mesh->FaceElements->Tag[sampleNo];
2115      break;     break;
2116    case(ReducedFaceElements):     case(ReducedFaceElements):
2117      out=mesh->FaceElements->Tag[sampleNo];     out=mesh->FaceElements->Tag[sampleNo];
2118      break;     break;
2119    case(Points):     case(Points):
2120      out=mesh->Points->Tag[sampleNo];     out=mesh->Points->Tag[sampleNo];
2121      break;     break;
2122    case(ContactElementsZero):     case(ContactElementsZero):
2123      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2124      break;     break;
2125    case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
2126      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2127      break;     break;
2128    case(ContactElementsOne):     case(ContactElementsOne):
2129      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2130      break;     break;
2131    case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
2132      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2133      break;     break;
2134    case(DegreesOfFreedom):     case(DegreesOfFreedom):
2135      throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2136      break;     break;
2137    case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2138      throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");     throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2139      break;     break;
2140    default:     default:
2141      stringstream temp;        stringstream temp;
2142      temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2143      throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2144      break;        break;
2145    }     }
2146    return out;     return out;
2147  }  }
2148    
2149    
2150  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
2151  {  {
2152       Finley_Mesh* mesh=m_finleyMesh.get();
2153       escriptDataC tmp=mask.getDataC();
2154       switch(functionSpaceType) {
2155       case(Nodes):
2156       Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2157       break;
2158       case(ReducedNodes):
2159       throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2160       break;
2161       case(DegreesOfFreedom):
2162       throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2163       break;
2164       case(ReducedDegreesOfFreedom):
2165       throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2166       break;
2167       case(Elements):
2168       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2169       break;
2170       case(ReducedElements):
2171       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2172       break;
2173       case(FaceElements):
2174       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2175       break;
2176       case(ReducedFaceElements):
2177       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2178       break;
2179       case(Points):
2180       Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2181       break;
2182       case(ContactElementsZero):
2183       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2184       break;
2185       case(ReducedContactElementsZero):
2186       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2187       break;
2188       case(ContactElementsOne):
2189       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2190       break;
2191       case(ReducedContactElementsOne):
2192       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2193       break;
2194       default:
2195          stringstream temp;
2196          temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2197          throw FinleyAdapterException(temp.str());
2198       }
2199       checkFinleyError();
2200       return;
2201    }
2202    
2203    void MeshAdapter::setTagMap(const string& name,  int tag)
2204    {
2205       Finley_Mesh* mesh=m_finleyMesh.get();
2206       Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2207       checkFinleyError();
2208       // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2209    }
2210    
2211    int MeshAdapter::getTag(const string& name) const
2212    {
2213       Finley_Mesh* mesh=m_finleyMesh.get();
2214       int tag=0;
2215       tag=Finley_Mesh_getTag(mesh, name.c_str());
2216       checkFinleyError();
2217       // throwStandardException("MeshAdapter::getTag is not implemented.");
2218       return tag;
2219    }
2220    
2221    bool MeshAdapter::isValidTagName(const string& name) const
2222    {
2223       Finley_Mesh* mesh=m_finleyMesh.get();
2224       return Finley_Mesh_isValidTagName(mesh,name.c_str());
2225    }
2226    
2227    string MeshAdapter::showTagNames() const
2228    {
2229       stringstream temp;
2230       Finley_Mesh* mesh=m_finleyMesh.get();
2231       Finley_TagMap* tag_map=mesh->TagMap;
2232       while (tag_map) {
2233          temp << tag_map->name;
2234          tag_map=tag_map->next;
2235          if (tag_map) temp << ", ";
2236       }
2237       return temp.str();
2238    }
2239    
2240    int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2241    {
2242    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2243    escriptDataC tmp=mask.getDataC();    dim_t numTags=0;
2244    switch(functionSpaceType) {    switch(functionSpaceCode) {
2245         case(Nodes):     case(Nodes):
2246            Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);            numTags=mesh->Nodes->numTagsInUse;
2247            break;            break;
2248         case(ReducedNodes):     case(ReducedNodes):
2249            throw FinleyAdapterException("Error - ReducedNodes does not support tags");            throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2250            break;            break;
2251         case(DegreesOfFreedom):     case(DegreesOfFreedom):
2252            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2253            break;            break;
2254         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2255            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2256            break;            break;
2257         case(Elements):     case(Elements):
2258            Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     case(ReducedElements):
2259            break;            numTags=mesh->Elements->numTagsInUse;
2260         case(ReducedElements):            break;
2261            Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     case(FaceElements):
2262            break;     case(ReducedFaceElements):
2263         case(FaceElements):            numTags=mesh->FaceElements->numTagsInUse;
2264            Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);            break;
2265       case(Points):
2266              numTags=mesh->Points->numTagsInUse;
2267              break;
2268       case(ContactElementsZero):
2269       case(ReducedContactElementsZero):
2270       case(ContactElementsOne):
2271       case(ReducedContactElementsOne):
2272              numTags=mesh->ContactElements->numTagsInUse;
2273            break;            break;
2274         case(ReducedFaceElements):     default:
2275            Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);        stringstream temp;
2276            break;        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2277         case(Points):        throw FinleyAdapterException(temp.str());
2278            Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);    }
2279      return numTags;
2280    }
2281    
2282    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2283    {
2284      Finley_Mesh* mesh=m_finleyMesh.get();
2285      index_t* tags=NULL;
2286      switch(functionSpaceCode) {
2287       case(Nodes):
2288              tags=mesh->Nodes->tagsInUse;
2289            break;            break;
2290         case(ContactElementsZero):     case(ReducedNodes):
2291            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2292            break;            break;
2293         case(ReducedContactElementsZero):     case(DegreesOfFreedom):
2294            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2295            break;            break;
2296         case(ContactElementsOne):     case(ReducedDegreesOfFreedom):
2297            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2298            break;            break;
2299         case(ReducedContactElementsOne):     case(Elements):
2300            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);     case(ReducedElements):
2301              tags=mesh->Elements->tagsInUse;
2302              break;
2303       case(FaceElements):
2304       case(ReducedFaceElements):
2305              tags=mesh->FaceElements->tagsInUse;
2306              break;
2307       case(Points):
2308              tags=mesh->Points->tagsInUse;
2309              break;
2310       case(ContactElementsZero):
2311       case(ReducedContactElementsZero):
2312       case(ContactElementsOne):
2313       case(ReducedContactElementsOne):
2314              tags=mesh->ContactElements->tagsInUse;
2315            break;            break;
2316         default:     default:
2317            stringstream temp;        stringstream temp;
2318            temp << "Error - Finley does not know anything about function space type " << functionSpaceType;        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2319            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2320    }    }
2321    checkFinleyError();    return tags;
   return;  
2322  }  }
2323    
2324  void MeshAdapter::setTagMap(const std::string& name,  int tag)  
2325    bool MeshAdapter::canTag(int functionSpaceCode) const
2326  {  {
2327    Finley_Mesh* mesh=m_finleyMesh.get();    switch(functionSpaceCode) {
2328    Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     case(Nodes):
2329    checkPasoError();     case(Elements):
2330    // throwStandardException("MeshAdapter::set TagMap is not implemented.");     case(ReducedElements):
2331       case(FaceElements):
2332       case(ReducedFaceElements):
2333       case(Points):
2334       case(ContactElementsZero):
2335       case(ReducedContactElementsZero):
2336       case(ContactElementsOne):
2337       case(ReducedContactElementsOne):
2338              return true;
2339       case(ReducedNodes):
2340       case(DegreesOfFreedom):
2341       case(ReducedDegreesOfFreedom):
2342          return false;
2343       default:
2344        return false;
2345      }
2346  }  }
2347    
2348  int MeshAdapter::getTag(const std::string& name) const  AbstractDomain::StatusType MeshAdapter::getStatus() const
2349  {  {
2350    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2351    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;  
2352  }  }
2353    
2354  bool MeshAdapter::isValidTagName(const std::string& name) const  int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2355  {  {
2356      
2357    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2358    return Finley_Mesh_isValidTagName(mesh,name.c_str());    int order =-1;
2359      switch(functionSpaceCode) {
2360       case(Nodes):
2361       case(DegreesOfFreedom):
2362              order=mesh->approximationOrder;
2363              break;
2364       case(ReducedNodes):
2365       case(ReducedDegreesOfFreedom):
2366              order=mesh->reducedApproximationOrder;
2367              break;
2368       case(Elements):
2369       case(FaceElements):
2370       case(Points):
2371       case(ContactElementsZero):
2372       case(ContactElementsOne):
2373              order=mesh->integrationOrder;
2374              break;
2375       case(ReducedElements):
2376       case(ReducedFaceElements):
2377       case(ReducedContactElementsZero):
2378       case(ReducedContactElementsOne):
2379              order=mesh->reducedIntegrationOrder;
2380              break;
2381       default:
2382          stringstream temp;
2383          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2384          throw FinleyAdapterException(temp.str());
2385      }
2386      return order;
2387  }  }
2388    
2389  std::string MeshAdapter::showTagNames() const  bool MeshAdapter::supportsContactElements() const
2390  {  {
2391    stringstream temp;    return true;
2392    Finley_Mesh* mesh=m_finleyMesh.get();  }
2393    Finley_TagMap* tag_map=mesh->TagMap;  
2394    while (tag_map) {  ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2395       temp << tag_map->name;  {
2396       tag_map=tag_map->next;    m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2397       if (tag_map) temp << ", ";  }
2398    }  
2399    return temp.str();  ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2400    {
2401      Finley_ReferenceElementSet_dealloc(m_refSet);
2402  }  }
2403    
2404    // points will be flattened
2405    void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2406    {
2407          const int dim = getDim();
2408          int numPoints=points.size()/dim;
2409          int numTags=tags.size();
2410          Finley_Mesh* mesh=m_finleyMesh.get();
2411          
2412          if ( points.size() % dim != 0 )
2413          {
2414        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2415          }
2416          
2417          if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2418         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2419          
2420          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2421          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2422          
2423          for (int i=0;i<numPoints;++i) {
2424           points_ptr[ i * dim     ] = points[i * dim ];
2425           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2426           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2427               tags_ptr[i]=tags[i];
2428          }
2429          
2430          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2431          checkFinleyError();
2432          
2433          TMPMEMFREE(points_ptr);
2434          TMPMEMFREE(tags_ptr);
2435    }
2436    
2437    
2438    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2439    // {
2440    //       const int dim = getDim();
2441    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2442    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2443    //       Finley_Mesh* mesh=m_finleyMesh.get();
2444    //      
2445    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2446    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2447    //      
2448    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2449    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2450    //      
2451    //       for (int i=0;i<numPoints;++i) {
2452    //     int tag_id=-1;
2453    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2454    //     if  ( numComps !=   dim ) {
2455    //                stringstream temp;            
2456    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2457    //                throw FinleyAdapterException(temp.str());
2458    //     }
2459    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2460    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2461    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2462    //    
2463    //     if ( numTags > 0) {
2464    //            boost::python::extract<string> ex_str(tags[i]);
2465    //        if  ( ex_str.check() ) {
2466    //            tag_id=getTag( ex_str());
2467    //        } else {
2468    //             boost::python::extract<int> ex_int(tags[i]);
2469    //             if ( ex_int.check() ) {
2470    //                 tag_id=ex_int();
2471    //             } else {
2472    //              stringstream temp;          
2473    //                  temp << "Error - unable to extract tag for point " << i;
2474    //              throw FinleyAdapterException(temp.str());
2475    //            }
2476    //        }
2477    //     }      
2478    //            tags_ptr[i]=tag_id;
2479    //       }
2480    //      
2481    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2482    //       checkPasoError();
2483    //      
2484    //       TMPMEMFREE(points_ptr);
2485    //       TMPMEMFREE(tags_ptr);
2486    // }
2487    
2488    /*
2489    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2490    {  
2491        boost::python::list points =  boost::python::list();
2492        boost::python::list tags =  boost::python::list();
2493        points.append(point);
2494        tags.append(tag);
2495        addDiracPoints(points, tags);
2496    }
2497    */
2498    
2499    /*
2500    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2501    {
2502            boost::python::list points =   boost::python::list();
2503            boost::python::list tags =   boost::python::list();
2504            points.append(point);
2505            tags.append(tag);
2506            addDiracPoints(points, tags);
2507    }
2508    */
2509  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1407  
changed lines
  Added in v.4346

  ViewVC Help
Powered by ViewVC 1.1.26