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

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

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

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

Legend:
Removed from v.1407  
changed lines
  Added in v.2635

  ViewVC Help
Powered by ViewVC 1.1.26