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

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

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

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

Legend:
Removed from v.1344  
changed lines
  Added in v.2641

  ViewVC Help
Powered by ViewVC 1.1.26