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

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

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

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

Legend:
Removed from v.1361  
changed lines
  Added in v.3355

  ViewVC Help
Powered by ViewVC 1.1.26