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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26