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

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

  ViewVC Help
Powered by ViewVC 1.1.26