/[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 1366 by gross, Tue Dec 18 05:49:17 2007 UTC revision 2748 by gross, Tue Nov 17 07:32:59 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[12];     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
136     NcVar *ids, *data;     NcVar *ids;
137     int *int_ptr;     int *int_ptr;
138     Finley_Mesh *mesh = m_finleyMesh.get();     Finley_Mesh *mesh = m_finleyMesh.get();
139     Finley_TagMap* tag_map;     Finley_TagMap* tag_map;
# Line 121  void MeshAdapter::dump(const std::string Line 149  void MeshAdapter::dump(const std::string
149     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
150     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
151     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
152     char *newFileName = Paso_MPI_appendRankToFileName(strdup(fileName.c_str()), mpi_size, mpi_rank);  #ifdef PASO_MPI
153       MPI_Status status;
154    #endif
155    
156    /* Incoming token indicates it's my turn to write */
157    #ifdef PASO_MPI
158       if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
159    #endif
160    
161       char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),
162                                                         mpi_size, mpi_rank);
163    
164     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
165     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
166       num_Tags = 0;
167     if (tag_map) {     if (tag_map) {
168       while (tag_map) {        while (tag_map) {
169      num_Tags++;           num_Tags++;
170          tag_map=tag_map->next;           tag_map=tag_map->next;
171       }        }
172     }     }
173    
174     // NetCDF error handler     // NetCDF error handler
# Line 138  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)     if (num_Tags>0)
211        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
212          throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);
213    
214     // Attributes: MPI size, MPI rank, Name, order, reduced_order     // Attributes: MPI size, MPI rank, Name, order, reduced_order
215     if (!dataFile.add_att("mpi_size", mpi_size) )     if (!dataFile.add_att("mpi_size", mpi_size) )
216          throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);        throw DataException("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) )     if (!dataFile.add_att("numDim",numDim) )
222          throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        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) )     if (!dataFile.add_att("numNodes",numNodes) )
228          throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);        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) )     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
238          throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);
239     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
240          throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);
241     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
242          throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);
243     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
244        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);
245     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
246        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);
247     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
248        throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);
249     if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
250        throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);
251     if (!dataFile.add_att("num_Tags", num_Tags) )     if (!dataFile.add_att("num_Tags", num_Tags) )
252        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);
# Line 217  void MeshAdapter::dump(const std::string Line 256  void MeshAdapter::dump(const std::string
256     // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)     // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)
257     if (numNodes>0) {     if (numNodes>0) {
258    
259     // Nodes Id        // Nodes Id
260     if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
261          throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);
262     int_ptr = &mesh->Nodes->Id[0];        int_ptr = &mesh->Nodes->Id[0];
263     if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
264          throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);
265    
266     // Nodes Tag        // Nodes Tag
267     if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
268          throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);
269     int_ptr = &mesh->Nodes->Tag[0];        int_ptr = &mesh->Nodes->Tag[0];
270     if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
271          throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);
272    
273     // Nodes gDOF        // Nodes gDOF
274     if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
275          throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);
276     int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
277     if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
278          throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);
279    
280     // Nodes global node index        // Nodes global node index
281     if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
282          throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);
283     int_ptr = &mesh->Nodes->globalNodesIndex[0];        int_ptr = &mesh->Nodes->globalNodesIndex[0];
284     if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
285          throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);
286    
287     // Nodes grDof        // Nodes grDof
288     if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
289          throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);
290     int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
291     if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
292          throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);
293    
294     // Nodes grNI        // Nodes grNI
295     if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
296          throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);
297     int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
298     if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
299          throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);
300    
301     // Nodes Coordinates        // Nodes Coordinates
302     if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )        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);           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)) )        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);           throw DataException("Error - MeshAdapter::dump: copy Nodes_Coordinates to netCDF buffer failed: " + *newFileName);
306    
307     // Nodes degreesOfFreedomDistribution        // Nodes degreesOfFreedomDistribution
308     if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )        if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
309          throw DataException("Error - MeshAdapter::dump: appending Nodes_DofDistribution to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Nodes_DofDistribution to netCDF file failed: " + *newFileName);
310     int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];        int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
311     if (! (ids->put(int_ptr, mpi_size+1)) )        if (! (ids->put(int_ptr, mpi_size+1)) )
312          throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);
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     }     }
322    
# Line 278  void MeshAdapter::dump(const std::string Line 324  void MeshAdapter::dump(const std::string
324    
325     if (num_Elements>0) {     if (num_Elements>0) {
326    
327       // Temp storage to gather node IDs        // Elements_Id
328       int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);        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       // Elements_Id        int_ptr = &mesh->Elements->Id[0];
331       if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! (ids->put(int_ptr, num_Elements)) )
332          throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);
333       int_ptr = &mesh->Elements->Id[0];  
334       if (! (ids->put(int_ptr, num_Elements)) )        // Elements_Tag
335          throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);        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       // Elements_Tag        int_ptr = &mesh->Elements->Tag[0];
338       if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )        if (! (ids->put(int_ptr, num_Elements)) )
339          throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);
340       int_ptr = &mesh->Elements->Tag[0];  
341       if (! (ids->put(int_ptr, num_Elements)) )        // Elements_Owner
342          throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);        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       // Elements_Owner        int_ptr = &mesh->Elements->Owner[0];
345       if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )        if (! (ids->put(int_ptr, num_Elements)) )
346          throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);
347       int_ptr = &mesh->Elements->Owner[0];  
348       if (! (ids->put(int_ptr, num_Elements)) )        // Elements_Color
349          throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);        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       // Elements_Color        int_ptr = &mesh->Elements->Color[0];
352       if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )        if (! (ids->put(int_ptr, num_Elements)) )
353          throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);
354       int_ptr = &mesh->Elements->Color[0];  
355       if (! (ids->put(int_ptr, num_Elements)) )        // Elements_Nodes
356          throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
357             throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);
358       // Elements_Nodes        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
359       for (int i=0; i<num_Elements; i++)           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);
        for (int j=0; j<num_Elements_numNodes; j++)  
          Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)] = mesh->Nodes->Id[mesh->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]];  
      if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )  
         throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);  
      if (! (ids->put(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes)) )  
         throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);  
   
      TMPMEMFREE(Elements_Nodes);  
360    
361     }     }
362    
# Line 326  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_Id        int_ptr = &mesh->FaceElements->Id[0];
371       if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! (ids->put(int_ptr, num_FaceElements)) )
372          throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);
373       int_ptr = &mesh->FaceElements->Id[0];  
374       if (! (ids->put(int_ptr, num_FaceElements)) )        // FaceElements_Tag
375          throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
376             throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);
377       // FaceElements_Tag        int_ptr = &mesh->FaceElements->Tag[0];
378       if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )        if (! (ids->put(int_ptr, num_FaceElements)) )
379          throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);
380       int_ptr = &mesh->FaceElements->Tag[0];  
381       if (! (ids->put(int_ptr, num_FaceElements)) )        // FaceElements_Owner
382          throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
383             throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);
384       // FaceElements_Owner        int_ptr = &mesh->FaceElements->Owner[0];
385       if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )        if (! (ids->put(int_ptr, num_FaceElements)) )
386          throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);
387       int_ptr = &mesh->FaceElements->Owner[0];  
388       if (! (ids->put(int_ptr, num_FaceElements)) )        // FaceElements_Color
389          throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
390             throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);
391       // FaceElements_Color        int_ptr = &mesh->FaceElements->Color[0];
392       if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )        if (! (ids->put(int_ptr, num_FaceElements)) )
393          throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);
394       int_ptr = &mesh->FaceElements->Color[0];  
395       if (! (ids->put(int_ptr, num_FaceElements)) )        // FaceElements_Nodes
396          throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
397             throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);
398       // FaceElements_Nodes        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
399       for (int i=0; i<num_FaceElements; i++)           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);
        for (int j=0; j<num_FaceElements_numNodes; j++)  
          FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = mesh->Nodes->Id[mesh->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]];  
      if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )  
         throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);  
      if (! (ids->put(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )  
         throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);  
   
      TMPMEMFREE(FaceElements_Nodes);  
400    
401     }     }
402    
# Line 374  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_Id        int_ptr = &mesh->ContactElements->Id[0];
411       if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )        if (! (ids->put(int_ptr, num_ContactElements)) )
412          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);
413       int_ptr = &mesh->ContactElements->Id[0];  
414       if (! (ids->put(int_ptr, num_ContactElements)) )        // ContactElements_Tag
415          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
416             throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);
417       // ContactElements_Tag        int_ptr = &mesh->ContactElements->Tag[0];
418       if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )        if (! (ids->put(int_ptr, num_ContactElements)) )
419          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);
420       int_ptr = &mesh->ContactElements->Tag[0];  
421       if (! (ids->put(int_ptr, num_ContactElements)) )        // ContactElements_Owner
422          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
423             throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);
424       // ContactElements_Owner        int_ptr = &mesh->ContactElements->Owner[0];
425       if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )        if (! (ids->put(int_ptr, num_ContactElements)) )
426          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);
427       int_ptr = &mesh->ContactElements->Owner[0];  
428       if (! (ids->put(int_ptr, num_ContactElements)) )        // ContactElements_Color
429          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
430             throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);
431       // ContactElements_Color        int_ptr = &mesh->ContactElements->Color[0];
432       if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )        if (! (ids->put(int_ptr, num_ContactElements)) )
433          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);
434       int_ptr = &mesh->ContactElements->Color[0];  
435       if (! (ids->put(int_ptr, num_ContactElements)) )        // ContactElements_Nodes
436          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
437             throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);
438       // ContactElements_Nodes        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
439       for (int i=0; i<num_ContactElements; i++)           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);
        for (int j=0; j<num_ContactElements_numNodes; j++)  
          ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)] = mesh->Nodes->Id[mesh->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]];  
      if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )  
         throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);  
      if (! (ids->put(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )  
         throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);  
   
      TMPMEMFREE(ContactElements_Nodes);  
440    
441     }     }
442    
# Line 422  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);  
448    
449       // Points_Id        // Points_Id
450       if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
451          throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);
452       int_ptr = &mesh->Points->Id[0];        int_ptr = &mesh->Points->Id[0];
453       if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
454          throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);
455    
456       // Points_Tag        // Points_Tag
457       if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
458          throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);
459       int_ptr = &mesh->Points->Tag[0];        int_ptr = &mesh->Points->Tag[0];
460       if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
461          throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);
462    
463       // Points_Owner        // Points_Owner
464       if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
465          throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);
466       int_ptr = &mesh->Points->Owner[0];        int_ptr = &mesh->Points->Owner[0];
467       if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
468          throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);
469    
470       // Points_Color        // Points_Color
471       if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
472          throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);
473       int_ptr = &mesh->Points->Color[0];        int_ptr = &mesh->Points->Color[0];
474       if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
475          throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);
476    
477       // Points_Nodes        // Points_Nodes
478       // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
479       for (int i=0; i<num_Points; i++)        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
480         Points_Nodes[i] = mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]];           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);
481       if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
482          throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);
      if (! (ids->put(&(Points_Nodes[0]), num_Points)) )  
         throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);  
   
      TMPMEMFREE(Points_Nodes);  
483    
484     }     }
485    
# Line 472  void MeshAdapter::dump(const std::string Line 487  void MeshAdapter::dump(const std::string
487    
488     if (num_Tags>0) {     if (num_Tags>0) {
489    
490       // Temp storage to gather node IDs        // Temp storage to gather node IDs
491       int *Tags_keys = TMPMEMALLOC(num_Tags, int);        int *Tags_keys = TMPMEMALLOC(num_Tags, int);
492       char name_temp[4096];        char name_temp[4096];
493    
494       /* Copy tag data into temp arrays */        /* Copy tag data into temp arrays */
495       tag_map = mesh->TagMap;        tag_map = mesh->TagMap;
496       if (tag_map) {        if (tag_map) {
497         int i = 0;           int i = 0;
498         while (tag_map) {           while (tag_map) {
499      Tags_keys[i++] = tag_map->tag_key;              Tags_keys[i++] = tag_map->tag_key;
500          tag_map=tag_map->next;              tag_map=tag_map->next;
501         }           }
502       }        }
503    
504       // Tags_keys        // Tags_keys
505       if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
506          throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);
507       int_ptr = &Tags_keys[0];        int_ptr = &Tags_keys[0];
508       if (! (ids->put(int_ptr, num_Tags)) )        if (! (ids->put(int_ptr, num_Tags)) )
509          throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);
510    
511       // Tags_names_*        // 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        // 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        // because the NetCDF manual doesn't tell how to do an array of strings
514       tag_map = mesh->TagMap;        tag_map = mesh->TagMap;
515       if (tag_map) {        if (tag_map) {
516         int i = 0;           int i = 0;
517         while (tag_map) {           while (tag_map) {
518           sprintf(name_temp, "Tags_name_%d", i);              sprintf(name_temp, "Tags_name_%d", i);
519           if (!dataFile.add_att(name_temp, tag_map->name) )              if (!dataFile.add_att(name_temp, tag_map->name) )
520             throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);                 throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);
521           tag_map=tag_map->next;              tag_map=tag_map->next;
522       i++;              i++;
523         }           }
524       }        }
525    
526       TMPMEMFREE(Tags_keys);        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 522  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 638  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 653  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->referenceElementSet->referenceElement->Parametrization->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->referenceElementSet->referenceElementReducedQuadrature->Parametrization->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->referenceElementSet->referenceElement->Parametrization->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->referenceElementSet->referenceElementReducedQuadrature->Parametrization->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->referenceElementSet->referenceElement->Parametrization->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->referenceElementSet->referenceElementReducedQuadrature->Parametrization->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->referenceElementSet->referenceElement->Parametrization->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->referenceElementSet->referenceElementReducedQuadrature->Parametrization->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 770  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 809  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  cout << "A\n";        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
993          if (target.getFunctionSpace().getTypeCode()==FaceElements) {           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
994  cout << "B\n";        } else {
995             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
996  cout << "C\n";        }
997          throw FinleyAdapterException("A");        break;
998          } else {     case(Points):
999             throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");        if (target.getFunctionSpace().getTypeCode()==Points) {
1000         }           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1001  cout << "A\n";        } else {
1002         break;           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1003       case(Points):        }
1004          if (target.getFunctionSpace().getTypeCode()==Points) {        break;
1005             Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);     case(ContactElementsZero):
1006          } else {     case(ContactElementsOne):
1007             throw FinleyAdapterException("Error - No interpolation with data on points possible.");        if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1008          }           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1009          break;        } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1010       case(ContactElementsZero):           Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1011       case(ContactElementsOne):        } else {
1012          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {           throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1013             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);        }
1014          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {        break;
1015             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);     case(ReducedContactElementsZero):
1016          } else {     case(ReducedContactElementsOne):
1017             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1018          }           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1019          break;        } else {
1020       case(ReducedContactElementsZero):           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1021       case(ReducedContactElementsOne):        }
1022          if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {        break;
1023             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);     case(DegreesOfFreedom):      
1024          } else {        switch(target.getFunctionSpace().getTypeCode()) {
1025             throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");        case(ReducedDegreesOfFreedom):
1026          }        case(DegreesOfFreedom):
1027          break;        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1028       case(DegreesOfFreedom):              break;
1029          switch(target.getFunctionSpace().getTypeCode()) {    
1030             case(ReducedDegreesOfFreedom):        case(Nodes):
1031             case(DegreesOfFreedom):        case(ReducedNodes):
1032                Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        if (getMPISize()>1) {
1033                break;           escript::Data temp=escript::Data(in);
1034             temp.expand();
1035             case(Nodes):           escriptDataC _in2 = temp.getDataC();
1036             case(ReducedNodes):           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1037                if (getMPISize()>1) {        } else {
1038                    escript::Data temp=escript::Data(in);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1039                    temp.expand();        }
1040                    escriptDataC _in2 = temp.getDataC();        break;
1041                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);        case(Elements):
1042                } else {        case(ReducedElements):
1043                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        if (getMPISize()>1) {
1044                }           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1045                break;           escriptDataC _in2 = temp.getDataC();
1046             case(Elements):           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1047             case(ReducedElements):        } else {
1048                if (getMPISize()>1) {           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1049                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );        }
1050                   escriptDataC _in2 = temp.getDataC();        break;
1051                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);        case(FaceElements):
1052                } else {        case(ReducedFaceElements):
1053                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        if (getMPISize()>1) {
1054                }           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1055                break;           escriptDataC _in2 = temp.getDataC();
1056             case(FaceElements):           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1057             case(ReducedFaceElements):    
1058                if (getMPISize()>1) {        } else {
1059                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1060                   escriptDataC _in2 = temp.getDataC();        }
1061                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);        break;
1062          case(Points):
1063                } else {        if (getMPISize()>1) {
1064                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1065                }           escriptDataC _in2 = temp.getDataC();
1066                break;        } else {
1067             case(Points):           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1068                if (getMPISize()>1) {        }
1069                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );        break;
1070                   escriptDataC _in2 = temp.getDataC();        case(ContactElementsZero):
1071                } else {        case(ContactElementsOne):
1072                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        case(ReducedContactElementsZero):
1073                }        case(ReducedContactElementsOne):
1074                break;        if (getMPISize()>1) {
1075             case(ContactElementsZero):           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1076             case(ContactElementsOne):           escriptDataC _in2 = temp.getDataC();
1077             case(ReducedContactElementsZero):           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1078             case(ReducedContactElementsOne):        } else {
1079                if (getMPISize()>1) {           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1080                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );        }
1081                   escriptDataC _in2 = temp.getDataC();        break;
1082                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);        default:
1083                } else {           stringstream temp;
1084                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1085                }           throw FinleyAdapterException(temp.str());
1086               break;           break;
1087             default:        }
1088               stringstream temp;        break;
1089               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();     case(ReducedDegreesOfFreedom):
1090               throw FinleyAdapterException(temp.str());        switch(target.getFunctionSpace().getTypeCode()) {
1091               break;        case(Nodes):
1092          }        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1093          break;        break;
1094       case(ReducedDegreesOfFreedom):        case(ReducedNodes):
1095         switch(target.getFunctionSpace().getTypeCode()) {        if (getMPISize()>1) {
1096            case(Nodes):           escript::Data temp=escript::Data(in);
1097               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           temp.expand();
1098               break;           escriptDataC _in2 = temp.getDataC();
1099            case(ReducedNodes):           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1100                if (getMPISize()>1) {        } else {
1101                    escript::Data temp=escript::Data(in);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1102                    temp.expand();        }
1103                    escriptDataC _in2 = temp.getDataC();        break;
1104                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);        case(DegreesOfFreedom):
1105                } else {        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1106                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        break;
1107                }        case(ReducedDegreesOfFreedom):
1108                break;        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1109            case(DegreesOfFreedom):        break;
1110               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");        case(Elements):
1111               break;        case(ReducedElements):
1112            case(ReducedDegreesOfFreedom):        if (getMPISize()>1) {
1113               Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1114               break;           escriptDataC _in2 = temp.getDataC();
1115            case(Elements):           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1116            case(ReducedElements):        } else {
1117                if (getMPISize()>1) {           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1118                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );        }
1119                   escriptDataC _in2 = temp.getDataC();        break;
1120                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);        case(FaceElements):
1121                } else {        case(ReducedFaceElements):
1122                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        if (getMPISize()>1) {
1123               }           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1124               break;           escriptDataC _in2 = temp.getDataC();
1125            case(FaceElements):           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1126            case(ReducedFaceElements):        } else {
1127                if (getMPISize()>1) {           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1128                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );        }
1129                   escriptDataC _in2 = temp.getDataC();        break;
1130                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);        case(Points):
1131                } else {        if (getMPISize()>1) {
1132                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1133                }           escriptDataC _in2 = temp.getDataC();
1134               break;           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1135            case(Points):        } else {
1136                if (getMPISize()>1) {           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1137                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );        }
1138                   escriptDataC _in2 = temp.getDataC();        break;
1139                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);        case(ContactElementsZero):
1140                } else {        case(ContactElementsOne):
1141                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        case(ReducedContactElementsZero):
1142               }        case(ReducedContactElementsOne):
1143               break;        if (getMPISize()>1) {
1144            case(ContactElementsZero):           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1145            case(ContactElementsOne):           escriptDataC _in2 = temp.getDataC();
1146            case(ReducedContactElementsZero):           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1147            case(ReducedContactElementsOne):        } else {
1148                if (getMPISize()>1) {           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1149                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );        }
1150                   escriptDataC _in2 = temp.getDataC();        break;
1151                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);        default:
1152                } else {           stringstream temp;
1153                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1154               }           throw FinleyAdapterException(temp.str());
1155               break;           break;
1156            default:        }
1157               stringstream temp;        break;
1158               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();     default:
1159               throw FinleyAdapterException(temp.str());        stringstream temp;
1160               break;        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1161         }        throw FinleyAdapterException(temp.str());
1162         break;        break;
1163       default:     }
1164          stringstream temp;     checkFinleyError();
         temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();  
         throw FinleyAdapterException(temp.str());  
         break;  
   }  
   checkFinleyError();  
1165  }  }
1166    
1167  //  //
# Line 1105  cout << "A\n"; Line 1169  cout << "A\n";
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 1128  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 1183  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 1195  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 1267  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 1345  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  // saves a data array in openDX format:  //
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    //
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;  
   }  
1542    
1543    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1544    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;     Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());
1545    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;     checkFinleyError();
   
     boost::python::list keys=arg.keys();  
     for (int i=0;i<num_data;++i) {  
          std::string n=boost::python::extract<std::string>(keys[i]);  
          escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);  
          if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)  
              throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");  
          data[i]=d.getDataC();  
          ptr_data[i]=&(data[i]);  
          c_names[i]=&(names[i][0]);  
          if (n.length()>MAX_namelength-1) {  
             strncpy(c_names[i],n.c_str(),MAX_namelength-1);  
             c_names[i][MAX_namelength-1]='\0';  
          } else {  
             strcpy(c_names[i],n.c_str());  
          }  
     }  
     Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);  
   
   checkFinleyError();  
   /* win32 refactor */  
   TMPMEMFREE(c_names);  
   TMPMEMFREE(data);  
   TMPMEMFREE(ptr_data);  
   for(int i=0;i<num_data;i++)  
   {  
     TMPMEMFREE(names[i]);  
   }  
   TMPMEMFREE(names);  
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                                                                                                                                                                            bool MeshAdapter::ownSample(int fs_code, index_t id) const
1558  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  {
1559  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  #ifdef PASO_MPI
1560                        const int row_blocksize,      index_t myFirstNode=0, myLastNode=0, k=0;
1561                        const escript::FunctionSpace& row_functionspace,      index_t* globalNodeIndex=0;
1562                        const int column_blocksize,      Finley_Mesh* mesh_p=m_finleyMesh.get();
1563                        const escript::FunctionSpace& column_functionspace,      if (fs_code == FINLEY_REDUCED_NODES)
1564                        const int type) const      {
1565  {      myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1566      int reduceRowOrder=0;      myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1567      int reduceColOrder=0;      globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
     // is the domain right?  
     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());  
     if (row_domain!=*this)  
           throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");  
     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());  
     if (col_domain!=*this)  
           throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");  
     // is the function space type right  
     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {  
         reduceRowOrder=0;  
     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {  
         reduceRowOrder=1;  
     } else {  
         throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");  
1568      }      }
1569      if (column_functionspace.getTypeCode()==DegreesOfFreedom) {      else
1570          reduceColOrder=0;      {
1571      } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1572          reduceColOrder=1;      myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1573      } else {      globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
         throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");  
1574      }      }
1575      // generate matrix:      k=globalNodeIndex[id];
1576            return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1577      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);  #endif
1578      checkFinleyError();      return true;
1579      Paso_SystemMatrix* fsystemMatrix;  }
1580      int trilinos = 0;  
1581      if (trilinos) {  
1582    
1583    //
1584    // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1585    //
1586    SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1587                                                     const int row_blocksize,
1588                                                     const escript::FunctionSpace& row_functionspace,
1589                                                     const int column_blocksize,
1590                                                     const escript::FunctionSpace& column_functionspace,
1591                                                     const int type) const
1592    {
1593       int reduceRowOrder=0;
1594       int reduceColOrder=0;
1595       // is the domain right?
1596       const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1597       if (row_domain!=*this)
1598          throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1599       const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1600       if (col_domain!=*this)
1601          throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1602       // is the function space type right
1603       if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1604          reduceRowOrder=0;
1605       } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1606          reduceRowOrder=1;
1607       } else {
1608          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1609       }
1610       if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1611          reduceColOrder=0;
1612       } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1613          reduceColOrder=1;
1614       } else {
1615          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1616       }
1617       // generate matrix:
1618    
1619       Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1620       checkFinleyError();
1621       Paso_SystemMatrix* fsystemMatrix;
1622       int trilinos = 0;
1623       if (trilinos) {
1624  #ifdef TRILINOS  #ifdef TRILINOS
1625        /* Allocation Epetra_VrbMatrix here */        /* Allocation Epetra_VrbMatrix here */
1626  #endif  #endif
1627      }     }
1628      else {     else {
1629        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1630      }     }
1631      checkPasoError();     checkPasoError();
1632      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1633      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1634  }  }
1635    
1636    //
1637  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1638    //
1639  TransportProblemAdapter MeshAdapter::newTransportProblem(  TransportProblemAdapter MeshAdapter::newTransportProblem(
1640                        const double theta,                                                           const double theta,
1641                        const double dt_max,                                                           const int blocksize,
1642                        const int blocksize,                                                           const escript::FunctionSpace& functionspace,
1643                        const escript::FunctionSpace& functionspace,                                                           const int type) const
1644                        const int type) const  {
1645  {     int reduceOrder=0;
1646      int reduceOrder=0;     // is the domain right?
1647      // is the domain right?     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1648      const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(functionspace.getDomain());     if (domain!=*this)
1649      if (domain!=*this)        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1650            throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");     // is the function space type right
1651      // is the function space type right     if (functionspace.getTypeCode()==DegreesOfFreedom) {
1652      if (functionspace.getTypeCode()==DegreesOfFreedom) {        reduceOrder=0;
1653          reduceOrder=0;     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1654      } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {        reduceOrder=1;
1655          reduceOrder=1;     } else {
1656      } else {        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1657          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");     }
1658      }     // generate matrix:
1659      // generate matrix:  
1660           Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1661      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     checkFinleyError();
1662      checkFinleyError();     Paso_FCTransportProblem* transportProblem;
1663      Paso_FCTransportProblem* transportProblem;     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);
1664      transportProblem=Paso_FCTransportProblem_alloc(theta,dt_max,fsystemMatrixPattern,blocksize);     checkPasoError();
1665      checkPasoError();     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1666      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);
     return TransportProblemAdapter(transportProblem,theta,dt_max,blocksize,functionspace);  
1667  }  }
1668    
1669  //  //
# Line 1593  TransportProblemAdapter MeshAdapter::new Line 1675  TransportProblemAdapter MeshAdapter::new
1675  // 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:
1676  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1677  {  {
1678    switch(functionSpaceCode) {     switch(functionSpaceCode) {
1679         case(Nodes):     case(Nodes):
1680         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1681         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1682            return false;     return false;
1683            break;     break;
1684         case(Elements):     case(Elements):
1685         case(FaceElements):     case(FaceElements):
1686         case(Points):     case(Points):
1687         case(ContactElementsZero):     case(ContactElementsZero):
1688         case(ContactElementsOne):     case(ContactElementsOne):
1689         case(ReducedElements):     case(ReducedElements):
1690         case(ReducedFaceElements):     case(ReducedFaceElements):
1691         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1692         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1693            return true;     return true;
1694            break;     break;
1695         default:     default:
1696            stringstream temp;        stringstream temp;
1697            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;
1698            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1699            break;        break;
1700    }     }
1701    checkFinleyError();     checkFinleyError();
1702    return false;     return false;
1703    }
1704    
1705    bool
1706    MeshAdapter::commonFunctionSpace(const std::vector<int>& fs, int& resultcode) const
1707    {
1708       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1709        class 1: DOF <-> Nodes
1710        class 2: ReducedDOF <-> ReducedNodes
1711        class 3: Points
1712        class 4: Elements
1713        class 5: ReducedElements
1714        class 6: FaceElements
1715        class 7: ReducedFaceElements
1716        class 8: ContactElementZero <-> ContactElementOne
1717        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1718    
1719       There is also a set of lines. Interpolation is possible down a line but not between lines.
1720       class 1 and 2 belong to all lines so aren't considered.
1721        line 0: class 3
1722        line 1: class 4,5
1723        line 2: class 6,7
1724        line 3: class 8,9
1725    
1726       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1727       eg hasnodes is true if we have at least one instance of Nodes.
1728       */
1729        if (fs.size()==0)
1730        {
1731        return false;
1732        }
1733        std::vector<int> hasclass(10);
1734        std::vector<int> hasline(4);    
1735        bool hasnodes=false;
1736        bool hasrednodes=false;
1737        bool hascez=false;
1738        bool hasrcez=false;
1739        for (int i=0;i<fs.size();++i)
1740        {
1741        switch(fs[i])
1742        {
1743        case(Nodes):   hasnodes=true;   // no break is deliberate
1744        case(DegreesOfFreedom):
1745            hasclass[1]=1;
1746            break;
1747        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1748        case(ReducedDegreesOfFreedom):
1749            hasclass[2]=1;
1750            break;
1751        case(Points):
1752            hasline[0]=1;
1753            hasclass[3]=1;
1754            break;
1755        case(Elements):
1756            hasclass[4]=1;
1757            hasline[1]=1;
1758            break;
1759        case(ReducedElements):
1760            hasclass[5]=1;
1761            hasline[1]=1;
1762            break;
1763        case(FaceElements):
1764            hasclass[6]=1;
1765            hasline[2]=1;
1766            break;
1767        case(ReducedFaceElements):
1768            hasclass[7]=1;
1769            hasline[2]=1;
1770            break;
1771        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1772        case(ContactElementsOne):
1773            hasclass[8]=1;
1774            hasline[3]=1;
1775            break;
1776        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1777        case(ReducedContactElementsOne):
1778            hasclass[9]=1;
1779            hasline[3]=1;
1780            break;
1781        default:
1782            return false;
1783        }
1784        }
1785        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1786        // fail if we have more than one leaf group
1787    
1788        if (totlines>1)
1789        {
1790        return false;   // there are at least two branches we can't interpolate between
1791        }
1792        else if (totlines==1)
1793        {
1794        if (hasline[0]==1)      // we have points
1795        {
1796            resultcode=Points;
1797        }
1798        else if (hasline[1]==1)
1799        {
1800            if (hasclass[5]==1)
1801            {
1802            resultcode=ReducedElements;
1803            }
1804            else
1805            {
1806            resultcode=Elements;
1807            }
1808        }
1809        else if (hasline[2]==1)
1810        {
1811            if (hasclass[7]==1)
1812            {
1813            resultcode=ReducedFaceElements;
1814            }
1815            else
1816            {
1817            resultcode=FaceElements;
1818            }
1819        }
1820        else    // so we must be in line3
1821        {
1822            if (hasclass[9]==1)
1823            {
1824            // need something from class 9
1825            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1826            }
1827            else
1828            {
1829            // something from class 8
1830            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1831            }
1832        }
1833        }
1834        else    // totlines==0
1835        {
1836        if (hasclass[2]==1)
1837        {
1838            // something from class 2
1839            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1840        }
1841        else
1842        {   // something from class 1
1843            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1844        }
1845        }
1846        return true;
1847  }  }
1848    
1849  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1850  {  {
1851    switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1852       case(Nodes):     case(Nodes):
1853          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1854             case(Nodes):      case(Nodes):
1855             case(ReducedNodes):      case(ReducedNodes):
1856             case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1857             case(DegreesOfFreedom):      case(DegreesOfFreedom):
1858             case(Elements):      case(Elements):
1859             case(ReducedElements):      case(ReducedElements):
1860             case(FaceElements):      case(FaceElements):
1861             case(ReducedFaceElements):      case(ReducedFaceElements):
1862             case(Points):      case(Points):
1863             case(ContactElementsZero):      case(ContactElementsZero):
1864             case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1865             case(ContactElementsOne):      case(ContactElementsOne):
1866             case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1867                 return true;      return true;
1868             default:      default:
1869                 stringstream temp;            stringstream temp;
1870                 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;
1871                 throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1872          }     }
1873          break;     break;
1874       case(ReducedNodes):     case(ReducedNodes):
1875          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1876             case(ReducedNodes):      case(ReducedNodes):
1877             case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1878             case(Elements):      case(Elements):
1879             case(ReducedElements):      case(ReducedElements):
1880             case(FaceElements):      case(FaceElements):
1881             case(ReducedFaceElements):      case(ReducedFaceElements):
1882             case(Points):      case(Points):
1883             case(ContactElementsZero):      case(ContactElementsZero):
1884             case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1885             case(ContactElementsOne):      case(ContactElementsOne):
1886             case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1887                 return true;      return true;
1888            case(Nodes):      case(Nodes):
1889            case(DegreesOfFreedom):      case(DegreesOfFreedom):
1890               return false;      return false;
1891             default:      default:
1892                 stringstream temp;          stringstream temp;
1893                 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;
1894                 throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1895          }     }
1896          break;     break;
1897       case(Elements):     case(Elements):
1898          if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1899             return true;        return true;
1900          } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1901             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;  
1902          } else {          } else {
1903             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());  
1904          }          }
1905          break;     case(ReducedElements):
1906       case(ReducedDegreesOfFreedom):      if (functionSpaceType_target==ReducedElements) {
1907         switch(functionSpaceType_target) {        return true;
1908            case(ReducedDegreesOfFreedom):      } else {
1909            case(ReducedNodes):            return false;
1910            case(Elements):      }
1911            case(ReducedElements):     case(FaceElements):
1912            case(FaceElements):      if (functionSpaceType_target==FaceElements) {
1913            case(ReducedFaceElements):              return true;
1914            case(Points):      } else if (functionSpaceType_target==ReducedFaceElements) {
1915            case(ContactElementsZero):              return true;
1916            case(ReducedContactElementsZero):      } else {
1917            case(ContactElementsOne):              return false;
1918            case(ReducedContactElementsOne):      }
1919                return true;     case(ReducedFaceElements):
1920            case(Nodes):      if (functionSpaceType_target==ReducedFaceElements) {
1921            case(DegreesOfFreedom):              return true;
1922               return false;      } else {
1923            default:          return false;
1924               stringstream temp;      }
1925               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;     case(Points):
1926               throw FinleyAdapterException(temp.str());      if (functionSpaceType_target==Points) {
1927         }              return true;
1928         break;      } else {
1929       default:              return false;
1930          stringstream temp;      }
1931          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;     case(ContactElementsZero):
1932          throw FinleyAdapterException(temp.str());     case(ContactElementsOne):
1933          break;      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1934    }              return true;
1935    checkFinleyError();      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1936    return false;              return true;
1937        } else {
1938                return false;
1939        }
1940       case(ReducedContactElementsZero):
1941       case(ReducedContactElementsOne):
1942        if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1943                return true;
1944        } else {
1945                return false;
1946        }
1947       case(DegreesOfFreedom):
1948        switch(functionSpaceType_target) {
1949        case(ReducedDegreesOfFreedom):
1950        case(DegreesOfFreedom):
1951        case(Nodes):
1952        case(ReducedNodes):
1953        case(Elements):
1954        case(ReducedElements):
1955        case(Points):
1956        case(FaceElements):
1957        case(ReducedFaceElements):
1958        case(ContactElementsZero):
1959        case(ReducedContactElementsZero):
1960        case(ContactElementsOne):
1961        case(ReducedContactElementsOne):
1962        return true;
1963        default:
1964            stringstream temp;
1965            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1966            throw FinleyAdapterException(temp.str());
1967        }
1968        break;
1969       case(ReducedDegreesOfFreedom):
1970       switch(functionSpaceType_target) {
1971        case(ReducedDegreesOfFreedom):
1972        case(ReducedNodes):
1973        case(Elements):
1974        case(ReducedElements):
1975        case(FaceElements):
1976        case(ReducedFaceElements):
1977        case(Points):
1978        case(ContactElementsZero):
1979        case(ReducedContactElementsZero):
1980        case(ContactElementsOne):
1981        case(ReducedContactElementsOne):
1982        return true;
1983        case(Nodes):
1984        case(DegreesOfFreedom):
1985        return false;
1986        default:
1987            stringstream temp;
1988            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1989            throw FinleyAdapterException(temp.str());
1990        }
1991        break;
1992       default:
1993          stringstream temp;
1994          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1995          throw FinleyAdapterException(temp.str());
1996          break;
1997       }
1998       checkFinleyError();
1999       return false;
2000  }  }
2001    
2002  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
2003  {  {
2004      return false;     return false;
2005  }  }
2006    
2007  bool MeshAdapter::operator==(const AbstractDomain& other) const  bool MeshAdapter::operator==(const AbstractDomain& other) const
2008  {  {
2009    const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);     const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2010    if (temp!=0) {     if (temp!=0) {
2011      return (m_finleyMesh==temp->m_finleyMesh);        return (m_finleyMesh==temp->m_finleyMesh);
2012    } else {     } else {
2013      return false;        return false;
2014    }     }
2015  }  }
2016    
2017  bool MeshAdapter::operator!=(const AbstractDomain& other) const  bool MeshAdapter::operator!=(const AbstractDomain& other) const
2018  {  {
2019    return !(operator==(other));     return !(operator==(other));
2020  }  }
2021    
2022  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
2023  {  {
2024     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2025       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2026       checkPasoError();
2027       return out;
2028    }
2029    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2030    {
2031       Finley_Mesh* mesh=m_finleyMesh.get();
2032       int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2033     checkPasoError();     checkPasoError();
2034     return out;     return out;
2035  }  }
2036    
2037  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2038  {  {
2039    return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(asAbstractContinuousDomain()).getX();
2040  }  }
2041    
2042  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2043  {  {
2044    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
2045  }  }
2046    
2047  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2048  {  {
2049    return function(asAbstractContinuousDomain()).getSize();     return escript::function(asAbstractContinuousDomain()).getSize();
2050  }  }
2051    
2052  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2053  {  {
2054    int *out=0,i;     int *out = NULL;
2055    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2056    switch (functionSpaceType) {     switch (functionSpaceType) {
2057      case(Nodes):     case(Nodes):
2058        out=mesh->Nodes->Id;     out=mesh->Nodes->Id;
2059        break;     break;
2060      case(ReducedNodes):     case(ReducedNodes):
2061        out=mesh->Nodes->reducedNodesId;     out=mesh->Nodes->reducedNodesId;
2062        break;     break;
2063      case(Elements):     case(Elements):
2064        out=mesh->Elements->Id;     out=mesh->Elements->Id;
2065        break;     break;
2066      case(ReducedElements):     case(ReducedElements):
2067        out=mesh->Elements->Id;     out=mesh->Elements->Id;
2068        break;     break;
2069      case(FaceElements):     case(FaceElements):
2070        out=mesh->FaceElements->Id;     out=mesh->FaceElements->Id;
2071        break;     break;
2072      case(ReducedFaceElements):     case(ReducedFaceElements):
2073        out=mesh->FaceElements->Id;     out=mesh->FaceElements->Id;
2074        break;     break;
2075      case(Points):     case(Points):
2076        out=mesh->Points->Id;     out=mesh->Points->Id;
2077        break;     break;
2078      case(ContactElementsZero):     case(ContactElementsZero):
2079        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2080        break;     break;
2081      case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
2082        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2083        break;     break;
2084      case(ContactElementsOne):     case(ContactElementsOne):
2085        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2086        break;     break;
2087      case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
2088        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2089        break;     break;
2090      case(DegreesOfFreedom):     case(DegreesOfFreedom):
2091        out=mesh->Nodes->degreesOfFreedomId;     out=mesh->Nodes->degreesOfFreedomId;
2092        break;     break;
2093      case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2094        out=mesh->Nodes->reducedDegreesOfFreedomId;     out=mesh->Nodes->reducedDegreesOfFreedomId;
2095        break;     break;
2096      default:     default:
2097        stringstream temp;        stringstream temp;
2098        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2099        throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2100        break;        break;
2101    }     }
2102    return out;     return out;
2103  }  }
2104  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2105  {  {
2106    int out=0;     int out=0;
2107    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2108    switch (functionSpaceType) {     switch (functionSpaceType) {
2109    case(Nodes):     case(Nodes):
2110      out=mesh->Nodes->Tag[sampleNo];     out=mesh->Nodes->Tag[sampleNo];
2111      break;     break;
2112    case(ReducedNodes):     case(ReducedNodes):
2113      throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2114      break;     break;
2115    case(Elements):     case(Elements):
2116      out=mesh->Elements->Tag[sampleNo];     out=mesh->Elements->Tag[sampleNo];
2117      break;     break;
2118    case(ReducedElements):     case(ReducedElements):
2119      out=mesh->Elements->Tag[sampleNo];     out=mesh->Elements->Tag[sampleNo];
2120      break;     break;
2121    case(FaceElements):     case(FaceElements):
2122      out=mesh->FaceElements->Tag[sampleNo];     out=mesh->FaceElements->Tag[sampleNo];
2123      break;     break;
2124    case(ReducedFaceElements):     case(ReducedFaceElements):
2125      out=mesh->FaceElements->Tag[sampleNo];     out=mesh->FaceElements->Tag[sampleNo];
2126      break;     break;
2127    case(Points):     case(Points):
2128      out=mesh->Points->Tag[sampleNo];     out=mesh->Points->Tag[sampleNo];
2129      break;     break;
2130    case(ContactElementsZero):     case(ContactElementsZero):
2131      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2132      break;     break;
2133    case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
2134      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2135      break;     break;
2136    case(ContactElementsOne):     case(ContactElementsOne):
2137      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2138      break;     break;
2139    case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
2140      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2141      break;     break;
2142    case(DegreesOfFreedom):     case(DegreesOfFreedom):
2143      throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2144      break;     break;
2145    case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2146      throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");     throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2147      break;     break;
2148    default:     default:
2149      stringstream temp;        stringstream temp;
2150      temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2151      throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2152      break;        break;
2153    }     }
2154    return out;     return out;
2155  }  }
2156    
2157    
2158  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
2159  {  {
2160       Finley_Mesh* mesh=m_finleyMesh.get();
2161       escriptDataC tmp=mask.getDataC();
2162       switch(functionSpaceType) {
2163       case(Nodes):
2164       Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2165       break;
2166       case(ReducedNodes):
2167       throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2168       break;
2169       case(DegreesOfFreedom):
2170       throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2171       break;
2172       case(ReducedDegreesOfFreedom):
2173       throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2174       break;
2175       case(Elements):
2176       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2177       break;
2178       case(ReducedElements):
2179       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2180       break;
2181       case(FaceElements):
2182       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2183       break;
2184       case(ReducedFaceElements):
2185       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2186       break;
2187       case(Points):
2188       Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2189       break;
2190       case(ContactElementsZero):
2191       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2192       break;
2193       case(ReducedContactElementsZero):
2194       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2195       break;
2196       case(ContactElementsOne):
2197       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2198       break;
2199       case(ReducedContactElementsOne):
2200       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2201       break;
2202       default:
2203          stringstream temp;
2204          temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2205          throw FinleyAdapterException(temp.str());
2206       }
2207       checkFinleyError();
2208       return;
2209    }
2210    
2211    void MeshAdapter::setTagMap(const std::string& name,  int tag)
2212    {
2213       Finley_Mesh* mesh=m_finleyMesh.get();
2214       Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2215       checkPasoError();
2216       // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2217    }
2218    
2219    int MeshAdapter::getTag(const std::string& name) const
2220    {
2221       Finley_Mesh* mesh=m_finleyMesh.get();
2222       int tag=0;
2223       tag=Finley_Mesh_getTag(mesh, name.c_str());
2224       checkPasoError();
2225       // throwStandardException("MeshAdapter::getTag is not implemented.");
2226       return tag;
2227    }
2228    
2229    bool MeshAdapter::isValidTagName(const std::string& name) const
2230    {
2231       Finley_Mesh* mesh=m_finleyMesh.get();
2232       return Finley_Mesh_isValidTagName(mesh,name.c_str());
2233    }
2234    
2235    std::string MeshAdapter::showTagNames() const
2236    {
2237       stringstream temp;
2238       Finley_Mesh* mesh=m_finleyMesh.get();
2239       Finley_TagMap* tag_map=mesh->TagMap;
2240       while (tag_map) {
2241          temp << tag_map->name;
2242          tag_map=tag_map->next;
2243          if (tag_map) temp << ", ";
2244       }
2245       return temp.str();
2246    }
2247    
2248    int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2249    {
2250    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2251    escriptDataC tmp=mask.getDataC();    dim_t numTags=0;
2252    switch(functionSpaceType) {    switch(functionSpaceCode) {
2253         case(Nodes):     case(Nodes):
2254            Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);            numTags=mesh->Nodes->numTagsInUse;
2255            break;            break;
2256         case(ReducedNodes):     case(ReducedNodes):
2257            throw FinleyAdapterException("Error - ReducedNodes does not support tags");            throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2258            break;            break;
2259         case(DegreesOfFreedom):     case(DegreesOfFreedom):
2260            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2261            break;            break;
2262         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2263            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2264            break;            break;
2265         case(Elements):     case(Elements):
2266            Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     case(ReducedElements):
2267            break;            numTags=mesh->Elements->numTagsInUse;
2268         case(ReducedElements):            break;
2269            Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     case(FaceElements):
2270            break;     case(ReducedFaceElements):
2271         case(FaceElements):            numTags=mesh->FaceElements->numTagsInUse;
2272            Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);            break;
2273       case(Points):
2274              numTags=mesh->Points->numTagsInUse;
2275              break;
2276       case(ContactElementsZero):
2277       case(ReducedContactElementsZero):
2278       case(ContactElementsOne):
2279       case(ReducedContactElementsOne):
2280              numTags=mesh->ContactElements->numTagsInUse;
2281            break;            break;
2282         case(ReducedFaceElements):     default:
2283            Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);        stringstream temp;
2284            break;        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2285         case(Points):        throw FinleyAdapterException(temp.str());
2286            Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);    }
2287      return numTags;
2288    }
2289    
2290    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2291    {
2292      Finley_Mesh* mesh=m_finleyMesh.get();
2293      index_t* tags=NULL;
2294      switch(functionSpaceCode) {
2295       case(Nodes):
2296              tags=mesh->Nodes->tagsInUse;
2297            break;            break;
2298         case(ContactElementsZero):     case(ReducedNodes):
2299            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2300            break;            break;
2301         case(ReducedContactElementsZero):     case(DegreesOfFreedom):
2302            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2303            break;            break;
2304         case(ContactElementsOne):     case(ReducedDegreesOfFreedom):
2305            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2306            break;            break;
2307         case(ReducedContactElementsOne):     case(Elements):
2308            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);     case(ReducedElements):
2309              tags=mesh->Elements->tagsInUse;
2310              break;
2311       case(FaceElements):
2312       case(ReducedFaceElements):
2313              tags=mesh->FaceElements->tagsInUse;
2314              break;
2315       case(Points):
2316              tags=mesh->Points->tagsInUse;
2317              break;
2318       case(ContactElementsZero):
2319       case(ReducedContactElementsZero):
2320       case(ContactElementsOne):
2321       case(ReducedContactElementsOne):
2322              tags=mesh->ContactElements->tagsInUse;
2323            break;            break;
2324         default:     default:
2325            stringstream temp;        stringstream temp;
2326            temp << "Error - Finley does not know anything about function space type " << functionSpaceType;        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2327            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2328    }    }
2329    checkFinleyError();    return tags;
   return;  
2330  }  }
2331    
 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.");  
 }  
2332    
2333  int MeshAdapter::getTag(const std::string& name) const  bool MeshAdapter::canTag(int functionSpaceCode) const
2334  {  {
2335    Finley_Mesh* mesh=m_finleyMesh.get();    switch(functionSpaceCode) {
2336    int tag=0;     case(Nodes):
2337    tag=Finley_Mesh_getTag(mesh, name.c_str());     case(Elements):
2338    checkPasoError();     case(ReducedElements):
2339    // throwStandardException("MeshAdapter::getTag is not implemented.");     case(FaceElements):
2340    return tag;     case(ReducedFaceElements):
2341       case(Points):
2342       case(ContactElementsZero):
2343       case(ReducedContactElementsZero):
2344       case(ContactElementsOne):
2345       case(ReducedContactElementsOne):
2346              return true;
2347       case(ReducedNodes):
2348       case(DegreesOfFreedom):
2349       case(ReducedDegreesOfFreedom):
2350          return false;
2351       default:
2352        return false;
2353      }
2354  }  }
2355    
2356  bool MeshAdapter::isValidTagName(const std::string& name) const  AbstractDomain::StatusType MeshAdapter::getStatus() const
2357  {  {
2358    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2359    return Finley_Mesh_isValidTagName(mesh,name.c_str());    return Finley_Mesh_getStatus(mesh);
2360  }  }
2361    
2362  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();  
 }  
2363    
2364  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1366  
changed lines
  Added in v.2748

  ViewVC Help
Powered by ViewVC 1.1.26