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

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

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

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

Legend:
Removed from v.1366  
changed lines
  Added in v.3379

  ViewVC Help
Powered by ViewVC 1.1.26