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

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

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

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

Legend:
Removed from v.1388  
changed lines
  Added in v.3981

  ViewVC Help
Powered by ViewVC 1.1.26