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

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

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

revision 2100 by gross, Wed Nov 26 08:13:00 2008 UTC revision 3379 by gross, Wed Nov 24 04:48:49 2010 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2010 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 18  Line 18 
18  #ifdef USE_NETCDF  #ifdef USE_NETCDF
19  #include <netcdfcpp.h>  #include <netcdfcpp.h>
20  #endif  #endif
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #include "paso/Paso_MPI.h"  
 #endif  
21  extern "C" {  extern "C" {
22  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
23  }  }
 #include <vector>  
24    
25  #define IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH  256  #include <boost/python/import.hpp>
26    
27  using namespace std;  using namespace std;
28  using namespace escript;  using namespace escript;
# Line 88  int MeshAdapter::getMPIRank() const Line 83  int MeshAdapter::getMPIRank() const
83  }  }
84  void MeshAdapter::MPIBarrier() const  void MeshAdapter::MPIBarrier() const
85  {  {
86  #ifdef PASO_MPI  #ifdef ESYS_MPI
87     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
88  #endif  #endif
89     return;     return;
# Line 99  bool MeshAdapter::onMasterProcessor() co Line 94  bool MeshAdapter::onMasterProcessor() co
94  }  }
95    
96    
97    #ifdef ESYS_MPI
98      MPI_Comm
99    #else
100      unsigned int
101    #endif
102    MeshAdapter::getMPIComm() const
103    {
104    #ifdef ESYS_MPI
105        return m_finleyMesh->MPIInfo->comm;
106    #else
107        return 0;
108    #endif
109    }
110    
111    
112  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
113     return m_finleyMesh.get();     return m_finleyMesh.get();
114  }  }
115    
116  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const string& fileName) const
117  {  {
118     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;
119     strcpy(fName,fileName.c_str());     strcpy(fName,fileName.c_str());
# Line 112  void MeshAdapter::write(const std::strin Line 122  void MeshAdapter::write(const std::strin
122     TMPMEMFREE(fName);     TMPMEMFREE(fName);
123  }  }
124    
125  void MeshAdapter::Print_Mesh_Info(const bool full=false) const  void MeshAdapter::Print_Mesh_Info(const bool full) const
126  {  {
127     Finley_PrintMesh_Info(m_finleyMesh.get(), full);     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
128  }  }
129    
130  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const string& fileName) const
131  {  {
132  #ifdef USE_NETCDF  #ifdef USE_NETCDF
133     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
# Line 137  void MeshAdapter::dump(const std::string Line 147  void MeshAdapter::dump(const std::string
147     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
148     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
149     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
150  #ifdef PASO_MPI  #ifdef ESYS_MPI
151     MPI_Status status;     MPI_Status status;
152  #endif  #endif
153    
154  /* Incoming token indicates it's my turn to write */  /* Incoming token indicates it's my turn to write */
155  #ifdef PASO_MPI  #ifdef ESYS_MPI
156     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
157  #endif  #endif
158    
159     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
160                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
161    
162     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
163     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
164     num_Tags = 0;     num_Tags = 0;
165     if (tag_map) {     while (tag_map) {
166        while (tag_map) {        num_Tags++;
167           num_Tags++;        tag_map=tag_map->next;
          tag_map=tag_map->next;  
       }  
168     }     }
169    
170     // NetCDF error handler     // NetCDF error handler
171     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
172     // Create the file.     // Create the file.
173     NcFile dataFile(newFileName, NcFile::Replace);     NcFile dataFile(newFileName, NcFile::Replace);
174       string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
175     // check if writing was successful     // check if writing was successful
176     if (!dataFile.is_valid())     if (!dataFile.is_valid())
177        throw DataException("Error - MeshAdapter::dump: opening of NetCDF file for output failed: " + *newFileName);        throw DataException(msgPrefix+"Open file for output");
178    
179     // Define dimensions (num_Elements and dim_Elements are identical, dim_Elements only appears if > 0)     // Define dimensions (num_Elements and dim_Elements are identical,
180       // dim_Elements only appears if > 0)
181     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
182        throw DataException("Error - MeshAdapter::dump: appending dimension numNodes to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numNodes)");
183     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
184        throw DataException("Error - MeshAdapter::dump: appending dimension numDim to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numDim)");
185     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)) )
186        throw DataException("Error - MeshAdapter::dump: appending dimension mpi_size to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(mpi_size)");
187     if (num_Elements>0)     if (num_Elements>0)
188        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
189           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements)");
190     if (num_FaceElements>0)     if (num_FaceElements>0)
191        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
192           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
193     if (num_ContactElements>0)     if (num_ContactElements>0)
194        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
195           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements)");
196     if (num_Points>0)     if (num_Points>0)
197        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
198           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Points to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Points)");
199     if (num_Elements>0)     if (num_Elements>0)
200        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
201           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
202     if (num_FaceElements>0)     if (num_FaceElements>0)
203        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
204           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
205     if (num_ContactElements>0)     if (num_ContactElements>0)
206        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
207           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements_numNodes)");
208     if (num_Tags>0)     if (num_Tags>0)
209        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
210           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Tags)");
211    
212     // Attributes: MPI size, MPI rank, Name, order, reduced_order     // Attributes: MPI size, MPI rank, Name, order, reduced_order
213     if (!dataFile.add_att("mpi_size", mpi_size) )     if (!dataFile.add_att("mpi_size", mpi_size) )
214        throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_size)");
215     if (!dataFile.add_att("mpi_rank", mpi_rank) )     if (!dataFile.add_att("mpi_rank", mpi_rank) )
216        throw DataException("Error - MeshAdapter::dump: appending mpi_rank to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_rank)");
217     if (!dataFile.add_att("Name",mesh->Name) )     if (!dataFile.add_att("Name",mesh->Name) )
218        throw DataException("Error - MeshAdapter::dump: appending Name to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Name)");
219     if (!dataFile.add_att("numDim",numDim) )     if (!dataFile.add_att("numDim",numDim) )
220        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
221     if (!dataFile.add_att("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
222        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
223     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
224        throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(reduced_order)");
225     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
226        throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(numNodes)");
227     if (!dataFile.add_att("num_Elements",num_Elements) )     if (!dataFile.add_att("num_Elements",num_Elements) )
228        throw DataException("Error - MeshAdapter::dump: appending num_Elements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements)");
229     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
230        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements)");
231     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
232        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements)");
233     if (!dataFile.add_att("num_Points",num_Points) )     if (!dataFile.add_att("num_Points",num_Points) )
234        throw DataException("Error - MeshAdapter::dump: appending num_Points to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Points)");
235     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
236        throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
237     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
238        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
239     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
240        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements_numNodes)");
241     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
242        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Elements_TypeId)");
243     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
244        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
245     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
246        throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(ContactElements_TypeId)");
247     if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
248        throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Points_TypeId)");
249     if (!dataFile.add_att("num_Tags", num_Tags) )     if (!dataFile.add_att("num_Tags", num_Tags) )
250        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Tags)");
251    
252     // // // // // Nodes // // // // //     // // // // // Nodes // // // // //
253    
254     // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)     // Nodes nodeDistribution
255       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
256          throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
257       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
258       if (! (ids->put(int_ptr, mpi_size+1)) )
259          throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
260    
261       // Nodes degreesOfFreedomDistribution
262       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
263          throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
264       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
265       if (! (ids->put(int_ptr, mpi_size+1)) )
266          throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
267    
268       // Only write nodes if non-empty because NetCDF doesn't like empty arrays
269       // (it treats them as NC_UNLIMITED)
270     if (numNodes>0) {     if (numNodes>0) {
271    
272        // Nodes Id        // Nodes Id
273        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
274           throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Id)");
275        int_ptr = &mesh->Nodes->Id[0];        int_ptr = &mesh->Nodes->Id[0];
276        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
277           throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Id)");
278    
279        // Nodes Tag        // Nodes Tag
280        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
281           throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Tag)");
282        int_ptr = &mesh->Nodes->Tag[0];        int_ptr = &mesh->Nodes->Tag[0];
283        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
284           throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Tag)");
285    
286        // Nodes gDOF        // Nodes gDOF
287        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
288           throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
289        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
290        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
291           throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gDOF)");
292    
293        // Nodes global node index        // Nodes global node index
294        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
295           throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gNI)");
296        int_ptr = &mesh->Nodes->globalNodesIndex[0];        int_ptr = &mesh->Nodes->globalNodesIndex[0];
297        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
298           throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gNI)");
299    
300        // Nodes grDof        // Nodes grDof
301        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
302           throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
303        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
304        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
305           throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grDfI)");
306    
307        // Nodes grNI        // Nodes grNI
308        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
309           throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grNI)");
310        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
311        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
312           throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grNI)");
313    
314        // Nodes Coordinates        // Nodes Coordinates
315        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
316           throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
317        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
318           throw DataException("Error - MeshAdapter::dump: copy Nodes_Coordinates to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Coordinates)");
   
       // Nodes degreesOfFreedomDistribution  
       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )  
          throw DataException("Error - MeshAdapter::dump: appending Nodes_DofDistribution to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];  
       if (! (ids->put(int_ptr, mpi_size+1)) )  
          throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);  
   
       // Nodes nodeDistribution  
       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )  
          throw DataException("Error - MeshAdapter::dump: appending Nodes_NodeDistribution to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];  
       if (! (ids->put(int_ptr, mpi_size+1)) )  
          throw DataException("Error - MeshAdapter::dump: copy Nodes_NodeDistribution to netCDF buffer failed: " + *newFileName);  
319    
320     }     }
321    
# Line 314  void MeshAdapter::dump(const std::string Line 325  void MeshAdapter::dump(const std::string
325    
326        // Elements_Id        // Elements_Id
327        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
328           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Id)");
329        int_ptr = &mesh->Elements->Id[0];        int_ptr = &mesh->Elements->Id[0];
330        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
331           throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Id)");
332    
333        // Elements_Tag        // Elements_Tag
334        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
335           throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Tag)");
336        int_ptr = &mesh->Elements->Tag[0];        int_ptr = &mesh->Elements->Tag[0];
337        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
338           throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Tag)");
339    
340        // Elements_Owner        // Elements_Owner
341        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
342           throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Owner)");
343        int_ptr = &mesh->Elements->Owner[0];        int_ptr = &mesh->Elements->Owner[0];
344        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
345           throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Owner)");
346    
347        // Elements_Color        // Elements_Color
348        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
349           throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Color)");
350        int_ptr = &mesh->Elements->Color[0];        int_ptr = &mesh->Elements->Color[0];
351        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
352           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Color)");
353    
354        // Elements_Nodes        // Elements_Nodes
355        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
356           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Nodes)");
357        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
358           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Nodes)");
359    
360     }     }
361    
# Line 354  void MeshAdapter::dump(const std::string Line 365  void MeshAdapter::dump(const std::string
365    
366        // FaceElements_Id        // FaceElements_Id
367        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
368           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Id)");
369        int_ptr = &mesh->FaceElements->Id[0];        int_ptr = &mesh->FaceElements->Id[0];
370        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
371           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Id)");
372    
373        // FaceElements_Tag        // FaceElements_Tag
374        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
375           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
376        int_ptr = &mesh->FaceElements->Tag[0];        int_ptr = &mesh->FaceElements->Tag[0];
377        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
378           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Tag)");
379    
380        // FaceElements_Owner        // FaceElements_Owner
381        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
382           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
383        int_ptr = &mesh->FaceElements->Owner[0];        int_ptr = &mesh->FaceElements->Owner[0];
384        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
385           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Owner)");
386    
387        // FaceElements_Color        // FaceElements_Color
388        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
389           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Color)");
390        int_ptr = &mesh->FaceElements->Color[0];        int_ptr = &mesh->FaceElements->Color[0];
391        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
392           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Color)");
393    
394        // FaceElements_Nodes        // FaceElements_Nodes
395        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
396           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
397        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
398           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Nodes)");
399    
400     }     }
401    
# Line 394  void MeshAdapter::dump(const std::string Line 405  void MeshAdapter::dump(const std::string
405    
406        // ContactElements_Id        // ContactElements_Id
407        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
408           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Id)");
409        int_ptr = &mesh->ContactElements->Id[0];        int_ptr = &mesh->ContactElements->Id[0];
410        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
411           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Id)");
412    
413        // ContactElements_Tag        // ContactElements_Tag
414        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
415           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Tag)");
416        int_ptr = &mesh->ContactElements->Tag[0];        int_ptr = &mesh->ContactElements->Tag[0];
417        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
418           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Tag)");
419    
420        // ContactElements_Owner        // ContactElements_Owner
421        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
422           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Owner)");
423        int_ptr = &mesh->ContactElements->Owner[0];        int_ptr = &mesh->ContactElements->Owner[0];
424        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
425           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Owner)");
426    
427        // ContactElements_Color        // ContactElements_Color
428        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
429           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Color)");
430        int_ptr = &mesh->ContactElements->Color[0];        int_ptr = &mesh->ContactElements->Color[0];
431        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
432           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Color)");
433    
434        // ContactElements_Nodes        // ContactElements_Nodes
435        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
436           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Nodes)");
437        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
438           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Nodes)");
439    
440     }     }
441    
# Line 436  void MeshAdapter::dump(const std::string Line 447  void MeshAdapter::dump(const std::string
447    
448        // Points_Id        // Points_Id
449        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
450           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Id)");
451        int_ptr = &mesh->Points->Id[0];        int_ptr = &mesh->Points->Id[0];
452        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
453           throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Id)");
454    
455        // Points_Tag        // Points_Tag
456        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
457           throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Tag)");
458        int_ptr = &mesh->Points->Tag[0];        int_ptr = &mesh->Points->Tag[0];
459        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
460           throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Tag)");
461    
462        // Points_Owner        // Points_Owner
463        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
464           throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Owner)");
465        int_ptr = &mesh->Points->Owner[0];        int_ptr = &mesh->Points->Owner[0];
466        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
467           throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Owner)");
468    
469        // Points_Color        // Points_Color
470        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
471           throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Color)");
472        int_ptr = &mesh->Points->Color[0];        int_ptr = &mesh->Points->Color[0];
473        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
474           throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Color)");
475    
476        // Points_Nodes        // Points_Nodes
477        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
478        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
479           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Nodes)");
480        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
481           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Nodes)");
482    
483     }     }
484    
# Line 491  void MeshAdapter::dump(const std::string Line 502  void MeshAdapter::dump(const std::string
502    
503        // Tags_keys        // Tags_keys
504        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
505           throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Tags_keys)");
506        int_ptr = &Tags_keys[0];        int_ptr = &Tags_keys[0];
507        if (! (ids->put(int_ptr, num_Tags)) )        if (! (ids->put(int_ptr, num_Tags)) )
508           throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Tags_keys)");
509    
510        // Tags_names_*        // Tags_names_*
511        // 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
512        // because the NetCDF manual doesn't tell how to do an array of strings        // instead I have hacked in one attribute per string because the NetCDF
513          // 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(msgPrefix+"add_att(Tags_names_XX)");
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  /* Send token to next MPI process so he can take his turn */  /* Send token to next MPI process so he can take his turn */
530  #ifdef PASO_MPI  #ifdef ESYS_MPI
531     if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);     if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
532  #endif  #endif
533    
# Line 681  pair<int,int> MeshAdapter::getDataShape( Line 692  pair<int,int> MeshAdapter::getDataShape(
692     case(Elements):     case(Elements):
693     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
694        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
695        numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
696     }     }
697     break;     break;
698     case(ReducedElements):     case(ReducedElements):
699     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
700        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
701        numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
702     }     }
703     break;     break;
704     case(FaceElements):     case(FaceElements):
705     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
706        numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
707        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
708     }     }
709     break;     break;
710     case(ReducedFaceElements):     case(ReducedFaceElements):
711     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
712        numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
713        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
714     }     }
715     break;     break;
# Line 710  pair<int,int> MeshAdapter::getDataShape( Line 721  pair<int,int> MeshAdapter::getDataShape(
721     break;     break;
722     case(ContactElementsZero):     case(ContactElementsZero):
723     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
724        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
725        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
726     }     }
727     break;     break;
728     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
729     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
730        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
731        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
732     }     }
733     break;     break;
734     case(ContactElementsOne):     case(ContactElementsOne):
735     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
736        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
737        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
738     }     }
739     break;     break;
740     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
741     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
742        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
743        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
744     }     }
745     break;     break;
# Line 757  pair<int,int> MeshAdapter::getDataShape( Line 768  pair<int,int> MeshAdapter::getDataShape(
768  // 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:
769  //  //
770  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
771                                   SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
772                                   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,
773                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
774                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact) const
775  {  {
776       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
777       if (smat==0)
778       {
779        throw FinleyAdapterException("finley only supports its own system matrices.");
780       }
781     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
782     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
783     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 776  void MeshAdapter::addPDEToSystem( Line 792  void MeshAdapter::addPDEToSystem(
792    
793     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
794    
795     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
796     checkFinleyError();     checkFinleyError();
797    
798     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
799     checkFinleyError();     checkFinleyError();
800    
801     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
802     checkFinleyError();     checkFinleyError();
803  }  }
804    
805  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
806                                          escript::Data& mat,                                          escript::Data& mat,
807                                          const escript::Data& D,                                          const escript::Data& D,
808                                          const escript::Data& d) const                                          const escript::Data& d,
809                                            const bool useHRZ) const
810  {  {
811     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
812     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
# Line 797  void  MeshAdapter::addPDEToLumpedSystem( Line 814  void  MeshAdapter::addPDEToLumpedSystem(
814    
815     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
816    
817     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
818     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     checkFinleyError();
819      
820       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
821     checkFinleyError();     checkFinleyError();
822    
823  }  }
824    
825    
# Line 830  void MeshAdapter::addPDEToRHS( escript:: Line 849  void MeshAdapter::addPDEToRHS( escript::
849  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
850  //  //
851  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
852                                             TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
853                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
854                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
855                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
856                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact) const
857  {  {
858       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
859       if (tpa==0)
860       {
861        throw FinleyAdapterException("finley only supports its own transport problems.");
862       }
863    
864    
865     DataTypes::ShapeType shape;     DataTypes::ShapeType shape;
866     source.expand();     source.expand();
867     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
# Line 852  void MeshAdapter::addPDEToTransportProbl Line 878  void MeshAdapter::addPDEToTransportProbl
878     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
879    
880     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
881     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
882    
883     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
884     checkFinleyError();     checkFinleyError();
# Line 884  void MeshAdapter::interpolateOnDomain(es Line 910  void MeshAdapter::interpolateOnDomain(es
910     escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
911     switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
912     case(Nodes):     case(Nodes):
913     switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
914     case(Nodes):        case(Nodes):
915     case(ReducedNodes):        case(ReducedNodes):
916     case(DegreesOfFreedom):        case(DegreesOfFreedom):
917     case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
918     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
    break;  
    case(Elements):  
    case(ReducedElements):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);  
    break;  
    case(FaceElements):  
    case(ReducedFaceElements):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);  
    break;  
    case(Points):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);  
    break;  
    case(ContactElementsZero):  
    case(ReducedContactElementsZero):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
919        break;        break;
920     }        case(Elements):
921     break;        case(ReducedElements):
922     case(ReducedNodes):        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
923     switch(target.getFunctionSpace().getTypeCode()) {        break;
924     case(Nodes):        case(FaceElements):
925     case(ReducedNodes):        case(ReducedFaceElements):
926     case(DegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
927     case(ReducedDegreesOfFreedom):        break;
928     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(Points):
929     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
930     case(Elements):        break;
931     case(ReducedElements):        case(ContactElementsZero):
932     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        case(ReducedContactElementsZero):
933     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
934     case(FaceElements):        break;
935     case(ReducedFaceElements):        case(ContactElementsOne):
936     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        case(ReducedContactElementsOne):
937     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
938     case(Points):        break;
939     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        default:
940     break;           stringstream temp;
941     case(ContactElementsZero):           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
942     case(ReducedContactElementsZero):           throw FinleyAdapterException(temp.str());
943     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           break;
944     break;        }
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
945        break;        break;
    }  
    break;  
    case(Elements):  
    if (target.getFunctionSpace().getTypeCode()==Elements) {  
       Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
       Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on elements possible.");  
    }  
    break;  
    case(ReducedElements):  
    if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
       Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");  
    }  
    break;  
    case(FaceElements):  
    if (target.getFunctionSpace().getTypeCode()==FaceElements) {  
       Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {  
       Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");  
    }  
    break;  
    case(ReducedFaceElements):  
    if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {  
       Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");  
    }  
    break;  
    case(Points):  
    if (target.getFunctionSpace().getTypeCode()==Points) {  
       Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on points possible.");  
    }  
    break;  
    case(ContactElementsZero):  
    case(ContactElementsOne):  
    if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {  
       Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {  
       Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");  
    }  
    break;  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {  
       Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");  
    }  
    break;  
    case(DegreesOfFreedom):        
    switch(target.getFunctionSpace().getTypeCode()) {  
    case(ReducedDegreesOfFreedom):  
    case(DegreesOfFreedom):  
    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    break;  
   
    case(Nodes):  
946     case(ReducedNodes):     case(ReducedNodes):
947     if (getMPISize()>1) {        switch(target.getFunctionSpace().getTypeCode()) {
948        escript::Data temp=escript::Data(in);        case(Nodes):
949        temp.expand();        case(ReducedNodes):
950        escriptDataC _in2 = temp.getDataC();        case(DegreesOfFreedom):
951        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);        case(ReducedDegreesOfFreedom):
    } else {  
952        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
953     }        break;
954     break;        case(Elements):
955     case(Elements):        case(ReducedElements):
    case(ReducedElements):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);  
    } else {  
956        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
957     }        break;
958     break;        case(FaceElements):
959     case(FaceElements):        case(ReducedFaceElements):
    case(ReducedFaceElements):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);  
   
    } else {  
960        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
961     }        break;
962     break;        case(Points):
    case(Points):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
    } else {  
963        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
964     }        break;
965     break;        case(ContactElementsZero):
966     case(ContactElementsZero):        case(ReducedContactElementsZero):
    case(ContactElementsOne):  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);  
    } else {  
967        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
    }  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
968        break;        break;
969     }        case(ContactElementsOne):
970     break;        case(ReducedContactElementsOne):
971     case(ReducedDegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
972     switch(target.getFunctionSpace().getTypeCode()) {        break;
973     case(Nodes):        default:
974     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           stringstream temp;
975     break;           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
976     case(ReducedNodes):           throw FinleyAdapterException(temp.str());
977     if (getMPISize()>1) {           break;
978        escript::Data temp=escript::Data(in);        }
979        temp.expand();        break;
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);  
    } else {  
       Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    }  
    break;  
    case(DegreesOfFreedom):  
    throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");  
    break;  
    case(ReducedDegreesOfFreedom):  
    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    break;  
980     case(Elements):     case(Elements):
981          if (target.getFunctionSpace().getTypeCode()==Elements) {
982             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
983          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
984             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
985          } else {
986             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
987          }
988          break;
989     case(ReducedElements):     case(ReducedElements):
990     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
991        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
992        escriptDataC _in2 = temp.getDataC();        } else {
993        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
994     } else {        }
995        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        break;
    }  
    break;  
996     case(FaceElements):     case(FaceElements):
997          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
998             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
999          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1000             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1001          } else {
1002             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1003          }
1004          break;
1005     case(ReducedFaceElements):     case(ReducedFaceElements):
1006     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1007        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1008        escriptDataC _in2 = temp.getDataC();        } else {
1009        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
1010     } else {        }
1011        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
    }  
    break;  
1012     case(Points):     case(Points):
1013     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==Points) {
1014        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1015        escriptDataC _in2 = temp.getDataC();        } else {
1016        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1017     } else {        }
1018        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        break;
    }  
    break;  
1019     case(ContactElementsZero):     case(ContactElementsZero):
1020     case(ContactElementsOne):     case(ContactElementsOne):
1021          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1022             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1023          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1024             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1025          } else {
1026             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1027          }
1028          break;
1029     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1030     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1031     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1032        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1033        escriptDataC _in2 = temp.getDataC();        } else {
1034        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1035     } else {        }
1036        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        break;
1037     }     case(DegreesOfFreedom):      
1038     break;        switch(target.getFunctionSpace().getTypeCode()) {
1039     default:        case(ReducedDegreesOfFreedom):
1040        stringstream temp;        case(DegreesOfFreedom):
1041        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1042        throw FinleyAdapterException(temp.str());        break;
1043      
1044          case(Nodes):
1045          case(ReducedNodes):
1046          if (getMPISize()>1) {
1047             escript::Data temp=escript::Data(in);
1048             temp.expand();
1049             escriptDataC _in2 = temp.getDataC();
1050             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1051          } else {
1052             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1053          }
1054          break;
1055          case(Elements):
1056          case(ReducedElements):
1057          if (getMPISize()>1) {
1058             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1059             escriptDataC _in2 = temp.getDataC();
1060             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1061          } else {
1062             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1063          }
1064          break;
1065          case(FaceElements):
1066          case(ReducedFaceElements):
1067          if (getMPISize()>1) {
1068             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1069             escriptDataC _in2 = temp.getDataC();
1070             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1071      
1072          } else {
1073             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1074          }
1075          break;
1076          case(Points):
1077          if (getMPISize()>1) {
1078             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1079             escriptDataC _in2 = temp.getDataC();
1080          } else {
1081             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1082          }
1083          break;
1084          case(ContactElementsZero):
1085          case(ContactElementsOne):
1086          case(ReducedContactElementsZero):
1087          case(ReducedContactElementsOne):
1088          if (getMPISize()>1) {
1089             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1090             escriptDataC _in2 = temp.getDataC();
1091             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1092          } else {
1093             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1094          }
1095          break;
1096          default:
1097             stringstream temp;
1098             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1099             throw FinleyAdapterException(temp.str());
1100             break;
1101          }
1102          break;
1103       case(ReducedDegreesOfFreedom):
1104          switch(target.getFunctionSpace().getTypeCode()) {
1105          case(Nodes):
1106          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1107          break;
1108          case(ReducedNodes):
1109          if (getMPISize()>1) {
1110             escript::Data temp=escript::Data(in);
1111             temp.expand();
1112             escriptDataC _in2 = temp.getDataC();
1113             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1114          } else {
1115             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1116          }
1117          break;
1118          case(DegreesOfFreedom):
1119          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1120          break;
1121          case(ReducedDegreesOfFreedom):
1122          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1123          break;
1124          case(Elements):
1125          case(ReducedElements):
1126          if (getMPISize()>1) {
1127             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1128             escriptDataC _in2 = temp.getDataC();
1129             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1130          } else {
1131             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1132          }
1133          break;
1134          case(FaceElements):
1135          case(ReducedFaceElements):
1136          if (getMPISize()>1) {
1137             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1138             escriptDataC _in2 = temp.getDataC();
1139             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1140          } else {
1141             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1142          }
1143          break;
1144          case(Points):
1145          if (getMPISize()>1) {
1146             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1147             escriptDataC _in2 = temp.getDataC();
1148             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1149          } else {
1150             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1151          }
1152          break;
1153          case(ContactElementsZero):
1154          case(ContactElementsOne):
1155          case(ReducedContactElementsZero):
1156          case(ReducedContactElementsOne):
1157          if (getMPISize()>1) {
1158             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1159             escriptDataC _in2 = temp.getDataC();
1160             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1161          } else {
1162             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1163          }
1164          break;
1165          default:
1166             stringstream temp;
1167             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1168             throw FinleyAdapterException(temp.str());
1169             break;
1170          }
1171        break;        break;
    }  
    break;  
1172     default:     default:
1173        stringstream temp;        stringstream temp;
1174        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
# Line 1166  void MeshAdapter::setToX(escript::Data& Line 1192  void MeshAdapter::setToX(escript::Data&
1192        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1193        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1194     } else {     } else {
1195        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1196        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1197        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1198        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1247  void MeshAdapter::interpolateACross(escr Line 1273  void MeshAdapter::interpolateACross(escr
1273  //  //
1274  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
1275  //  //
1276  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1277  {  {
1278     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1279     if (argDomain!=*this)     if (argDomain!=*this)
# Line 1260  void MeshAdapter::setToIntegrals(std::ve Line 1286  void MeshAdapter::setToIntegrals(std::ve
1286     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1287     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1288     case(Nodes):     case(Nodes):
1289     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1290     _temp=temp.getDataC();     _temp=temp.getDataC();
1291     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1292     break;     break;
1293     case(ReducedNodes):     case(ReducedNodes):
1294     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1295     _temp=temp.getDataC();     _temp=temp.getDataC();
1296     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1297     break;     break;
# Line 1297  void MeshAdapter::setToIntegrals(std::ve Line 1323  void MeshAdapter::setToIntegrals(std::ve
1323     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1324     break;     break;
1325     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1326     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1327     _temp=temp.getDataC();     _temp=temp.getDataC();
1328     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1329     break;     break;
1330     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1331     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1332     _temp=temp.getDataC();     _temp=temp.getDataC();
1333     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1334     break;     break;
# Line 1334  void MeshAdapter::setToGradient(escript: Line 1360  void MeshAdapter::setToGradient(escript:
1360     escript::Data temp;     escript::Data temp;
1361     if (getMPISize()>1) {     if (getMPISize()>1) {
1362        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1363           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1364           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1365        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1366           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1367           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1368        } else {        } else {
1369           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1446  void MeshAdapter::setToSize(escript::Dat Line 1472  void MeshAdapter::setToSize(escript::Dat
1472     checkFinleyError();     checkFinleyError();
1473  }  }
1474    
1475  // sets the location of nodes:  //
1476    // sets the location of nodes
1477    //
1478  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1479  {  {
1480     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1454  void MeshAdapter::setNewX(const escript: Line 1482  void MeshAdapter::setNewX(const escript:
1482     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1483     if (newDomain!=*this)     if (newDomain!=*this)
1484        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1485     tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1486     Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1487           Finley_Mesh_setCoordinates(mesh,&tmp);
1488       } else {
1489           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1490           tmp = new_x_inter.getDataC();
1491           Finley_Mesh_setCoordinates(mesh,&tmp);
1492       }
1493     checkFinleyError();     checkFinleyError();
1494  }  }
1495    
1496  // saves a data array in openDX format:  //
1497  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  // Helper for the save* methods. Extracts optional data variable names and the
1498    // corresponding pointers from python dictionary. Caller must free arrays.
1499    //
1500    void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1501  {  {
1502     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;     numData = boost::python::extract<int>(arg.attr("__len__")());
    const int num_data=boost::python::extract<int>(arg.attr("__len__")());  
1503     /* win32 refactor */     /* win32 refactor */
1504     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1505     for(int i=0;i<num_data;i++)     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1506     {     dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
       names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
    }  
   
    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  
    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  
1507    
1508     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1509     for (int i=0;i<num_data;++i) {     for (int i=0; i<numData; ++i) {
1510        std::string n=boost::python::extract<std::string>(keys[i]);        string n=boost::python::extract<string>(keys[i]);
1511        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1512        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1513           throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");           throw FinleyAdapterException("Error: Data must be defined on same Domain");
1514        data[i]=d.getDataC();        data[i] = d.getDataC();
1515        ptr_data[i]= &(data[i]);        dataPtr[i] = &(data[i]);
1516        if (n.length()>MAX_namelength-1) {        names[i] = TMPMEMALLOC(n.length()+1, char);
1517           strncpy(names[i],n.c_str(),MAX_namelength-1);        strcpy(names[i], n.c_str());
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
1518     }     }
1519     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  }
1520    
1521    //
1522    // saves mesh and optionally data arrays in openDX format
1523    //
1524    void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1525    {
1526       int num_data;
1527       char **names;
1528       escriptDataC *data;
1529       escriptDataC **ptr_data;
1530    
1531       extractArgsFromDict(arg, num_data, names, data, ptr_data);
1532       Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1533     checkFinleyError();     checkFinleyError();
1534    
1535     /* win32 refactor */     /* win32 refactor */
1536     TMPMEMFREE(data);     TMPMEMFREE(data);
1537     TMPMEMFREE(ptr_data);     TMPMEMFREE(ptr_data);
1538     for(int i=0;i<num_data;i++)     for(int i=0; i<num_data; i++)
1539     {     {
1540        TMPMEMFREE(names[i]);        TMPMEMFREE(names[i]);
1541     }     }
# Line 1504  void MeshAdapter::saveDX(const std::stri Line 1544  void MeshAdapter::saveDX(const std::stri
1544     return;     return;
1545  }  }
1546    
1547  // saves a data array in openVTK format:  //
1548  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  // saves mesh and optionally data arrays in VTK format
1549    //
1550    void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1551  {  {
1552     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1553     const int num_data=boost::python::extract<int>(arg.attr("__len__")());      pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1554     /* win32 refactor */                metadata, metadata_schema, arg);
1555     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;  }
1556     for(int i=0;i<num_data;i++)  
1557     {  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1558        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  {
1559     }  #ifdef ESYS_MPI
1560        index_t myFirstNode=0, myLastNode=0, k=0;
1561        index_t* globalNodeIndex=0;
1562        Finley_Mesh* mesh_p=m_finleyMesh.get();
1563        if (fs_code == FINLEY_REDUCED_NODES)
1564        {
1565        myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1566        myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1567        globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1568        }
1569        else
1570        {
1571        myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1572        myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1573        globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1574        }
1575        k=globalNodeIndex[id];
1576        return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1577    #endif
1578        return true;
1579    }
1580    
    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  
    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  
1581    
    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]);  
       if (n.length()>MAX_namelength-1) {  
          strncpy(names[i],n.c_str(),MAX_namelength-1);  
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
    }  
    Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  
1582    
1583     checkFinleyError();  //
1584     /* win32 refactor */  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1585     TMPMEMFREE(data);  //
1586     TMPMEMFREE(ptr_data);  ASM_ptr MeshAdapter::newSystemMatrix(
    for(int i=0;i<num_data;i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
   
    return;  
 }  
                                                                                                                                                                       
                                                                                                                                                                       
 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  
 SystemMatrixAdapter MeshAdapter::newSystemMatrix(  
1587                                                   const int row_blocksize,                                                   const int row_blocksize,
1588                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1589                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1594  SystemMatrixAdapter MeshAdapter::newSyst Line 1626  SystemMatrixAdapter MeshAdapter::newSyst
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);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1634       return ASM_ptr(sma);
1635    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1636  }  }
1637    
1638    //
1639  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1640  TransportProblemAdapter MeshAdapter::newTransportProblem(  //
1641                                                           const double theta,  ATP_ptr MeshAdapter::newTransportProblem(
1642                                                             const bool useBackwardEuler,
1643                                                           const int blocksize,                                                           const int blocksize,
1644                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1645                                                           const int type) const                                                           const int type) const
# Line 1624  TransportProblemAdapter MeshAdapter::new Line 1661  TransportProblemAdapter MeshAdapter::new
1661    
1662     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1663     checkFinleyError();     checkFinleyError();
1664     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1665     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1666     checkPasoError();     checkPasoError();
1667     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1668     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1669       return ATP_ptr(tpa);
1670    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1671  }  }
1672    
1673  //  //
# Line 1667  bool MeshAdapter::isCellOriented(int fun Line 1706  bool MeshAdapter::isCellOriented(int fun
1706     return false;     return false;
1707  }  }
1708    
1709    bool
1710    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1711    {
1712       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1713        class 1: DOF <-> Nodes
1714        class 2: ReducedDOF <-> ReducedNodes
1715        class 3: Points
1716        class 4: Elements
1717        class 5: ReducedElements
1718        class 6: FaceElements
1719        class 7: ReducedFaceElements
1720        class 8: ContactElementZero <-> ContactElementOne
1721        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1722    
1723       There is also a set of lines. Interpolation is possible down a line but not between lines.
1724       class 1 and 2 belong to all lines so aren't considered.
1725        line 0: class 3
1726        line 1: class 4,5
1727        line 2: class 6,7
1728        line 3: class 8,9
1729    
1730       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1731       eg hasnodes is true if we have at least one instance of Nodes.
1732       */
1733        if (fs.empty())
1734        {
1735            return false;
1736        }
1737        vector<int> hasclass(10);
1738        vector<int> hasline(4);
1739        bool hasnodes=false;
1740        bool hasrednodes=false;
1741        bool hascez=false;
1742        bool hasrcez=false;
1743        for (int i=0;i<fs.size();++i)
1744        {
1745        switch(fs[i])
1746        {
1747        case(Nodes):   hasnodes=true;   // no break is deliberate
1748        case(DegreesOfFreedom):
1749            hasclass[1]=1;
1750            break;
1751        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1752        case(ReducedDegreesOfFreedom):
1753            hasclass[2]=1;
1754            break;
1755        case(Points):
1756            hasline[0]=1;
1757            hasclass[3]=1;
1758            break;
1759        case(Elements):
1760            hasclass[4]=1;
1761            hasline[1]=1;
1762            break;
1763        case(ReducedElements):
1764            hasclass[5]=1;
1765            hasline[1]=1;
1766            break;
1767        case(FaceElements):
1768            hasclass[6]=1;
1769            hasline[2]=1;
1770            break;
1771        case(ReducedFaceElements):
1772            hasclass[7]=1;
1773            hasline[2]=1;
1774            break;
1775        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1776        case(ContactElementsOne):
1777            hasclass[8]=1;
1778            hasline[3]=1;
1779            break;
1780        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1781        case(ReducedContactElementsOne):
1782            hasclass[9]=1;
1783            hasline[3]=1;
1784            break;
1785        default:
1786            return false;
1787        }
1788        }
1789        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1790        // fail if we have more than one leaf group
1791    
1792        if (totlines>1)
1793        {
1794        return false;   // there are at least two branches we can't interpolate between
1795        }
1796        else if (totlines==1)
1797        {
1798        if (hasline[0]==1)      // we have points
1799        {
1800            resultcode=Points;
1801        }
1802        else if (hasline[1]==1)
1803        {
1804            if (hasclass[5]==1)
1805            {
1806            resultcode=ReducedElements;
1807            }
1808            else
1809            {
1810            resultcode=Elements;
1811            }
1812        }
1813        else if (hasline[2]==1)
1814        {
1815            if (hasclass[7]==1)
1816            {
1817            resultcode=ReducedFaceElements;
1818            }
1819            else
1820            {
1821            resultcode=FaceElements;
1822            }
1823        }
1824        else    // so we must be in line3
1825        {
1826            if (hasclass[9]==1)
1827            {
1828            // need something from class 9
1829            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1830            }
1831            else
1832            {
1833            // something from class 8
1834            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1835            }
1836        }
1837        }
1838        else    // totlines==0
1839        {
1840        if (hasclass[2]==1)
1841        {
1842            // something from class 2
1843            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1844        }
1845        else
1846        {   // something from class 1
1847            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1848        }
1849        }
1850        return true;
1851    }
1852    
1853  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1854  {  {
1855     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1856     case(Nodes):     case(Nodes):
1857     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1858     case(Nodes):      case(Nodes):
1859     case(ReducedNodes):      case(ReducedNodes):
1860     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1861     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1862     case(Elements):      case(Elements):
1863     case(ReducedElements):      case(ReducedElements):
1864     case(FaceElements):      case(FaceElements):
1865     case(ReducedFaceElements):      case(ReducedFaceElements):
1866     case(Points):      case(Points):
1867     case(ContactElementsZero):      case(ContactElementsZero):
1868     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1869     case(ContactElementsOne):      case(ContactElementsOne):
1870     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1871     return true;      return true;
1872     default:      default:
1873        stringstream temp;            stringstream temp;
1874        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;
1875        throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1876     }     }
1877     break;     break;
1878     case(ReducedNodes):     case(ReducedNodes):
1879     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1880     case(ReducedNodes):      case(ReducedNodes):
1881     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1882     case(Elements):      case(Elements):
1883     case(ReducedElements):      case(ReducedElements):
1884     case(FaceElements):      case(FaceElements):
1885     case(ReducedFaceElements):      case(ReducedFaceElements):
1886     case(Points):      case(Points):
1887     case(ContactElementsZero):      case(ContactElementsZero):
1888     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1889     case(ContactElementsOne):      case(ContactElementsOne):
1890     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1891     return true;      return true;
1892     case(Nodes):      case(Nodes):
1893     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1894     return false;      return false;
1895     default:      default:
1896        stringstream temp;          stringstream temp;
1897        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;
1898        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1899     }     }
1900     break;     break;
1901     case(Elements):     case(Elements):
1902     if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1903        return true;        return true;
1904     } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1905        return true;        return true;
1906     } else {          } else {
1907        return false;            return false;
1908     }          }
1909     case(ReducedElements):     case(ReducedElements):
1910     if (functionSpaceType_target==ReducedElements) {      if (functionSpaceType_target==ReducedElements) {
1911        return true;        return true;
1912     } else {      } else {
1913        return false;            return false;
1914     }      }
1915     case(FaceElements):     case(FaceElements):
1916     if (functionSpaceType_target==FaceElements) {      if (functionSpaceType_target==FaceElements) {
1917        return true;              return true;
1918     } else if (functionSpaceType_target==ReducedFaceElements) {      } else if (functionSpaceType_target==ReducedFaceElements) {
1919        return true;              return true;
1920     } else {      } else {
1921        return false;              return false;
1922     }      }
1923     case(ReducedFaceElements):     case(ReducedFaceElements):
1924     if (functionSpaceType_target==ReducedFaceElements) {      if (functionSpaceType_target==ReducedFaceElements) {
1925        return true;              return true;
1926     } else {      } else {
1927        return false;          return false;
1928     }      }
1929     case(Points):     case(Points):
1930     if (functionSpaceType_target==Points) {      if (functionSpaceType_target==Points) {
1931        return true;              return true;
1932     } else {      } else {
1933        return false;              return false;
1934     }      }
1935     case(ContactElementsZero):     case(ContactElementsZero):
1936     case(ContactElementsOne):     case(ContactElementsOne):
1937     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1938        return true;              return true;
1939     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1940        return true;              return true;
1941     } else {      } else {
1942        return false;              return false;
1943     }      }
1944     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1945     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1946     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1947        return true;              return true;
1948     } else {      } else {
1949        return false;              return false;
1950     }      }
1951     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1952     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1953     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1954     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1955     case(Nodes):      case(Nodes):
1956     case(ReducedNodes):      case(ReducedNodes):
1957     case(Elements):      case(Elements):
1958     case(ReducedElements):      case(ReducedElements):
1959     case(Points):      case(Points):
1960     case(FaceElements):      case(FaceElements):
1961     case(ReducedFaceElements):      case(ReducedFaceElements):
1962     case(ContactElementsZero):      case(ContactElementsZero):
1963     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1964     case(ContactElementsOne):      case(ContactElementsOne):
1965     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1966     return true;      return true;
1967     default:      default:
1968        stringstream temp;          stringstream temp;
1969        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;
1970        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1971     }      }
1972     break;      break;
1973     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1974     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
1975     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1976     case(ReducedNodes):      case(ReducedNodes):
1977     case(Elements):      case(Elements):
1978     case(ReducedElements):      case(ReducedElements):
1979     case(FaceElements):      case(FaceElements):
1980     case(ReducedFaceElements):      case(ReducedFaceElements):
1981     case(Points):      case(Points):
1982     case(ContactElementsZero):      case(ContactElementsZero):
1983     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1984     case(ContactElementsOne):      case(ContactElementsOne):
1985     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1986     return true;      return true;
1987     case(Nodes):      case(Nodes):
1988     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1989     return false;      return false;
1990     default:      default:
1991        stringstream temp;          stringstream temp;
1992        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;
1993        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1994     }      }
1995     break;      break;
1996     default:     default:
1997        stringstream temp;        stringstream temp;
1998        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
# Line 1842  bool MeshAdapter::operator!=(const Abstr Line 2025  bool MeshAdapter::operator!=(const Abstr
2025    
2026  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2027  {  {
2028     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2029       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2030     checkPasoError();     checkPasoError();
2031     return out;     return out;
2032  }  }
2033  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2034  {  {
2035     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2036       int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2037     checkPasoError();     checkPasoError();
2038     return out;     return out;
2039  }  }
2040    
2041  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2042  {  {
2043     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2044  }  }
2045    
2046  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2047  {  {
2048     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2049  }  }
2050    
2051  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2052  {  {
2053     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2054  }  }
2055    
2056  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2057  {  {
2058     int *out = NULL;     int *out = NULL;
2059     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2027  void MeshAdapter::setTags(const int func Line 2212  void MeshAdapter::setTags(const int func
2212     return;     return;
2213  }  }
2214    
2215  void MeshAdapter::setTagMap(const std::string& name,  int tag)  void MeshAdapter::setTagMap(const string& name,  int tag)
2216  {  {
2217     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2218     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
# Line 2035  void MeshAdapter::setTagMap(const std::s Line 2220  void MeshAdapter::setTagMap(const std::s
2220     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2221  }  }
2222    
2223  int MeshAdapter::getTag(const std::string& name) const  int MeshAdapter::getTag(const string& name) const
2224  {  {
2225     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2226     int tag=0;     int tag=0;
# Line 2045  int MeshAdapter::getTag(const std::strin Line 2230  int MeshAdapter::getTag(const std::strin
2230     return tag;     return tag;
2231  }  }
2232    
2233  bool MeshAdapter::isValidTagName(const std::string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2234  {  {
2235     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2236     return Finley_Mesh_isValidTagName(mesh,name.c_str());     return Finley_Mesh_isValidTagName(mesh,name.c_str());
2237  }  }
2238    
2239  std::string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2240  {  {
2241     stringstream temp;     stringstream temp;
2242     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2105  int MeshAdapter::getNumberOfTagsInUse(in Line 2290  int MeshAdapter::getNumberOfTagsInUse(in
2290    }    }
2291    return numTags;    return numTags;
2292  }  }
2293  int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  
2294    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2295  {  {
2296    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2297    index_t* tags=NULL;    index_t* tags=NULL;
# Line 2171  bool MeshAdapter::canTag(int functionSpa Line 2357  bool MeshAdapter::canTag(int functionSpa
2357    }    }
2358  }  }
2359    
2360    AbstractDomain::StatusType MeshAdapter::getStatus() const
2361    {
2362      Finley_Mesh* mesh=m_finleyMesh.get();
2363      return Finley_Mesh_getStatus(mesh);
2364    }
2365    
2366    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2367    {
2368      
2369      Finley_Mesh* mesh=m_finleyMesh.get();
2370      int order =-1;
2371      switch(functionSpaceCode) {
2372       case(Nodes):
2373       case(DegreesOfFreedom):
2374              order=mesh->approximationOrder;
2375              break;
2376       case(ReducedNodes):
2377       case(ReducedDegreesOfFreedom):
2378              order=mesh->reducedApproximationOrder;
2379              break;
2380       case(Elements):
2381       case(FaceElements):
2382       case(Points):
2383       case(ContactElementsZero):
2384       case(ContactElementsOne):
2385              order=mesh->integrationOrder;
2386              break;
2387       case(ReducedElements):
2388       case(ReducedFaceElements):
2389       case(ReducedContactElementsZero):
2390       case(ReducedContactElementsOne):
2391              order=mesh->reducedIntegrationOrder;
2392              break;
2393       default:
2394          stringstream temp;
2395          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2396          throw FinleyAdapterException(temp.str());
2397      }
2398      return order;
2399    }
2400    
2401    bool MeshAdapter::supportsContactElements() const
2402    {
2403      return true;
2404    }
2405    
2406    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2407    {
2408      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2409    }
2410    
2411    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2412    {
2413      Finley_ReferenceElementSet_dealloc(m_refSet);
2414    }
2415    
2416  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.2100  
changed lines
  Added in v.3379

  ViewVC Help
Powered by ViewVC 1.1.26