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

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

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

revision 1990 by ksteube, Fri Nov 7 04:19:07 2008 UTC revision 3344 by caltinay, Thu Nov 11 23:26:52 2010 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2010 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 18  Line 18 
18  #ifdef USE_NETCDF  #ifdef USE_NETCDF
19  #include <netcdfcpp.h>  #include <netcdfcpp.h>
20  #endif  #endif
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #include "paso/Paso_MPI.h"  
 #endif  
21  extern "C" {  extern "C" {
22  #include "escript/blocktimer.h"  #include "esysUtils/blocktimer.h"
23  }  }
 #include <vector>  
24    
25  #define IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH  256  #include <boost/python/import.hpp>
26    #include <boost/python/tuple.hpp>
27    
28  using namespace std;  using namespace std;
29  using namespace escript;  using namespace escript;
# Line 88  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 99  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 112  void MeshAdapter::write(const std::strin Line 123  void MeshAdapter::write(const std::strin
123     TMPMEMFREE(fName);     TMPMEMFREE(fName);
124  }  }
125    
126  void MeshAdapter::Print_Mesh_Info(const bool full=false) const  void MeshAdapter::Print_Mesh_Info(const bool full) const
127  {  {
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 137  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 314  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 354  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 394  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 436  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 491  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 646  int MeshAdapter::getDiracDeltaFunctionCo Line 658  int MeshAdapter::getDiracDeltaFunctionCo
658  //  //
659  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
660  {  {
661     int numDim=Finley_Mesh_getDim(m_finleyMesh.get());     Finley_Mesh* mesh=m_finleyMesh.get();
662       int numDim=Finley_Mesh_getDim(mesh);
663     checkFinleyError();     checkFinleyError();
664     return numDim;     return numDim;
665  }  }
# Line 680  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 709  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 756  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 775  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 797  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 829  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 851  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 883  void MeshAdapter::interpolateOnDomain(es Line 910  void MeshAdapter::interpolateOnDomain(es
910     escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
911     switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
912     case(Nodes):     case(Nodes):
913     switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
914     case(Nodes):        case(Nodes):
915     case(ReducedNodes):        case(ReducedNodes):
916     case(DegreesOfFreedom):        case(DegreesOfFreedom):
917     case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
918     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
    break;  
    case(Elements):  
    case(ReducedElements):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);  
    break;  
    case(FaceElements):  
    case(ReducedFaceElements):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);  
    break;  
    case(Points):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);  
    break;  
    case(ContactElementsZero):  
    case(ReducedContactElementsZero):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
919        break;        break;
920     }        case(Elements):
921     break;        case(ReducedElements):
922     case(ReducedNodes):        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
923     switch(target.getFunctionSpace().getTypeCode()) {        break;
924     case(Nodes):        case(FaceElements):
925     case(ReducedNodes):        case(ReducedFaceElements):
926     case(DegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
927     case(ReducedDegreesOfFreedom):        break;
928     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(Points):
929     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
930     case(Elements):        break;
931     case(ReducedElements):        case(ContactElementsZero):
932     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        case(ReducedContactElementsZero):
933     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
934     case(FaceElements):        break;
935     case(ReducedFaceElements):        case(ContactElementsOne):
936     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        case(ReducedContactElementsOne):
937     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
938     case(Points):        break;
939     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        default:
940     break;           stringstream temp;
941     case(ContactElementsZero):           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
942     case(ReducedContactElementsZero):           throw FinleyAdapterException(temp.str());
943     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           break;
944     break;        }
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
945        break;        break;
    }  
    break;  
    case(Elements):  
    if (target.getFunctionSpace().getTypeCode()==Elements) {  
       Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
       Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on elements possible.");  
    }  
    break;  
    case(ReducedElements):  
    if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
       Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");  
    }  
    break;  
    case(FaceElements):  
    if (target.getFunctionSpace().getTypeCode()==FaceElements) {  
       Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {  
       Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");  
    }  
    break;  
    case(ReducedFaceElements):  
    if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {  
       Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");  
    }  
    break;  
    case(Points):  
    if (target.getFunctionSpace().getTypeCode()==Points) {  
       Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on points possible.");  
    }  
    break;  
    case(ContactElementsZero):  
    case(ContactElementsOne):  
    if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {  
       Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {  
       Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");  
    }  
    break;  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {  
       Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");  
    }  
    break;  
    case(DegreesOfFreedom):        
    switch(target.getFunctionSpace().getTypeCode()) {  
    case(ReducedDegreesOfFreedom):  
    case(DegreesOfFreedom):  
    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    break;  
   
    case(Nodes):  
946     case(ReducedNodes):     case(ReducedNodes):
947     if (getMPISize()>1) {        switch(target.getFunctionSpace().getTypeCode()) {
948        escript::Data temp=escript::Data(in);        case(Nodes):
949        temp.expand();        case(ReducedNodes):
950        escriptDataC _in2 = temp.getDataC();        case(DegreesOfFreedom):
951        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);        case(ReducedDegreesOfFreedom):
    } else {  
952        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
953     }        break;
954     break;        case(Elements):
955     case(Elements):        case(ReducedElements):
    case(ReducedElements):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);  
    } else {  
956        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
957     }        break;
958     break;        case(FaceElements):
959     case(FaceElements):        case(ReducedFaceElements):
    case(ReducedFaceElements):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);  
   
    } else {  
960        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
961     }        break;
962     break;        case(Points):
    case(Points):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
    } else {  
963        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
964     }        break;
965     break;        case(ContactElementsZero):
966     case(ContactElementsZero):        case(ReducedContactElementsZero):
    case(ContactElementsOne):  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);  
    } else {  
967        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
    }  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
968        break;        break;
969     }        case(ContactElementsOne):
970     break;        case(ReducedContactElementsOne):
971     case(ReducedDegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
972     switch(target.getFunctionSpace().getTypeCode()) {        break;
973     case(Nodes):        default:
974     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           stringstream temp;
975     break;           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
976     case(ReducedNodes):           throw FinleyAdapterException(temp.str());
977     if (getMPISize()>1) {           break;
978        escript::Data temp=escript::Data(in);        }
979        temp.expand();        break;
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);  
    } else {  
       Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    }  
    break;  
    case(DegreesOfFreedom):  
    throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");  
    break;  
    case(ReducedDegreesOfFreedom):  
    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    break;  
980     case(Elements):     case(Elements):
981          if (target.getFunctionSpace().getTypeCode()==Elements) {
982             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
983          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
984             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
985          } else {
986             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
987          }
988          break;
989     case(ReducedElements):     case(ReducedElements):
990     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
991        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
992        escriptDataC _in2 = temp.getDataC();        } else {
993        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
994     } else {        }
995        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        break;
    }  
    break;  
996     case(FaceElements):     case(FaceElements):
997          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
998             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
999          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1000             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1001          } else {
1002             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1003          }
1004          break;
1005     case(ReducedFaceElements):     case(ReducedFaceElements):
1006     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1007        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1008        escriptDataC _in2 = temp.getDataC();        } else {
1009        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
1010     } else {        }
1011        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
    }  
    break;  
1012     case(Points):     case(Points):
1013     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==Points) {
1014        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1015        escriptDataC _in2 = temp.getDataC();        } else {
1016        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1017     } else {        }
1018        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        break;
    }  
    break;  
1019     case(ContactElementsZero):     case(ContactElementsZero):
1020     case(ContactElementsOne):     case(ContactElementsOne):
1021          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1022             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1023          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1024             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1025          } else {
1026             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1027          }
1028          break;
1029     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1030     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1031     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1032        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1033        escriptDataC _in2 = temp.getDataC();        } else {
1034        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1035     } else {        }
1036        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        break;
1037     }     case(DegreesOfFreedom):      
1038     break;        switch(target.getFunctionSpace().getTypeCode()) {
1039     default:        case(ReducedDegreesOfFreedom):
1040        stringstream temp;        case(DegreesOfFreedom):
1041        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1042        throw FinleyAdapterException(temp.str());        break;
1043      
1044          case(Nodes):
1045          case(ReducedNodes):
1046          if (getMPISize()>1) {
1047             escript::Data temp=escript::Data(in);
1048             temp.expand();
1049             escriptDataC _in2 = temp.getDataC();
1050             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1051          } else {
1052             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1053          }
1054          break;
1055          case(Elements):
1056          case(ReducedElements):
1057          if (getMPISize()>1) {
1058             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1059             escriptDataC _in2 = temp.getDataC();
1060             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1061          } else {
1062             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1063          }
1064          break;
1065          case(FaceElements):
1066          case(ReducedFaceElements):
1067          if (getMPISize()>1) {
1068             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1069             escriptDataC _in2 = temp.getDataC();
1070             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1071      
1072          } else {
1073             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1074          }
1075          break;
1076          case(Points):
1077          if (getMPISize()>1) {
1078             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1079             escriptDataC _in2 = temp.getDataC();
1080          } else {
1081             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1082          }
1083          break;
1084          case(ContactElementsZero):
1085          case(ContactElementsOne):
1086          case(ReducedContactElementsZero):
1087          case(ReducedContactElementsOne):
1088          if (getMPISize()>1) {
1089             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1090             escriptDataC _in2 = temp.getDataC();
1091             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1092          } else {
1093             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1094          }
1095          break;
1096          default:
1097             stringstream temp;
1098             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1099             throw FinleyAdapterException(temp.str());
1100             break;
1101          }
1102          break;
1103       case(ReducedDegreesOfFreedom):
1104          switch(target.getFunctionSpace().getTypeCode()) {
1105          case(Nodes):
1106          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1107          break;
1108          case(ReducedNodes):
1109          if (getMPISize()>1) {
1110             escript::Data temp=escript::Data(in);
1111             temp.expand();
1112             escriptDataC _in2 = temp.getDataC();
1113             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1114          } else {
1115             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1116          }
1117          break;
1118          case(DegreesOfFreedom):
1119          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1120          break;
1121          case(ReducedDegreesOfFreedom):
1122          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1123          break;
1124          case(Elements):
1125          case(ReducedElements):
1126          if (getMPISize()>1) {
1127             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1128             escriptDataC _in2 = temp.getDataC();
1129             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1130          } else {
1131             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1132          }
1133          break;
1134          case(FaceElements):
1135          case(ReducedFaceElements):
1136          if (getMPISize()>1) {
1137             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1138             escriptDataC _in2 = temp.getDataC();
1139             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1140          } else {
1141             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1142          }
1143          break;
1144          case(Points):
1145          if (getMPISize()>1) {
1146             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1147             escriptDataC _in2 = temp.getDataC();
1148             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1149          } else {
1150             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1151          }
1152          break;
1153          case(ContactElementsZero):
1154          case(ContactElementsOne):
1155          case(ReducedContactElementsZero):
1156          case(ReducedContactElementsOne):
1157          if (getMPISize()>1) {
1158             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1159             escriptDataC _in2 = temp.getDataC();
1160             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1161          } else {
1162             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1163          }
1164          break;
1165          default:
1166             stringstream temp;
1167             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1168             throw FinleyAdapterException(temp.str());
1169             break;
1170          }
1171        break;        break;
    }  
    break;  
1172     default:     default:
1173        stringstream temp;        stringstream temp;
1174        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
# Line 1165  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 1235  void MeshAdapter::setToNormal(escript::D Line 1262  void MeshAdapter::setToNormal(escript::D
1262  //  //
1263  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1264  {  {
1265     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1266     if (targetDomain!=*this)     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1267       if (targetDomain!=this)
1268        throw FinleyAdapterException("Error - Illegal domain of interpolation target");        throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1269    
1270     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
# Line 1245  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 1258  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 1295  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 1332  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 1444  void MeshAdapter::setToSize(escript::Dat Line 1472  void MeshAdapter::setToSize(escript::Dat
1472     checkFinleyError();     checkFinleyError();
1473  }  }
1474    
1475  // sets the location of nodes:  //
1476    // sets the location of nodes
1477    //
1478  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1479  {  {
1480     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1452  void MeshAdapter::setNewX(const escript: Line 1482  void MeshAdapter::setNewX(const escript:
1482     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1483     if (newDomain!=*this)     if (newDomain!=*this)
1484        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1485     tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1486     Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1487           Finley_Mesh_setCoordinates(mesh,&tmp);
1488       } else {
1489           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1490           tmp = new_x_inter.getDataC();
1491           Finley_Mesh_setCoordinates(mesh,&tmp);
1492       }
1493     checkFinleyError();     checkFinleyError();
1494  }  }
1495    
1496  // saves a data array in openDX format:  //
1497  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  // Helper for the save* methods. Extracts optional data variable names and the
1498    // corresponding pointers from python dictionary. Caller must free arrays.
1499    //
1500    void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1501  {  {
1502     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;     numData = boost::python::extract<int>(arg.attr("__len__")());
    const int num_data=boost::python::extract<int>(arg.attr("__len__")());  
1503     /* win32 refactor */     /* win32 refactor */
1504     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1505     for(int i=0;i<num_data;i++)     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1506     {     dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
       names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
    }  
   
    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  
    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  
1507    
1508     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1509     for (int i=0;i<num_data;++i) {     for (int i=0; i<numData; ++i) {
1510        std::string n=boost::python::extract<std::string>(keys[i]);        string n=boost::python::extract<string>(keys[i]);
1511        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1512        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1513           throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");           throw FinleyAdapterException("Error: Data must be defined on same Domain");
1514        data[i]=d.getDataC();        data[i] = d.getDataC();
1515        ptr_data[i]= &(data[i]);        dataPtr[i] = &(data[i]);
1516        if (n.length()>MAX_namelength-1) {        names[i] = TMPMEMALLOC(n.length()+1, char);
1517           strncpy(names[i],n.c_str(),MAX_namelength-1);        strcpy(names[i], n.c_str());
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
1518     }     }
1519     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  }
1520    
1521    //
1522    // saves mesh and optionally data arrays in openDX format
1523    //
1524    void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1525    {
1526       int num_data;
1527       char **names;
1528       escriptDataC *data;
1529       escriptDataC **ptr_data;
1530    
1531       extractArgsFromDict(arg, num_data, names, data, ptr_data);
1532       Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1533     checkFinleyError();     checkFinleyError();
1534    
1535     /* win32 refactor */     /* win32 refactor */
1536     TMPMEMFREE(data);     TMPMEMFREE(data);
1537     TMPMEMFREE(ptr_data);     TMPMEMFREE(ptr_data);
1538     for(int i=0;i<num_data;i++)     for(int i=0; i<num_data; i++)
1539     {     {
1540        TMPMEMFREE(names[i]);        TMPMEMFREE(names[i]);
1541     }     }
# Line 1502  void MeshAdapter::saveDX(const std::stri Line 1544  void MeshAdapter::saveDX(const std::stri
1544     return;     return;
1545  }  }
1546    
1547  // saves a data array in openVTK format:  //
1548  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  // saves mesh and optionally data arrays in VTK format
1549    //
1550    void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1551  {  {
1552     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("saveVTK");
1553     const int num_data=boost::python::extract<int>(arg.attr("__len__")());      pySaveVTK(*boost::python::make_tuple(filename,
1554     /* win32 refactor */                 const_cast<MeshAdapter*>(this)->getPtr(),
1555     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;                 metadata, metadata_schema), **arg);
1556     for(int i=0;i<num_data;i++)  }
1557     {  
1558        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1559     }  {
1560    #ifdef ESYS_MPI
1561     escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;      index_t myFirstNode=0, myLastNode=0, k=0;
1562     escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;      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    
    boost::python::list keys=arg.keys();  
    for (int i=0;i<num_data;++i) {  
       std::string n=boost::python::extract<std::string>(keys[i]);  
       escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);  
       if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)  
          throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");  
       data[i]=d.getDataC();  
       ptr_data[i]=&(data[i]);  
       if (n.length()>MAX_namelength-1) {  
          strncpy(names[i],n.c_str(),MAX_namelength-1);  
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
    }  
    Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  
1582    
    checkFinleyError();  
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0;i<num_data;i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
1583    
1584     return;  //
1585  }  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1586                                                                                                                                                                        //
1587                                                                                                                                                                        ASM_ptr MeshAdapter::newSystemMatrix(
 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  
 SystemMatrixAdapter 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 1592  SystemMatrixAdapter MeshAdapter::newSyst Line 1627  SystemMatrixAdapter MeshAdapter::newSyst
1627  #endif  #endif
1628     }     }
1629     else {     else {
1630        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1631     }     }
1632     checkPasoError();     checkPasoError();
1633     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1634     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     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  TransportProblemAdapter MeshAdapter::newTransportProblem(  //
1642                                                           const double theta,  ATP_ptr MeshAdapter::newTransportProblem(
1643                                                             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 1622  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 1665  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 1840  bool MeshAdapter::operator!=(const Abstr Line 2026  bool MeshAdapter::operator!=(const Abstr
2026    
2027  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2028  {  {
2029     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2030       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2031     checkPasoError();     checkPasoError();
2032     return out;     return out;
2033  }  }
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     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2037       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  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2058  {  {
2059     int *out = NULL;     int *out = NULL;
2060     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2025  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 2033  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 2043  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 2103  int MeshAdapter::getNumberOfTagsInUse(in Line 2291  int MeshAdapter::getNumberOfTagsInUse(in
2291    }    }
2292    return numTags;    return numTags;
2293  }  }
2294  int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  
2295    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2296  {  {
2297    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2298    index_t* tags=NULL;    index_t* tags=NULL;
# Line 2169  bool MeshAdapter::canTag(int functionSpa Line 2358  bool MeshAdapter::canTag(int functionSpa
2358    }    }
2359  }  }
2360    
2361    AbstractDomain::StatusType MeshAdapter::getStatus() const
2362    {
2363      Finley_Mesh* mesh=m_finleyMesh.get();
2364      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.1990  
changed lines
  Added in v.3344

  ViewVC Help
Powered by ViewVC 1.1.26