/[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 2551 by gross, Thu Jul 23 09:19:15 2009 UTC revision 3344 by caltinay, Thu Nov 11 23:26:52 2010 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2009 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  }  }
24    
25    #include <boost/python/import.hpp>
26    #include <boost/python/tuple.hpp>
27    
28  using namespace std;  using namespace std;
29  using namespace escript;  using namespace escript;
30    
# Line 85  int MeshAdapter::getMPIRank() const Line 84  int MeshAdapter::getMPIRank() const
84  }  }
85  void MeshAdapter::MPIBarrier() const  void MeshAdapter::MPIBarrier() const
86  {  {
87  #ifdef PASO_MPI  #ifdef ESYS_MPI
88     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
89  #endif  #endif
90     return;     return;
# Line 96  bool MeshAdapter::onMasterProcessor() co Line 95  bool MeshAdapter::onMasterProcessor() co
95  }  }
96    
97    
98    #ifdef ESYS_MPI
99      MPI_Comm
100    #else
101      unsigned int
102    #endif
103    MeshAdapter::getMPIComm() const
104    {
105    #ifdef ESYS_MPI
106        return m_finleyMesh->MPIInfo->comm;
107    #else
108        return 0;
109    #endif
110    }
111    
112    
113  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
114     return m_finleyMesh.get();     return m_finleyMesh.get();
115  }  }
116    
117  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const string& fileName) const
118  {  {
119     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;
120     strcpy(fName,fileName.c_str());     strcpy(fName,fileName.c_str());
# Line 114  void MeshAdapter::Print_Mesh_Info(const Line 128  void MeshAdapter::Print_Mesh_Info(const
128     Finley_PrintMesh_Info(m_finleyMesh.get(), full);     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
129  }  }
130    
131  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const string& fileName) const
132  {  {
133  #ifdef USE_NETCDF  #ifdef USE_NETCDF
134     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 134  void MeshAdapter::dump(const std::string Line 148  void MeshAdapter::dump(const std::string
148     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
149     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
150     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
151  #ifdef PASO_MPI  #ifdef ESYS_MPI
152     MPI_Status status;     MPI_Status status;
153  #endif  #endif
154    
155  /* Incoming token indicates it's my turn to write */  /* Incoming token indicates it's my turn to write */
156  #ifdef PASO_MPI  #ifdef ESYS_MPI
157     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);
158  #endif  #endif
159    
160     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
161                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
162    
163     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
164     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
165     num_Tags = 0;     num_Tags = 0;
166     if (tag_map) {     while (tag_map) {
167        while (tag_map) {        num_Tags++;
168           num_Tags++;        tag_map=tag_map->next;
          tag_map=tag_map->next;  
       }  
169     }     }
170    
171     // NetCDF error handler     // NetCDF error handler
172     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
173     // Create the file.     // Create the file.
174     NcFile dataFile(newFileName, NcFile::Replace);     NcFile dataFile(newFileName, NcFile::Replace);
175       string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
176     // check if writing was successful     // check if writing was successful
177     if (!dataFile.is_valid())     if (!dataFile.is_valid())
178        throw DataException("Error - MeshAdapter::dump: opening of NetCDF file for output failed: " + *newFileName);        throw DataException(msgPrefix+"Open file for output");
179    
180     // Define dimensions (num_Elements and dim_Elements are identical, dim_Elements only appears if > 0)     // Define dimensions (num_Elements and dim_Elements are identical,
181       // dim_Elements only appears if > 0)
182     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
183        throw DataException("Error - MeshAdapter::dump: appending dimension numNodes to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numNodes)");
184     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
185        throw DataException("Error - MeshAdapter::dump: appending dimension numDim to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numDim)");
186     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)) )
187        throw DataException("Error - MeshAdapter::dump: appending dimension mpi_size to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(mpi_size)");
188     if (num_Elements>0)     if (num_Elements>0)
189        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
190           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements)");
191     if (num_FaceElements>0)     if (num_FaceElements>0)
192        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
193           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
194     if (num_ContactElements>0)     if (num_ContactElements>0)
195        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
196           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements)");
197     if (num_Points>0)     if (num_Points>0)
198        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
199           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Points to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Points)");
200     if (num_Elements>0)     if (num_Elements>0)
201        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
202           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
203     if (num_FaceElements>0)     if (num_FaceElements>0)
204        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
205           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
206     if (num_ContactElements>0)     if (num_ContactElements>0)
207        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
208           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements_numNodes)");
209     if (num_Tags>0)     if (num_Tags>0)
210        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
211           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Tags)");
212    
213     // Attributes: MPI size, MPI rank, Name, order, reduced_order     // Attributes: MPI size, MPI rank, Name, order, reduced_order
214     if (!dataFile.add_att("mpi_size", mpi_size) )     if (!dataFile.add_att("mpi_size", mpi_size) )
215        throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_size)");
216     if (!dataFile.add_att("mpi_rank", mpi_rank) )     if (!dataFile.add_att("mpi_rank", mpi_rank) )
217        throw DataException("Error - MeshAdapter::dump: appending mpi_rank to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_rank)");
218     if (!dataFile.add_att("Name",mesh->Name) )     if (!dataFile.add_att("Name",mesh->Name) )
219        throw DataException("Error - MeshAdapter::dump: appending Name to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Name)");
220     if (!dataFile.add_att("numDim",numDim) )     if (!dataFile.add_att("numDim",numDim) )
221        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
222     if (!dataFile.add_att("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
223        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
224     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
225        throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(reduced_order)");
226     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
227        throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(numNodes)");
228     if (!dataFile.add_att("num_Elements",num_Elements) )     if (!dataFile.add_att("num_Elements",num_Elements) )
229        throw DataException("Error - MeshAdapter::dump: appending num_Elements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements)");
230     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
231        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements)");
232     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
233        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements)");
234     if (!dataFile.add_att("num_Points",num_Points) )     if (!dataFile.add_att("num_Points",num_Points) )
235        throw DataException("Error - MeshAdapter::dump: appending num_Points to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Points)");
236     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
237        throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
238     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
239        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
240     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
241        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements_numNodes)");
242     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
243        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Elements_TypeId)");
244     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
245        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
246     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
247        throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(ContactElements_TypeId)");
248     if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
249        throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Points_TypeId)");
250     if (!dataFile.add_att("num_Tags", num_Tags) )     if (!dataFile.add_att("num_Tags", num_Tags) )
251        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Tags)");
252    
253     // // // // // Nodes // // // // //     // // // // // Nodes // // // // //
254    
255     // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)     // Nodes nodeDistribution
256       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
257          throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
258       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
259       if (! (ids->put(int_ptr, mpi_size+1)) )
260          throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
261    
262       // Nodes degreesOfFreedomDistribution
263       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
264          throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
265       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
266       if (! (ids->put(int_ptr, mpi_size+1)) )
267          throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
268    
269       // Only write nodes if non-empty because NetCDF doesn't like empty arrays
270       // (it treats them as NC_UNLIMITED)
271     if (numNodes>0) {     if (numNodes>0) {
272    
273        // Nodes Id        // Nodes Id
274        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
275           throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Id)");
276        int_ptr = &mesh->Nodes->Id[0];        int_ptr = &mesh->Nodes->Id[0];
277        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
278           throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Id)");
279    
280        // Nodes Tag        // Nodes Tag
281        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
282           throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Tag)");
283        int_ptr = &mesh->Nodes->Tag[0];        int_ptr = &mesh->Nodes->Tag[0];
284        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
285           throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Tag)");
286    
287        // Nodes gDOF        // Nodes gDOF
288        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
289           throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
290        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
291        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
292           throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gDOF)");
293    
294        // Nodes global node index        // Nodes global node index
295        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
296           throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gNI)");
297        int_ptr = &mesh->Nodes->globalNodesIndex[0];        int_ptr = &mesh->Nodes->globalNodesIndex[0];
298        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
299           throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gNI)");
300    
301        // Nodes grDof        // Nodes grDof
302        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
303           throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
304        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
305        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
306           throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grDfI)");
307    
308        // Nodes grNI        // Nodes grNI
309        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
310           throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grNI)");
311        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
312        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
313           throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grNI)");
314    
315        // Nodes Coordinates        // Nodes Coordinates
316        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
317           throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
318        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
319           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);  
320    
321     }     }
322    
# Line 311  void MeshAdapter::dump(const std::string Line 326  void MeshAdapter::dump(const std::string
326    
327        // Elements_Id        // Elements_Id
328        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
329           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Id)");
330        int_ptr = &mesh->Elements->Id[0];        int_ptr = &mesh->Elements->Id[0];
331        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
332           throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Id)");
333    
334        // Elements_Tag        // Elements_Tag
335        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
336           throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Tag)");
337        int_ptr = &mesh->Elements->Tag[0];        int_ptr = &mesh->Elements->Tag[0];
338        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
339           throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Tag)");
340    
341        // Elements_Owner        // Elements_Owner
342        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
343           throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Owner)");
344        int_ptr = &mesh->Elements->Owner[0];        int_ptr = &mesh->Elements->Owner[0];
345        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
346           throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Owner)");
347    
348        // Elements_Color        // Elements_Color
349        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
350           throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Color)");
351        int_ptr = &mesh->Elements->Color[0];        int_ptr = &mesh->Elements->Color[0];
352        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
353           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Color)");
354    
355        // Elements_Nodes        // Elements_Nodes
356        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
357           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Nodes)");
358        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
359           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Nodes)");
360    
361     }     }
362    
# Line 351  void MeshAdapter::dump(const std::string Line 366  void MeshAdapter::dump(const std::string
366    
367        // FaceElements_Id        // FaceElements_Id
368        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
369           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Id)");
370        int_ptr = &mesh->FaceElements->Id[0];        int_ptr = &mesh->FaceElements->Id[0];
371        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
372           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Id)");
373    
374        // FaceElements_Tag        // FaceElements_Tag
375        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
376           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
377        int_ptr = &mesh->FaceElements->Tag[0];        int_ptr = &mesh->FaceElements->Tag[0];
378        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
379           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Tag)");
380    
381        // FaceElements_Owner        // FaceElements_Owner
382        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
383           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
384        int_ptr = &mesh->FaceElements->Owner[0];        int_ptr = &mesh->FaceElements->Owner[0];
385        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
386           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Owner)");
387    
388        // FaceElements_Color        // FaceElements_Color
389        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
390           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Color)");
391        int_ptr = &mesh->FaceElements->Color[0];        int_ptr = &mesh->FaceElements->Color[0];
392        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
393           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Color)");
394    
395        // FaceElements_Nodes        // FaceElements_Nodes
396        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
397           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
398        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
399           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Nodes)");
400    
401     }     }
402    
# Line 391  void MeshAdapter::dump(const std::string Line 406  void MeshAdapter::dump(const std::string
406    
407        // ContactElements_Id        // ContactElements_Id
408        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
409           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Id)");
410        int_ptr = &mesh->ContactElements->Id[0];        int_ptr = &mesh->ContactElements->Id[0];
411        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
412           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Id)");
413    
414        // ContactElements_Tag        // ContactElements_Tag
415        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
416           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Tag)");
417        int_ptr = &mesh->ContactElements->Tag[0];        int_ptr = &mesh->ContactElements->Tag[0];
418        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
419           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Tag)");
420    
421        // ContactElements_Owner        // ContactElements_Owner
422        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
423           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Owner)");
424        int_ptr = &mesh->ContactElements->Owner[0];        int_ptr = &mesh->ContactElements->Owner[0];
425        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
426           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Owner)");
427    
428        // ContactElements_Color        // ContactElements_Color
429        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
430           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Color)");
431        int_ptr = &mesh->ContactElements->Color[0];        int_ptr = &mesh->ContactElements->Color[0];
432        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
433           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Color)");
434    
435        // ContactElements_Nodes        // ContactElements_Nodes
436        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
437           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Nodes)");
438        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
439           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Nodes)");
440    
441     }     }
442    
# Line 433  void MeshAdapter::dump(const std::string Line 448  void MeshAdapter::dump(const std::string
448    
449        // Points_Id        // Points_Id
450        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
451           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Id)");
452        int_ptr = &mesh->Points->Id[0];        int_ptr = &mesh->Points->Id[0];
453        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
454           throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Id)");
455    
456        // Points_Tag        // Points_Tag
457        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
458           throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Tag)");
459        int_ptr = &mesh->Points->Tag[0];        int_ptr = &mesh->Points->Tag[0];
460        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
461           throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Tag)");
462    
463        // Points_Owner        // Points_Owner
464        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
465           throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Owner)");
466        int_ptr = &mesh->Points->Owner[0];        int_ptr = &mesh->Points->Owner[0];
467        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
468           throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Owner)");
469    
470        // Points_Color        // Points_Color
471        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
472           throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Color)");
473        int_ptr = &mesh->Points->Color[0];        int_ptr = &mesh->Points->Color[0];
474        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
475           throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Color)");
476    
477        // Points_Nodes        // Points_Nodes
478        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
479        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
480           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Nodes)");
481        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
482           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Nodes)");
483    
484     }     }
485    
# Line 488  void MeshAdapter::dump(const std::string Line 503  void MeshAdapter::dump(const std::string
503    
504        // Tags_keys        // Tags_keys
505        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
506           throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Tags_keys)");
507        int_ptr = &Tags_keys[0];        int_ptr = &Tags_keys[0];
508        if (! (ids->put(int_ptr, num_Tags)) )        if (! (ids->put(int_ptr, num_Tags)) )
509           throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Tags_keys)");
510    
511        // Tags_names_*        // Tags_names_*
512        // This is an array of strings, it should be stored as an array but instead I have hacked in one attribute per string        // This is an array of strings, it should be stored as an array but
513        // 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
514          // manual doesn't tell how to do an array of strings
515        tag_map = mesh->TagMap;        tag_map = mesh->TagMap;
516        if (tag_map) {        if (tag_map) {
517           int i = 0;           int i = 0;
518           while (tag_map) {           while (tag_map) {
519              sprintf(name_temp, "Tags_name_%d", i);              sprintf(name_temp, "Tags_name_%d", i);
520              if (!dataFile.add_att(name_temp, tag_map->name) )              if (!dataFile.add_att(name_temp, tag_map->name) )
521                 throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);                 throw DataException(msgPrefix+"add_att(Tags_names_XX)");
522              tag_map=tag_map->next;              tag_map=tag_map->next;
523              i++;              i++;
524           }           }
525        }        }
526    
527        TMPMEMFREE(Tags_keys);        TMPMEMFREE(Tags_keys);
   
528     }     }
529    
530  /* Send token to next MPI process so he can take his turn */  /* Send token to next MPI process so he can take his turn */
531  #ifdef PASO_MPI  #ifdef ESYS_MPI
532     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);
533  #endif  #endif
534    
# Line 678  pair<int,int> MeshAdapter::getDataShape( Line 693  pair<int,int> MeshAdapter::getDataShape(
693     case(Elements):     case(Elements):
694     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
695        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
696        numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
697     }     }
698     break;     break;
699     case(ReducedElements):     case(ReducedElements):
700     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
701        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
702        numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
703     }     }
704     break;     break;
705     case(FaceElements):     case(FaceElements):
706     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
707        numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
708        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
709     }     }
710     break;     break;
711     case(ReducedFaceElements):     case(ReducedFaceElements):
712     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
713        numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
714        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
715     }     }
716     break;     break;
# Line 707  pair<int,int> MeshAdapter::getDataShape( Line 722  pair<int,int> MeshAdapter::getDataShape(
722     break;     break;
723     case(ContactElementsZero):     case(ContactElementsZero):
724     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
725        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
726        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
727     }     }
728     break;     break;
729     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
730     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
731        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
732        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
733     }     }
734     break;     break;
735     case(ContactElementsOne):     case(ContactElementsOne):
736     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
737        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
738        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
739     }     }
740     break;     break;
741     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
742     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
743        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
744        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
745     }     }
746     break;     break;
# Line 754  pair<int,int> MeshAdapter::getDataShape( Line 769  pair<int,int> MeshAdapter::getDataShape(
769  // adds linear PDE of second order into a given stiffness matrix and right hand side:  // adds linear PDE of second order into a given stiffness matrix and right hand side:
770  //  //
771  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
772                                   SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
773                                   const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                   const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
774                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
775                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact) const
776  {  {
777       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
778       if (smat==0)
779       {
780        throw FinleyAdapterException("finley only supports its own system matrices.");
781       }
782     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
783     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
784     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 773  void MeshAdapter::addPDEToSystem( Line 793  void MeshAdapter::addPDEToSystem(
793    
794     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
795    
796     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 );
797     checkFinleyError();     checkFinleyError();
798    
799     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 );
800     checkFinleyError();     checkFinleyError();
801    
802     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 );
803     checkFinleyError();     checkFinleyError();
804  }  }
805    
# Line 795  void  MeshAdapter::addPDEToLumpedSystem( Line 815  void  MeshAdapter::addPDEToLumpedSystem(
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);
818       checkFinleyError();
819      
820     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
   
821     checkFinleyError();     checkFinleyError();
822    
823  }  }
824    
825    
# Line 827  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 849  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 1026  void MeshAdapter::interpolateOnDomain(es Line 1055  void MeshAdapter::interpolateOnDomain(es
1055        case(Elements):        case(Elements):
1056        case(ReducedElements):        case(ReducedElements):
1057        if (getMPISize()>1) {        if (getMPISize()>1) {
1058           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1059           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1060           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1061        } else {        } else {
# Line 1036  void MeshAdapter::interpolateOnDomain(es Line 1065  void MeshAdapter::interpolateOnDomain(es
1065        case(FaceElements):        case(FaceElements):
1066        case(ReducedFaceElements):        case(ReducedFaceElements):
1067        if (getMPISize()>1) {        if (getMPISize()>1) {
1068           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1069           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1070           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1071        
# Line 1046  void MeshAdapter::interpolateOnDomain(es Line 1075  void MeshAdapter::interpolateOnDomain(es
1075        break;        break;
1076        case(Points):        case(Points):
1077        if (getMPISize()>1) {        if (getMPISize()>1) {
1078           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1079           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1080        } else {        } else {
1081           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
# Line 1057  void MeshAdapter::interpolateOnDomain(es Line 1086  void MeshAdapter::interpolateOnDomain(es
1086        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1087        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1088        if (getMPISize()>1) {        if (getMPISize()>1) {
1089           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1090           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1091           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1092        } else {        } else {
# Line 1095  void MeshAdapter::interpolateOnDomain(es Line 1124  void MeshAdapter::interpolateOnDomain(es
1124        case(Elements):        case(Elements):
1125        case(ReducedElements):        case(ReducedElements):
1126        if (getMPISize()>1) {        if (getMPISize()>1) {
1127           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1128           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1129           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1130        } else {        } else {
# Line 1105  void MeshAdapter::interpolateOnDomain(es Line 1134  void MeshAdapter::interpolateOnDomain(es
1134        case(FaceElements):        case(FaceElements):
1135        case(ReducedFaceElements):        case(ReducedFaceElements):
1136        if (getMPISize()>1) {        if (getMPISize()>1) {
1137           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1138           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1139           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1140        } else {        } else {
# Line 1114  void MeshAdapter::interpolateOnDomain(es Line 1143  void MeshAdapter::interpolateOnDomain(es
1143        break;        break;
1144        case(Points):        case(Points):
1145        if (getMPISize()>1) {        if (getMPISize()>1) {
1146           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1147           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1148           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1149        } else {        } else {
# Line 1126  void MeshAdapter::interpolateOnDomain(es Line 1155  void MeshAdapter::interpolateOnDomain(es
1155        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1156        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1157        if (getMPISize()>1) {        if (getMPISize()>1) {
1158           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1159           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1160           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1161        } else {        } else {
# Line 1163  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 1244  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 1257  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 1294  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 1331  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 1453  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     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1486         tmp = new_x.getDataC();         tmp = new_x.getDataC();
1487         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1488     } else {     } else {
1489         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1490         tmp = new_x_inter.getDataC();         tmp = new_x_inter.getDataC();
1491         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1492     }     }
# Line 1478  void MeshAdapter::extractArgsFromDict(co Line 1507  void MeshAdapter::extractArgsFromDict(co
1507    
1508     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1509     for (int i=0; i<numData; ++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: Data must be defined on same Domain");           throw FinleyAdapterException("Error: Data must be defined on same Domain");
# Line 1492  void MeshAdapter::extractArgsFromDict(co Line 1521  void MeshAdapter::extractArgsFromDict(co
1521  //  //
1522  // saves mesh and optionally data arrays in openDX format  // saves mesh and optionally data arrays in openDX format
1523  //  //
1524  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1525  {  {
1526     int num_data;     int num_data;
1527     char **names;     char **names;
# Line 1518  void MeshAdapter::saveDX(const std::stri Line 1547  void MeshAdapter::saveDX(const std::stri
1547  //  //
1548  // saves mesh and optionally data arrays in VTK format  // saves mesh and optionally data arrays in VTK format
1549  //  //
1550  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg,  const std::string& metadata, const std::string& metadata_schema) const  void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1551  {  {
1552     int num_data;      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("saveVTK");
1553     char **names;      pySaveVTK(*boost::python::make_tuple(filename,
1554     escriptDataC *data;                 const_cast<MeshAdapter*>(this)->getPtr(),
1555     escriptDataC **ptr_data;                 metadata, metadata_schema), **arg);
1556    }
1557    
1558    bool MeshAdapter::ownSample(int fs_code, index_t id) const
1559    {
1560    #ifdef ESYS_MPI
1561        index_t myFirstNode=0, myLastNode=0, k=0;
1562        index_t* globalNodeIndex=0;
1563        Finley_Mesh* mesh_p=m_finleyMesh.get();
1564        if (fs_code == FINLEY_REDUCED_NODES)
1565        {
1566        myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1567        myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1568        globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1569        }
1570        else
1571        {
1572        myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1573        myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1574        globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1575        }
1576        k=globalNodeIndex[id];
1577        return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1578    #endif
1579        return true;
1580    }
1581    
    extractArgsFromDict(arg, num_data, names, data, ptr_data);  
    Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());  
    checkFinleyError();  
1582    
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0; i<num_data; i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
 }  
1583    
1584  //  //
1585  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1586  //  //
1587  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  ASM_ptr MeshAdapter::newSystemMatrix(
1588                                                   const int row_blocksize,                                                   const int row_blocksize,
1589                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1590                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1589  SystemMatrixAdapter MeshAdapter::newSyst Line 1631  SystemMatrixAdapter MeshAdapter::newSyst
1631     }     }
1632     checkPasoError();     checkPasoError();
1633     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1634     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1635       return ASM_ptr(sma);
1636    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1637  }  }
1638    
1639  //  //
1640  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1641  //  //
1642  TransportProblemAdapter MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
1643                                                           const double theta,                                                           const bool useBackwardEuler,
1644                                                           const int blocksize,                                                           const int blocksize,
1645                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1646                                                           const int type) const                                                           const int type) const
# Line 1618  TransportProblemAdapter MeshAdapter::new Line 1662  TransportProblemAdapter MeshAdapter::new
1662    
1663     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1664     checkFinleyError();     checkFinleyError();
1665     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1666     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1667     checkPasoError();     checkPasoError();
1668     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1669     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1670       return ATP_ptr(tpa);
1671    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1672  }  }
1673    
1674  //  //
# Line 1661  bool MeshAdapter::isCellOriented(int fun Line 1707  bool MeshAdapter::isCellOriented(int fun
1707     return false;     return false;
1708  }  }
1709    
1710    bool
1711    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1712    {
1713       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1714        class 1: DOF <-> Nodes
1715        class 2: ReducedDOF <-> ReducedNodes
1716        class 3: Points
1717        class 4: Elements
1718        class 5: ReducedElements
1719        class 6: FaceElements
1720        class 7: ReducedFaceElements
1721        class 8: ContactElementZero <-> ContactElementOne
1722        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1723    
1724       There is also a set of lines. Interpolation is possible down a line but not between lines.
1725       class 1 and 2 belong to all lines so aren't considered.
1726        line 0: class 3
1727        line 1: class 4,5
1728        line 2: class 6,7
1729        line 3: class 8,9
1730    
1731       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1732       eg hasnodes is true if we have at least one instance of Nodes.
1733       */
1734        if (fs.empty())
1735        {
1736            return false;
1737        }
1738        vector<int> hasclass(10);
1739        vector<int> hasline(4);
1740        bool hasnodes=false;
1741        bool hasrednodes=false;
1742        bool hascez=false;
1743        bool hasrcez=false;
1744        for (int i=0;i<fs.size();++i)
1745        {
1746        switch(fs[i])
1747        {
1748        case(Nodes):   hasnodes=true;   // no break is deliberate
1749        case(DegreesOfFreedom):
1750            hasclass[1]=1;
1751            break;
1752        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1753        case(ReducedDegreesOfFreedom):
1754            hasclass[2]=1;
1755            break;
1756        case(Points):
1757            hasline[0]=1;
1758            hasclass[3]=1;
1759            break;
1760        case(Elements):
1761            hasclass[4]=1;
1762            hasline[1]=1;
1763            break;
1764        case(ReducedElements):
1765            hasclass[5]=1;
1766            hasline[1]=1;
1767            break;
1768        case(FaceElements):
1769            hasclass[6]=1;
1770            hasline[2]=1;
1771            break;
1772        case(ReducedFaceElements):
1773            hasclass[7]=1;
1774            hasline[2]=1;
1775            break;
1776        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1777        case(ContactElementsOne):
1778            hasclass[8]=1;
1779            hasline[3]=1;
1780            break;
1781        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1782        case(ReducedContactElementsOne):
1783            hasclass[9]=1;
1784            hasline[3]=1;
1785            break;
1786        default:
1787            return false;
1788        }
1789        }
1790        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1791        // fail if we have more than one leaf group
1792    
1793        if (totlines>1)
1794        {
1795        return false;   // there are at least two branches we can't interpolate between
1796        }
1797        else if (totlines==1)
1798        {
1799        if (hasline[0]==1)      // we have points
1800        {
1801            resultcode=Points;
1802        }
1803        else if (hasline[1]==1)
1804        {
1805            if (hasclass[5]==1)
1806            {
1807            resultcode=ReducedElements;
1808            }
1809            else
1810            {
1811            resultcode=Elements;
1812            }
1813        }
1814        else if (hasline[2]==1)
1815        {
1816            if (hasclass[7]==1)
1817            {
1818            resultcode=ReducedFaceElements;
1819            }
1820            else
1821            {
1822            resultcode=FaceElements;
1823            }
1824        }
1825        else    // so we must be in line3
1826        {
1827            if (hasclass[9]==1)
1828            {
1829            // need something from class 9
1830            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1831            }
1832            else
1833            {
1834            // something from class 8
1835            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1836            }
1837        }
1838        }
1839        else    // totlines==0
1840        {
1841        if (hasclass[2]==1)
1842        {
1843            // something from class 2
1844            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1845        }
1846        else
1847        {   // something from class 1
1848            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1849        }
1850        }
1851        return true;
1852    }
1853    
1854  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1855  {  {
1856     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1857     case(Nodes):     case(Nodes):
1858     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1859     case(Nodes):      case(Nodes):
1860     case(ReducedNodes):      case(ReducedNodes):
1861     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1862     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1863     case(Elements):      case(Elements):
1864     case(ReducedElements):      case(ReducedElements):
1865     case(FaceElements):      case(FaceElements):
1866     case(ReducedFaceElements):      case(ReducedFaceElements):
1867     case(Points):      case(Points):
1868     case(ContactElementsZero):      case(ContactElementsZero):
1869     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1870     case(ContactElementsOne):      case(ContactElementsOne):
1871     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1872     return true;      return true;
1873     default:      default:
1874        stringstream temp;            stringstream temp;
1875        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;
1876        throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1877     }     }
1878     break;     break;
1879     case(ReducedNodes):     case(ReducedNodes):
1880     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1881     case(ReducedNodes):      case(ReducedNodes):
1882     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1883     case(Elements):      case(Elements):
1884     case(ReducedElements):      case(ReducedElements):
1885     case(FaceElements):      case(FaceElements):
1886     case(ReducedFaceElements):      case(ReducedFaceElements):
1887     case(Points):      case(Points):
1888     case(ContactElementsZero):      case(ContactElementsZero):
1889     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1890     case(ContactElementsOne):      case(ContactElementsOne):
1891     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1892     return true;      return true;
1893     case(Nodes):      case(Nodes):
1894     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1895     return false;      return false;
1896     default:      default:
1897        stringstream temp;          stringstream temp;
1898        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;
1899        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1900     }     }
1901     break;     break;
1902     case(Elements):     case(Elements):
1903     if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1904        return true;        return true;
1905     } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1906        return true;        return true;
1907     } else {          } else {
1908        return false;            return false;
1909     }          }
1910     case(ReducedElements):     case(ReducedElements):
1911     if (functionSpaceType_target==ReducedElements) {      if (functionSpaceType_target==ReducedElements) {
1912        return true;        return true;
1913     } else {      } else {
1914        return false;            return false;
1915     }      }
1916     case(FaceElements):     case(FaceElements):
1917     if (functionSpaceType_target==FaceElements) {      if (functionSpaceType_target==FaceElements) {
1918        return true;              return true;
1919     } else if (functionSpaceType_target==ReducedFaceElements) {      } else if (functionSpaceType_target==ReducedFaceElements) {
1920        return true;              return true;
1921     } else {      } else {
1922        return false;              return false;
1923     }      }
1924     case(ReducedFaceElements):     case(ReducedFaceElements):
1925     if (functionSpaceType_target==ReducedFaceElements) {      if (functionSpaceType_target==ReducedFaceElements) {
1926        return true;              return true;
1927     } else {      } else {
1928        return false;          return false;
1929     }      }
1930     case(Points):     case(Points):
1931     if (functionSpaceType_target==Points) {      if (functionSpaceType_target==Points) {
1932        return true;              return true;
1933     } else {      } else {
1934        return false;              return false;
1935     }      }
1936     case(ContactElementsZero):     case(ContactElementsZero):
1937     case(ContactElementsOne):     case(ContactElementsOne):
1938     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1939        return true;              return true;
1940     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1941        return true;              return true;
1942     } else {      } else {
1943        return false;              return false;
1944     }      }
1945     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1946     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1947     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1948        return true;              return true;
1949     } else {      } else {
1950        return false;              return false;
1951     }      }
1952     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1953     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1954     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1955     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1956     case(Nodes):      case(Nodes):
1957     case(ReducedNodes):      case(ReducedNodes):
1958     case(Elements):      case(Elements):
1959     case(ReducedElements):      case(ReducedElements):
1960     case(Points):      case(Points):
1961     case(FaceElements):      case(FaceElements):
1962     case(ReducedFaceElements):      case(ReducedFaceElements):
1963     case(ContactElementsZero):      case(ContactElementsZero):
1964     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1965     case(ContactElementsOne):      case(ContactElementsOne):
1966     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1967     return true;      return true;
1968     default:      default:
1969        stringstream temp;          stringstream temp;
1970        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;
1971        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1972     }      }
1973     break;      break;
1974     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1975     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
1976     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1977     case(ReducedNodes):      case(ReducedNodes):
1978     case(Elements):      case(Elements):
1979     case(ReducedElements):      case(ReducedElements):
1980     case(FaceElements):      case(FaceElements):
1981     case(ReducedFaceElements):      case(ReducedFaceElements):
1982     case(Points):      case(Points):
1983     case(ContactElementsZero):      case(ContactElementsZero):
1984     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1985     case(ContactElementsOne):      case(ContactElementsOne):
1986     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1987     return true;      return true;
1988     case(Nodes):      case(Nodes):
1989     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1990     return false;      return false;
1991     default:      default:
1992        stringstream temp;          stringstream temp;
1993        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;
1994        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1995     }      }
1996     break;      break;
1997     default:     default:
1998        stringstream temp;        stringstream temp;
1999        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 1844  int MeshAdapter::getSystemMatrixTypeId(c Line 2034  int MeshAdapter::getSystemMatrixTypeId(c
2034  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
2035  {  {
2036     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2037     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2038     checkPasoError();     checkPasoError();
2039     return out;     return out;
2040  }  }
2041    
2042  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2043  {  {
2044     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2045  }  }
2046    
2047  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2048  {  {
2049     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2050  }  }
2051    
2052  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2053  {  {
2054     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2055  }  }
2056    
2057  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
# Line 2023  void MeshAdapter::setTags(const int func Line 2213  void MeshAdapter::setTags(const int func
2213     return;     return;
2214  }  }
2215    
2216  void MeshAdapter::setTagMap(const std::string& name,  int tag)  void MeshAdapter::setTagMap(const string& name,  int tag)
2217  {  {
2218     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2219     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
# Line 2031  void MeshAdapter::setTagMap(const std::s Line 2221  void MeshAdapter::setTagMap(const std::s
2221     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2222  }  }
2223    
2224  int MeshAdapter::getTag(const std::string& name) const  int MeshAdapter::getTag(const string& name) const
2225  {  {
2226     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2227     int tag=0;     int tag=0;
# Line 2041  int MeshAdapter::getTag(const std::strin Line 2231  int MeshAdapter::getTag(const std::strin
2231     return tag;     return tag;
2232  }  }
2233    
2234  bool MeshAdapter::isValidTagName(const std::string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2235  {  {
2236     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2237     return Finley_Mesh_isValidTagName(mesh,name.c_str());     return Finley_Mesh_isValidTagName(mesh,name.c_str());
2238  }  }
2239    
2240  std::string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2241  {  {
2242     stringstream temp;     stringstream temp;
2243     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2174  AbstractDomain::StatusType MeshAdapter:: Line 2364  AbstractDomain::StatusType MeshAdapter::
2364    return Finley_Mesh_getStatus(mesh);    return Finley_Mesh_getStatus(mesh);
2365  }  }
2366    
2367    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2368    {
2369      
2370      Finley_Mesh* mesh=m_finleyMesh.get();
2371      int order =-1;
2372      switch(functionSpaceCode) {
2373       case(Nodes):
2374       case(DegreesOfFreedom):
2375              order=mesh->approximationOrder;
2376              break;
2377       case(ReducedNodes):
2378       case(ReducedDegreesOfFreedom):
2379              order=mesh->reducedApproximationOrder;
2380              break;
2381       case(Elements):
2382       case(FaceElements):
2383       case(Points):
2384       case(ContactElementsZero):
2385       case(ContactElementsOne):
2386              order=mesh->integrationOrder;
2387              break;
2388       case(ReducedElements):
2389       case(ReducedFaceElements):
2390       case(ReducedContactElementsZero):
2391       case(ReducedContactElementsOne):
2392              order=mesh->reducedIntegrationOrder;
2393              break;
2394       default:
2395          stringstream temp;
2396          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2397          throw FinleyAdapterException(temp.str());
2398      }
2399      return order;
2400    }
2401    
2402    bool MeshAdapter::supportsContactElements() const
2403    {
2404      return true;
2405    }
2406    
2407    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2408    {
2409      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2410    }
2411    
2412    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2413    {
2414      Finley_ReferenceElementSet_dealloc(m_refSet);
2415    }
2416    
2417  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.2551  
changed lines
  Added in v.3344

  ViewVC Help
Powered by ViewVC 1.1.26