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

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

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

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

Legend:
Removed from v.1375  
changed lines
  Added in v.2898

  ViewVC Help
Powered by ViewVC 1.1.26