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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26