/[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 3522 by gross, Tue May 24 00:57:58 2011 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2010 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 18  Line 18 
18  #ifdef USE_NETCDF  #ifdef USE_NETCDF
19  #include <netcdfcpp.h>  #include <netcdfcpp.h>
20  #endif  #endif
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #include "paso/Paso_MPI.h"  
 #endif  
21  extern "C" {  extern "C" {
22  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
23  }  }
 #include <vector>  
24    
25  #define IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH  256  #include <boost/python/import.hpp>
26    
27  using namespace std;  using namespace std;
28  using namespace escript;  using namespace escript;
# Line 88  int MeshAdapter::getMPIRank() const Line 83  int MeshAdapter::getMPIRank() const
83  }  }
84  void MeshAdapter::MPIBarrier() const  void MeshAdapter::MPIBarrier() const
85  {  {
86  #ifdef PASO_MPI  #ifdef ESYS_MPI
87     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
88  #endif  #endif
89     return;     return;
# Line 99  bool MeshAdapter::onMasterProcessor() co Line 94  bool MeshAdapter::onMasterProcessor() co
94  }  }
95    
96    
97    #ifdef ESYS_MPI
98      MPI_Comm
99    #else
100      unsigned int
101    #endif
102    MeshAdapter::getMPIComm() const
103    {
104    #ifdef ESYS_MPI
105        return m_finleyMesh->MPIInfo->comm;
106    #else
107        return 0;
108    #endif
109    }
110    
111    
112  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
113     return m_finleyMesh.get();     return m_finleyMesh.get();
114  }  }
115    
116  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const string& fileName) const
117  {  {
118     char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;     char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
119     strcpy(fName,fileName.c_str());     strcpy(fName,fileName.c_str());
# Line 112  void MeshAdapter::write(const std::strin Line 122  void MeshAdapter::write(const std::strin
122     TMPMEMFREE(fName);     TMPMEMFREE(fName);
123  }  }
124    
125  void MeshAdapter::Print_Mesh_Info(const bool full=false) const  void MeshAdapter::Print_Mesh_Info(const bool full) const
126  {  {
127     Finley_PrintMesh_Info(m_finleyMesh.get(), full);     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
128  }  }
129    
130  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const string& fileName) const
131  {  {
132  #ifdef USE_NETCDF  #ifdef USE_NETCDF
133     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
# Line 137  void MeshAdapter::dump(const std::string Line 147  void MeshAdapter::dump(const std::string
147     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
148     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
149     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
150  #ifdef PASO_MPI  #ifdef ESYS_MPI
151     MPI_Status status;     MPI_Status status;
152  #endif  #endif
153    
154  /* Incoming token indicates it's my turn to write */  /* Incoming token indicates it's my turn to write */
155  #ifdef PASO_MPI  #ifdef ESYS_MPI
156     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
157  #endif  #endif
158    
159     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
160                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
161    
162     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
163     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
164     num_Tags = 0;     num_Tags = 0;
165     if (tag_map) {     while (tag_map) {
166        while (tag_map) {        num_Tags++;
167           num_Tags++;        tag_map=tag_map->next;
          tag_map=tag_map->next;  
       }  
168     }     }
169    
170     // NetCDF error handler     // NetCDF error handler
171     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
172     // Create the file.     // Create the file.
173     NcFile dataFile(newFileName, NcFile::Replace);     NcFile dataFile(newFileName, NcFile::Replace);
174       string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
175     // check if writing was successful     // check if writing was successful
176     if (!dataFile.is_valid())     if (!dataFile.is_valid())
177        throw DataException("Error - MeshAdapter::dump: opening of NetCDF file for output failed: " + *newFileName);        throw DataException(msgPrefix+"Open file for output");
178    
179     // Define dimensions (num_Elements and dim_Elements are identical, dim_Elements only appears if > 0)     // Define dimensions (num_Elements and dim_Elements are identical,
180       // dim_Elements only appears if > 0)
181     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
182        throw DataException("Error - MeshAdapter::dump: appending dimension numNodes to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numNodes)");
183     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
184        throw DataException("Error - MeshAdapter::dump: appending dimension numDim to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numDim)");
185     if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )     if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )
186        throw DataException("Error - MeshAdapter::dump: appending dimension mpi_size to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(mpi_size)");
187     if (num_Elements>0)     if (num_Elements>0)
188        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
189           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements)");
190     if (num_FaceElements>0)     if (num_FaceElements>0)
191        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
192           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
193     if (num_ContactElements>0)     if (num_ContactElements>0)
194        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
195           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements)");
196     if (num_Points>0)     if (num_Points>0)
197        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
198           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Points to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Points)");
199     if (num_Elements>0)     if (num_Elements>0)
200        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
201           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
202     if (num_FaceElements>0)     if (num_FaceElements>0)
203        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
204           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
205     if (num_ContactElements>0)     if (num_ContactElements>0)
206        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
207           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements_numNodes)");
208     if (num_Tags>0)     if (num_Tags>0)
209        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
210           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Tags)");
211    
212     // Attributes: MPI size, MPI rank, Name, order, reduced_order     // Attributes: MPI size, MPI rank, Name, order, reduced_order
213     if (!dataFile.add_att("mpi_size", mpi_size) )     if (!dataFile.add_att("mpi_size", mpi_size) )
214        throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_size)");
215     if (!dataFile.add_att("mpi_rank", mpi_rank) )     if (!dataFile.add_att("mpi_rank", mpi_rank) )
216        throw DataException("Error - MeshAdapter::dump: appending mpi_rank to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_rank)");
217     if (!dataFile.add_att("Name",mesh->Name) )     if (!dataFile.add_att("Name",mesh->Name) )
218        throw DataException("Error - MeshAdapter::dump: appending Name to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Name)");
219     if (!dataFile.add_att("numDim",numDim) )     if (!dataFile.add_att("numDim",numDim) )
220        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
221     if (!dataFile.add_att("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
222        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
223     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
224        throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(reduced_order)");
225     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
226        throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(numNodes)");
227     if (!dataFile.add_att("num_Elements",num_Elements) )     if (!dataFile.add_att("num_Elements",num_Elements) )
228        throw DataException("Error - MeshAdapter::dump: appending num_Elements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements)");
229     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
230        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements)");
231     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
232        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements)");
233     if (!dataFile.add_att("num_Points",num_Points) )     if (!dataFile.add_att("num_Points",num_Points) )
234        throw DataException("Error - MeshAdapter::dump: appending num_Points to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Points)");
235     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
236        throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
237     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
238        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
239     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
240        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements_numNodes)");
241     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
242        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Elements_TypeId)");
243     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
244        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
245     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
246        throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(ContactElements_TypeId)");
247     if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
248        throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Points_TypeId)");
249     if (!dataFile.add_att("num_Tags", num_Tags) )     if (!dataFile.add_att("num_Tags", num_Tags) )
250        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Tags)");
251    
252     // // // // // Nodes // // // // //     // // // // // Nodes // // // // //
253    
254     // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)     // Nodes nodeDistribution
255       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
256          throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
257       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
258       if (! (ids->put(int_ptr, mpi_size+1)) )
259          throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
260    
261       // Nodes degreesOfFreedomDistribution
262       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
263          throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
264       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
265       if (! (ids->put(int_ptr, mpi_size+1)) )
266          throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
267    
268       // Only write nodes if non-empty because NetCDF doesn't like empty arrays
269       // (it treats them as NC_UNLIMITED)
270     if (numNodes>0) {     if (numNodes>0) {
271    
272        // Nodes Id        // Nodes Id
273        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
274           throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Id)");
275        int_ptr = &mesh->Nodes->Id[0];        int_ptr = &mesh->Nodes->Id[0];
276        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
277           throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Id)");
278    
279        // Nodes Tag        // Nodes Tag
280        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
281           throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Tag)");
282        int_ptr = &mesh->Nodes->Tag[0];        int_ptr = &mesh->Nodes->Tag[0];
283        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
284           throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Tag)");
285    
286        // Nodes gDOF        // Nodes gDOF
287        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
288           throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
289        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
290        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
291           throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gDOF)");
292    
293        // Nodes global node index        // Nodes global node index
294        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
295           throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gNI)");
296        int_ptr = &mesh->Nodes->globalNodesIndex[0];        int_ptr = &mesh->Nodes->globalNodesIndex[0];
297        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
298           throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gNI)");
299    
300        // Nodes grDof        // Nodes grDof
301        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
302           throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
303        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
304        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
305           throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grDfI)");
306    
307        // Nodes grNI        // Nodes grNI
308        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
309           throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grNI)");
310        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
311        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
312           throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grNI)");
313    
314        // Nodes Coordinates        // Nodes Coordinates
315        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
316           throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
317        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
318           throw DataException("Error - MeshAdapter::dump: copy Nodes_Coordinates to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Coordinates)");
   
       // Nodes degreesOfFreedomDistribution  
       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )  
          throw DataException("Error - MeshAdapter::dump: appending Nodes_DofDistribution to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];  
       if (! (ids->put(int_ptr, mpi_size+1)) )  
          throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);  
   
       // Nodes nodeDistribution  
       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )  
          throw DataException("Error - MeshAdapter::dump: appending Nodes_NodeDistribution to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];  
       if (! (ids->put(int_ptr, mpi_size+1)) )  
          throw DataException("Error - MeshAdapter::dump: copy Nodes_NodeDistribution to netCDF buffer failed: " + *newFileName);  
319    
320     }     }
321    
# Line 314  void MeshAdapter::dump(const std::string Line 325  void MeshAdapter::dump(const std::string
325    
326        // Elements_Id        // Elements_Id
327        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
328           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Id)");
329        int_ptr = &mesh->Elements->Id[0];        int_ptr = &mesh->Elements->Id[0];
330        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
331           throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Id)");
332    
333        // Elements_Tag        // Elements_Tag
334        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
335           throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Tag)");
336        int_ptr = &mesh->Elements->Tag[0];        int_ptr = &mesh->Elements->Tag[0];
337        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
338           throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Tag)");
339    
340        // Elements_Owner        // Elements_Owner
341        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
342           throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Owner)");
343        int_ptr = &mesh->Elements->Owner[0];        int_ptr = &mesh->Elements->Owner[0];
344        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
345           throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Owner)");
346    
347        // Elements_Color        // Elements_Color
348        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
349           throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Color)");
350        int_ptr = &mesh->Elements->Color[0];        int_ptr = &mesh->Elements->Color[0];
351        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
352           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Color)");
353    
354        // Elements_Nodes        // Elements_Nodes
355        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
356           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Nodes)");
357        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
358           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Nodes)");
359    
360     }     }
361    
# Line 354  void MeshAdapter::dump(const std::string Line 365  void MeshAdapter::dump(const std::string
365    
366        // FaceElements_Id        // FaceElements_Id
367        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
368           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Id)");
369        int_ptr = &mesh->FaceElements->Id[0];        int_ptr = &mesh->FaceElements->Id[0];
370        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
371           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Id)");
372    
373        // FaceElements_Tag        // FaceElements_Tag
374        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
375           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
376        int_ptr = &mesh->FaceElements->Tag[0];        int_ptr = &mesh->FaceElements->Tag[0];
377        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
378           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Tag)");
379    
380        // FaceElements_Owner        // FaceElements_Owner
381        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
382           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
383        int_ptr = &mesh->FaceElements->Owner[0];        int_ptr = &mesh->FaceElements->Owner[0];
384        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
385           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Owner)");
386    
387        // FaceElements_Color        // FaceElements_Color
388        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
389           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Color)");
390        int_ptr = &mesh->FaceElements->Color[0];        int_ptr = &mesh->FaceElements->Color[0];
391        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
392           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Color)");
393    
394        // FaceElements_Nodes        // FaceElements_Nodes
395        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
396           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
397        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
398           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Nodes)");
399    
400     }     }
401    
# Line 394  void MeshAdapter::dump(const std::string Line 405  void MeshAdapter::dump(const std::string
405    
406        // ContactElements_Id        // ContactElements_Id
407        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
408           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Id)");
409        int_ptr = &mesh->ContactElements->Id[0];        int_ptr = &mesh->ContactElements->Id[0];
410        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
411           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Id)");
412    
413        // ContactElements_Tag        // ContactElements_Tag
414        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
415           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Tag)");
416        int_ptr = &mesh->ContactElements->Tag[0];        int_ptr = &mesh->ContactElements->Tag[0];
417        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
418           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Tag)");
419    
420        // ContactElements_Owner        // ContactElements_Owner
421        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
422           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Owner)");
423        int_ptr = &mesh->ContactElements->Owner[0];        int_ptr = &mesh->ContactElements->Owner[0];
424        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
425           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Owner)");
426    
427        // ContactElements_Color        // ContactElements_Color
428        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
429           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Color)");
430        int_ptr = &mesh->ContactElements->Color[0];        int_ptr = &mesh->ContactElements->Color[0];
431        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
432           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Color)");
433    
434        // ContactElements_Nodes        // ContactElements_Nodes
435        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
436           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Nodes)");
437        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
438           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Nodes)");
439    
440     }     }
441    
# Line 436  void MeshAdapter::dump(const std::string Line 447  void MeshAdapter::dump(const std::string
447    
448        // Points_Id        // Points_Id
449        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
450           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Id)");
451        int_ptr = &mesh->Points->Id[0];        int_ptr = &mesh->Points->Id[0];
452        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
453           throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Id)");
454    
455        // Points_Tag        // Points_Tag
456        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
457           throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Tag)");
458        int_ptr = &mesh->Points->Tag[0];        int_ptr = &mesh->Points->Tag[0];
459        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
460           throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Tag)");
461    
462        // Points_Owner        // Points_Owner
463        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
464           throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Owner)");
465        int_ptr = &mesh->Points->Owner[0];        int_ptr = &mesh->Points->Owner[0];
466        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
467           throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Owner)");
468    
469        // Points_Color        // Points_Color
470        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
471           throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Color)");
472        int_ptr = &mesh->Points->Color[0];        int_ptr = &mesh->Points->Color[0];
473        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
474           throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Color)");
475    
476        // Points_Nodes        // Points_Nodes
477        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
478        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
479           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Nodes)");
480        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
481           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Nodes)");
482    
483     }     }
484    
# Line 491  void MeshAdapter::dump(const std::string Line 502  void MeshAdapter::dump(const std::string
502    
503        // Tags_keys        // Tags_keys
504        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
505           throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Tags_keys)");
506        int_ptr = &Tags_keys[0];        int_ptr = &Tags_keys[0];
507        if (! (ids->put(int_ptr, num_Tags)) )        if (! (ids->put(int_ptr, num_Tags)) )
508           throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Tags_keys)");
509    
510        // Tags_names_*        // Tags_names_*
511        // This is an array of strings, it should be stored as an array but instead I have hacked in one attribute per string        // This is an array of strings, it should be stored as an array but
512        // because the NetCDF manual doesn't tell how to do an array of strings        // instead I have hacked in one attribute per string because the NetCDF
513          // manual doesn't tell how to do an array of strings
514        tag_map = mesh->TagMap;        tag_map = mesh->TagMap;
515        if (tag_map) {        if (tag_map) {
516           int i = 0;           int i = 0;
517           while (tag_map) {           while (tag_map) {
518              sprintf(name_temp, "Tags_name_%d", i);              sprintf(name_temp, "Tags_name_%d", i);
519              if (!dataFile.add_att(name_temp, tag_map->name) )              if (!dataFile.add_att(name_temp, tag_map->name) )
520                 throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);                 throw DataException(msgPrefix+"add_att(Tags_names_XX)");
521              tag_map=tag_map->next;              tag_map=tag_map->next;
522              i++;              i++;
523           }           }
524        }        }
525    
526        TMPMEMFREE(Tags_keys);        TMPMEMFREE(Tags_keys);
   
527     }     }
528    
529  /* Send token to next MPI process so he can take his turn */  /* Send token to next MPI process so he can take his turn */
530  #ifdef PASO_MPI  #ifdef ESYS_MPI
531     if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);     if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
532  #endif  #endif
533    
# Line 636  int MeshAdapter::getReducedSolutionCode( Line 647  int MeshAdapter::getReducedSolutionCode(
647     return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
648  }  }
649    
650  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
651  {  {
652     return Points;     return Points;
653  }  }
# Line 681  pair<int,int> MeshAdapter::getDataShape( Line 692  pair<int,int> MeshAdapter::getDataShape(
692     case(Elements):     case(Elements):
693     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
694        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
695        numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
696     }     }
697     break;     break;
698     case(ReducedElements):     case(ReducedElements):
699     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
700        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
701        numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
702     }     }
703     break;     break;
704     case(FaceElements):     case(FaceElements):
705     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
706        numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
707        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
708     }     }
709     break;     break;
710     case(ReducedFaceElements):     case(ReducedFaceElements):
711     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
712        numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
713        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
714     }     }
715     break;     break;
# Line 710  pair<int,int> MeshAdapter::getDataShape( Line 721  pair<int,int> MeshAdapter::getDataShape(
721     break;     break;
722     case(ContactElementsZero):     case(ContactElementsZero):
723     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
724        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
725        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
726     }     }
727     break;     break;
728     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
729     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
730        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
731        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
732     }     }
733     break;     break;
734     case(ContactElementsOne):     case(ContactElementsOne):
735     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
736        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
737        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
738     }     }
739     break;     break;
740     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
741     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
742        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
743        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
744     }     }
745     break;     break;
# Line 757  pair<int,int> MeshAdapter::getDataShape( Line 768  pair<int,int> MeshAdapter::getDataShape(
768  // adds linear PDE of second order into a given stiffness matrix and right hand side:  // adds linear PDE of second order into a given stiffness matrix and right hand side:
769  //  //
770  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
771                                   SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
772                                   const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                   const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
773                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
774                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact,
775                                     const escript::Data& d_dirac,const escript::Data& y_dirac) const
776  {  {
777       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
778       if (smat==0)
779       {
780        throw FinleyAdapterException("finley only supports its own system matrices.");
781       }
782     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
783     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
784     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 773  void MeshAdapter::addPDEToSystem( Line 790  void MeshAdapter::addPDEToSystem(
790     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
791     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
792     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
793       escriptDataC _d_dirac=d_dirac.getDataC();
794       escriptDataC _y_dirac=y_dirac.getDataC();
795    
796     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
797    
798     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 );
799       checkFinleyError();
800    
801       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
802     checkFinleyError();     checkFinleyError();
803    
804     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
805     checkFinleyError();     checkFinleyError();
806    
807     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->Points, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
808     checkFinleyError();     checkFinleyError();
809  }  }
810    
811  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
812                                          escript::Data& mat,                                          escript::Data& mat,
813                                          const escript::Data& D,                                          const escript::Data& D,
814                                          const escript::Data& d) const                                          const escript::Data& d,
815                                            const escript::Data& d_dirac,
816                                            const bool useHRZ) const
817  {  {
818     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
819     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
820     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
821       escriptDataC _d_dirac=d_dirac.getDataC();
822    
823     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
824    
825     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
826     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     checkFinleyError();
827      
828       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
829       checkFinleyError();
830    
831       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Points,&_mat, &_d_dirac, useHRZ);
832     checkFinleyError();     checkFinleyError();
833    
834  }  }
835    
836    
837  //  //
838  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
839  //  //
840  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact,  const escript::Data& y_dirac) const
841  {  {
842     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
843    
# Line 816  void MeshAdapter::addPDEToRHS( escript:: Line 846  void MeshAdapter::addPDEToRHS( escript::
846     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
847     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
848     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
849       escriptDataC _y_dirac=y_dirac.getDataC();
850    
851     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
852     checkFinleyError();     checkFinleyError();
# Line 825  void MeshAdapter::addPDEToRHS( escript:: Line 856  void MeshAdapter::addPDEToRHS( escript::
856    
857     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
858     checkFinleyError();     checkFinleyError();
859    
860       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs , 0, 0, 0, 0, 0, &_y_dirac );
861       checkFinleyError();
862    
863  }  }
864  //  //
865  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
866  //  //
867  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
868                                             TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
869                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
870                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
871                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
872                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact,
873                                               const escript::Data& d_dirac, const escript::Data& y_dirac) const
874  {  {
875       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
876       if (tpa==0)
877       {
878        throw FinleyAdapterException("finley only supports its own transport problems.");
879       }
880    
881    
882     DataTypes::ShapeType shape;     DataTypes::ShapeType shape;
883     source.expand();     source.expand();
884     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
# Line 850  void MeshAdapter::addPDEToTransportProbl Line 893  void MeshAdapter::addPDEToTransportProbl
893     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
894     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
895     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
896       escriptDataC _d_dirac=d_dirac.getDataC();
897       escriptDataC _y_dirac=y_dirac.getDataC();
898    
899    
900     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
901     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
902    
903     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 );
904     checkFinleyError();     checkFinleyError();
# Line 865  void MeshAdapter::addPDEToTransportProbl Line 911  void MeshAdapter::addPDEToTransportProbl
911    
912     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
913     checkFinleyError();     checkFinleyError();
914    
915       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
916       checkFinleyError();
917    
918  }  }
919    
920  //  //
# Line 884  void MeshAdapter::interpolateOnDomain(es Line 934  void MeshAdapter::interpolateOnDomain(es
934     escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
935     switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
936     case(Nodes):     case(Nodes):
937     switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
938     case(Nodes):        case(Nodes):
939     case(ReducedNodes):        case(ReducedNodes):
940     case(DegreesOfFreedom):        case(DegreesOfFreedom):
941     case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
942     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());  
943        break;        break;
944     }        case(Elements):
945     break;        case(ReducedElements):
946     case(ReducedNodes):        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
947     switch(target.getFunctionSpace().getTypeCode()) {        break;
948     case(Nodes):        case(FaceElements):
949     case(ReducedNodes):        case(ReducedFaceElements):
950     case(DegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
951     case(ReducedDegreesOfFreedom):        break;
952     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(Points):
953     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
954     case(Elements):        break;
955     case(ReducedElements):        case(ContactElementsZero):
956     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        case(ReducedContactElementsZero):
957     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
958     case(FaceElements):        break;
959     case(ReducedFaceElements):        case(ContactElementsOne):
960     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        case(ReducedContactElementsOne):
961     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
962     case(Points):        break;
963     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        default:
964     break;           stringstream temp;
965     case(ContactElementsZero):           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
966     case(ReducedContactElementsZero):           throw FinleyAdapterException(temp.str());
967     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           break;
968     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());  
969        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):  
970     case(ReducedNodes):     case(ReducedNodes):
971     if (getMPISize()>1) {        switch(target.getFunctionSpace().getTypeCode()) {
972        escript::Data temp=escript::Data(in);        case(Nodes):
973        temp.expand();        case(ReducedNodes):
974        escriptDataC _in2 = temp.getDataC();        case(DegreesOfFreedom):
975        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);        case(ReducedDegreesOfFreedom):
    } else {  
976        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
977     }        break;
978     break;        case(Elements):
979     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 {  
980        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
981     }        break;
982     break;        case(FaceElements):
983     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 {  
984        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
985     }        break;
986     break;        case(Points):
    case(Points):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
    } else {  
987        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
988     }        break;
989     break;        case(ContactElementsZero):
990     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 {  
991        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());  
992        break;        break;
993     }        case(ContactElementsOne):
994     break;        case(ReducedContactElementsOne):
995     case(ReducedDegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
996     switch(target.getFunctionSpace().getTypeCode()) {        break;
997     case(Nodes):        default:
998     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           stringstream temp;
999     break;           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1000     case(ReducedNodes):           throw FinleyAdapterException(temp.str());
1001     if (getMPISize()>1) {           break;
1002        escript::Data temp=escript::Data(in);        }
1003        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;  
1004     case(Elements):     case(Elements):
1005          if (target.getFunctionSpace().getTypeCode()==Elements) {
1006             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1007          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1008             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
1009          } else {
1010             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
1011          }
1012          break;
1013     case(ReducedElements):     case(ReducedElements):
1014     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1015        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1016        escriptDataC _in2 = temp.getDataC();        } else {
1017        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
1018     } else {        }
1019        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        break;
    }  
    break;  
1020     case(FaceElements):     case(FaceElements):
1021          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
1022             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1023          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1024             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1025          } else {
1026             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1027          }
1028          break;
1029     case(ReducedFaceElements):     case(ReducedFaceElements):
1030     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1031        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1032        escriptDataC _in2 = temp.getDataC();        } else {
1033        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
1034     } else {        }
1035        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
    }  
    break;  
1036     case(Points):     case(Points):
1037     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==Points) {
1038        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1039        escriptDataC _in2 = temp.getDataC();        } else {
1040        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1041     } else {        }
1042        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        break;
    }  
    break;  
1043     case(ContactElementsZero):     case(ContactElementsZero):
1044     case(ContactElementsOne):     case(ContactElementsOne):
1045          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1046             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1047          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1048             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1049          } else {
1050             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1051          }
1052          break;
1053     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1054     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1055     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1056        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1057        escriptDataC _in2 = temp.getDataC();        } else {
1058        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1059     } else {        }
1060        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        break;
1061     }     case(DegreesOfFreedom):      
1062     break;        switch(target.getFunctionSpace().getTypeCode()) {
1063     default:        case(ReducedDegreesOfFreedom):
1064        stringstream temp;        case(DegreesOfFreedom):
1065        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1066        throw FinleyAdapterException(temp.str());        break;
1067      
1068          case(Nodes):
1069          case(ReducedNodes):
1070          if (getMPISize()>1) {
1071             escript::Data temp=escript::Data(in);
1072             temp.expand();
1073             escriptDataC _in2 = temp.getDataC();
1074             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1075          } else {
1076             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1077          }
1078          break;
1079          case(Elements):
1080          case(ReducedElements):
1081          if (getMPISize()>1) {
1082             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1083             escriptDataC _in2 = temp.getDataC();
1084             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1085          } else {
1086             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1087          }
1088          break;
1089          case(FaceElements):
1090          case(ReducedFaceElements):
1091          if (getMPISize()>1) {
1092             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1093             escriptDataC _in2 = temp.getDataC();
1094             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1095      
1096          } else {
1097             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1098          }
1099          break;
1100          case(Points):
1101          if (getMPISize()>1) {
1102             //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1103             //escriptDataC _in2 = temp.getDataC();
1104          } else {
1105             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1106          }
1107          break;
1108          case(ContactElementsZero):
1109          case(ContactElementsOne):
1110          case(ReducedContactElementsZero):
1111          case(ReducedContactElementsOne):
1112          if (getMPISize()>1) {
1113             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1114             escriptDataC _in2 = temp.getDataC();
1115             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1116          } else {
1117             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1118          }
1119          break;
1120          default:
1121             stringstream temp;
1122             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1123             throw FinleyAdapterException(temp.str());
1124             break;
1125          }
1126          break;
1127       case(ReducedDegreesOfFreedom):
1128          switch(target.getFunctionSpace().getTypeCode()) {
1129          case(Nodes):
1130          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1131          break;
1132          case(ReducedNodes):
1133          if (getMPISize()>1) {
1134             escript::Data temp=escript::Data(in);
1135             temp.expand();
1136             escriptDataC _in2 = temp.getDataC();
1137             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1138          } else {
1139             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1140          }
1141          break;
1142          case(DegreesOfFreedom):
1143          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1144          break;
1145          case(ReducedDegreesOfFreedom):
1146          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1147          break;
1148          case(Elements):
1149          case(ReducedElements):
1150          if (getMPISize()>1) {
1151             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1152             escriptDataC _in2 = temp.getDataC();
1153             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1154          } else {
1155             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1156          }
1157          break;
1158          case(FaceElements):
1159          case(ReducedFaceElements):
1160          if (getMPISize()>1) {
1161             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1162             escriptDataC _in2 = temp.getDataC();
1163             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1164          } else {
1165             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1166          }
1167          break;
1168          case(Points):
1169          if (getMPISize()>1) {
1170             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1171             escriptDataC _in2 = temp.getDataC();
1172             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1173          } else {
1174             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1175          }
1176          break;
1177          case(ContactElementsZero):
1178          case(ContactElementsOne):
1179          case(ReducedContactElementsZero):
1180          case(ReducedContactElementsOne):
1181          if (getMPISize()>1) {
1182             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1183             escriptDataC _in2 = temp.getDataC();
1184             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1185          } else {
1186             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1187          }
1188          break;
1189          default:
1190             stringstream temp;
1191             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1192             throw FinleyAdapterException(temp.str());
1193             break;
1194          }
1195        break;        break;
    }  
    break;  
1196     default:     default:
1197        stringstream temp;        stringstream temp;
1198        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 1216  void MeshAdapter::setToX(escript::Data&
1216        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1217        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1218     } else {     } else {
1219        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1220        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1221        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1222        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1247  void MeshAdapter::interpolateACross(escr Line 1297  void MeshAdapter::interpolateACross(escr
1297  //  //
1298  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
1299  //  //
1300  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1301  {  {
1302     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1303     if (argDomain!=*this)     if (argDomain!=*this)
# Line 1260  void MeshAdapter::setToIntegrals(std::ve Line 1310  void MeshAdapter::setToIntegrals(std::ve
1310     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1311     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1312     case(Nodes):     case(Nodes):
1313     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1314     _temp=temp.getDataC();     _temp=temp.getDataC();
1315     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1316     break;     break;
1317     case(ReducedNodes):     case(ReducedNodes):
1318     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1319     _temp=temp.getDataC();     _temp=temp.getDataC();
1320     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1321     break;     break;
# Line 1297  void MeshAdapter::setToIntegrals(std::ve Line 1347  void MeshAdapter::setToIntegrals(std::ve
1347     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1348     break;     break;
1349     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1350     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1351     _temp=temp.getDataC();     _temp=temp.getDataC();
1352     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1353     break;     break;
1354     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1355     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1356     _temp=temp.getDataC();     _temp=temp.getDataC();
1357     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1358     break;     break;
# Line 1334  void MeshAdapter::setToGradient(escript: Line 1384  void MeshAdapter::setToGradient(escript:
1384     escript::Data temp;     escript::Data temp;
1385     if (getMPISize()>1) {     if (getMPISize()>1) {
1386        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1387           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1388           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1389        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1390           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1391           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1392        } else {        } else {
1393           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1446  void MeshAdapter::setToSize(escript::Dat Line 1496  void MeshAdapter::setToSize(escript::Dat
1496     checkFinleyError();     checkFinleyError();
1497  }  }
1498    
1499  // sets the location of nodes:  //
1500    // sets the location of nodes
1501    //
1502  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1503  {  {
1504     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1454  void MeshAdapter::setNewX(const escript: Line 1506  void MeshAdapter::setNewX(const escript:
1506     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1507     if (newDomain!=*this)     if (newDomain!=*this)
1508        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1509     tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1510     Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1511           Finley_Mesh_setCoordinates(mesh,&tmp);
1512       } else {
1513           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1514           tmp = new_x_inter.getDataC();
1515           Finley_Mesh_setCoordinates(mesh,&tmp);
1516       }
1517     checkFinleyError();     checkFinleyError();
1518  }  }
1519    
1520  // saves a data array in openDX format:  //
1521  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
1522    // corresponding pointers from python dictionary. Caller must free arrays.
1523    //
1524    void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1525  {  {
1526     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__")());  
1527     /* win32 refactor */     /* win32 refactor */
1528     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1529     for(int i=0;i<num_data;i++)     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1530     {     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;  
1531    
1532     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1533     for (int i=0;i<num_data;++i) {     for (int i=0; i<numData; ++i) {
1534        std::string n=boost::python::extract<std::string>(keys[i]);        string n=boost::python::extract<string>(keys[i]);
1535        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1536        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1537           throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");           throw FinleyAdapterException("Error: Data must be defined on same Domain");
1538        data[i]=d.getDataC();        data[i] = d.getDataC();
1539        ptr_data[i]= &(data[i]);        dataPtr[i] = &(data[i]);
1540        if (n.length()>MAX_namelength-1) {        names[i] = TMPMEMALLOC(n.length()+1, char);
1541           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());  
       }  
1542     }     }
1543     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  }
1544    
1545    //
1546    // saves mesh and optionally data arrays in openDX format
1547    //
1548    void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1549    {
1550       int num_data;
1551       char **names;
1552       escriptDataC *data;
1553       escriptDataC **ptr_data;
1554    
1555       extractArgsFromDict(arg, num_data, names, data, ptr_data);
1556       Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1557     checkFinleyError();     checkFinleyError();
1558    
1559     /* win32 refactor */     /* win32 refactor */
1560     TMPMEMFREE(data);     TMPMEMFREE(data);
1561     TMPMEMFREE(ptr_data);     TMPMEMFREE(ptr_data);
1562     for(int i=0;i<num_data;i++)     for(int i=0; i<num_data; i++)
1563     {     {
1564        TMPMEMFREE(names[i]);        TMPMEMFREE(names[i]);
1565     }     }
# Line 1504  void MeshAdapter::saveDX(const std::stri Line 1568  void MeshAdapter::saveDX(const std::stri
1568     return;     return;
1569  }  }
1570    
1571  // saves a data array in openVTK format:  //
1572  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  // saves mesh and optionally data arrays in VTK format
1573    //
1574    void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1575  {  {
1576     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1577     const int num_data=boost::python::extract<int>(arg.attr("__len__")());      pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1578     /* win32 refactor */                metadata, metadata_schema, arg);
1579     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;  }
1580     for(int i=0;i<num_data;i++)  
1581     {  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1582        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  {
1583     }  #ifdef ESYS_MPI
1584        index_t myFirstNode=0, myLastNode=0, k=0;
1585        index_t* globalNodeIndex=0;
1586        Finley_Mesh* mesh_p=m_finleyMesh.get();
1587        if (fs_code == FINLEY_REDUCED_NODES)
1588        {
1589        myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1590        myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1591        globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1592        }
1593        else
1594        {
1595        myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1596        myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1597        globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1598        }
1599        k=globalNodeIndex[id];
1600        return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1601    #endif
1602        return true;
1603    }
1604    
    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  
    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  
1605    
    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);  
   
    checkFinleyError();  
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0;i<num_data;i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
1606    
1607     return;  //
1608  }  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1609                                                                                                                                                                        //
1610                                                                                                                                                                        ASM_ptr MeshAdapter::newSystemMatrix(
 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  
 SystemMatrixAdapter MeshAdapter::newSystemMatrix(  
1611                                                   const int row_blocksize,                                                   const int row_blocksize,
1612                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1613                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1594  SystemMatrixAdapter MeshAdapter::newSyst Line 1650  SystemMatrixAdapter MeshAdapter::newSyst
1650  #endif  #endif
1651     }     }
1652     else {     else {
1653        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1654     }     }
1655     checkPasoError();     checkPasoError();
1656     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1657     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1658       return ASM_ptr(sma);
1659    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1660  }  }
1661    
1662    //
1663  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1664  TransportProblemAdapter MeshAdapter::newTransportProblem(  //
1665                                                           const double theta,  ATP_ptr MeshAdapter::newTransportProblem(
1666                                                             const bool useBackwardEuler,
1667                                                           const int blocksize,                                                           const int blocksize,
1668                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1669                                                           const int type) const                                                           const int type) const
# Line 1624  TransportProblemAdapter MeshAdapter::new Line 1685  TransportProblemAdapter MeshAdapter::new
1685    
1686     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1687     checkFinleyError();     checkFinleyError();
1688     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1689     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1690     checkPasoError();     checkPasoError();
1691     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1692     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1693       return ATP_ptr(tpa);
1694    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1695  }  }
1696    
1697  //  //
# Line 1667  bool MeshAdapter::isCellOriented(int fun Line 1730  bool MeshAdapter::isCellOriented(int fun
1730     return false;     return false;
1731  }  }
1732    
1733    bool
1734    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1735    {
1736       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1737        class 1: DOF <-> Nodes
1738        class 2: ReducedDOF <-> ReducedNodes
1739        class 3: Points
1740        class 4: Elements
1741        class 5: ReducedElements
1742        class 6: FaceElements
1743        class 7: ReducedFaceElements
1744        class 8: ContactElementZero <-> ContactElementOne
1745        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1746    
1747       There is also a set of lines. Interpolation is possible down a line but not between lines.
1748       class 1 and 2 belong to all lines so aren't considered.
1749        line 0: class 3
1750        line 1: class 4,5
1751        line 2: class 6,7
1752        line 3: class 8,9
1753    
1754       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1755       eg hasnodes is true if we have at least one instance of Nodes.
1756       */
1757        if (fs.empty())
1758        {
1759            return false;
1760        }
1761        vector<int> hasclass(10);
1762        vector<int> hasline(4);
1763        bool hasnodes=false;
1764        bool hasrednodes=false;
1765        bool hascez=false;
1766        bool hasrcez=false;
1767        for (int i=0;i<fs.size();++i)
1768        {
1769        switch(fs[i])
1770        {
1771        case(Nodes):   hasnodes=true;   // no break is deliberate
1772        case(DegreesOfFreedom):
1773            hasclass[1]=1;
1774            break;
1775        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1776        case(ReducedDegreesOfFreedom):
1777            hasclass[2]=1;
1778            break;
1779        case(Points):
1780            hasline[0]=1;
1781            hasclass[3]=1;
1782            break;
1783        case(Elements):
1784            hasclass[4]=1;
1785            hasline[1]=1;
1786            break;
1787        case(ReducedElements):
1788            hasclass[5]=1;
1789            hasline[1]=1;
1790            break;
1791        case(FaceElements):
1792            hasclass[6]=1;
1793            hasline[2]=1;
1794            break;
1795        case(ReducedFaceElements):
1796            hasclass[7]=1;
1797            hasline[2]=1;
1798            break;
1799        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1800        case(ContactElementsOne):
1801            hasclass[8]=1;
1802            hasline[3]=1;
1803            break;
1804        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1805        case(ReducedContactElementsOne):
1806            hasclass[9]=1;
1807            hasline[3]=1;
1808            break;
1809        default:
1810            return false;
1811        }
1812        }
1813        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1814        // fail if we have more than one leaf group
1815    
1816        if (totlines>1)
1817        {
1818        return false;   // there are at least two branches we can't interpolate between
1819        }
1820        else if (totlines==1)
1821        {
1822        if (hasline[0]==1)      // we have points
1823        {
1824            resultcode=Points;
1825        }
1826        else if (hasline[1]==1)
1827        {
1828            if (hasclass[5]==1)
1829            {
1830            resultcode=ReducedElements;
1831            }
1832            else
1833            {
1834            resultcode=Elements;
1835            }
1836        }
1837        else if (hasline[2]==1)
1838        {
1839            if (hasclass[7]==1)
1840            {
1841            resultcode=ReducedFaceElements;
1842            }
1843            else
1844            {
1845            resultcode=FaceElements;
1846            }
1847        }
1848        else    // so we must be in line3
1849        {
1850            if (hasclass[9]==1)
1851            {
1852            // need something from class 9
1853            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1854            }
1855            else
1856            {
1857            // something from class 8
1858            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1859            }
1860        }
1861        }
1862        else    // totlines==0
1863        {
1864        if (hasclass[2]==1)
1865        {
1866            // something from class 2
1867            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1868        }
1869        else
1870        {   // something from class 1
1871            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1872        }
1873        }
1874        return true;
1875    }
1876    
1877  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1878  {  {
1879     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1880     case(Nodes):     case(Nodes):
1881     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1882     case(Nodes):      case(Nodes):
1883     case(ReducedNodes):      case(ReducedNodes):
1884     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1885     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1886     case(Elements):      case(Elements):
1887     case(ReducedElements):      case(ReducedElements):
1888     case(FaceElements):      case(FaceElements):
1889     case(ReducedFaceElements):      case(ReducedFaceElements):
1890     case(Points):      case(Points):
1891     case(ContactElementsZero):      case(ContactElementsZero):
1892     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1893     case(ContactElementsOne):      case(ContactElementsOne):
1894     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1895     return true;      return true;
1896     default:      default:
1897        stringstream temp;            stringstream temp;
1898        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1899        throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1900     }     }
1901     break;     break;
1902     case(ReducedNodes):     case(ReducedNodes):
1903     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1904     case(ReducedNodes):      case(ReducedNodes):
1905     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1906     case(Elements):      case(Elements):
1907     case(ReducedElements):      case(ReducedElements):
1908     case(FaceElements):      case(FaceElements):
1909     case(ReducedFaceElements):      case(ReducedFaceElements):
1910     case(Points):      case(Points):
1911     case(ContactElementsZero):      case(ContactElementsZero):
1912     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1913     case(ContactElementsOne):      case(ContactElementsOne):
1914     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1915     return true;      return true;
1916     case(Nodes):      case(Nodes):
1917     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1918     return false;      return false;
1919     default:      default:
1920        stringstream temp;          stringstream temp;
1921        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;
1922        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1923     }     }
1924     break;     break;
1925     case(Elements):     case(Elements):
1926     if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1927        return true;        return true;
1928     } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1929        return true;        return true;
1930     } else {          } else {
1931        return false;            return false;
1932     }          }
1933     case(ReducedElements):     case(ReducedElements):
1934     if (functionSpaceType_target==ReducedElements) {      if (functionSpaceType_target==ReducedElements) {
1935        return true;        return true;
1936     } else {      } else {
1937        return false;            return false;
1938     }      }
1939     case(FaceElements):     case(FaceElements):
1940     if (functionSpaceType_target==FaceElements) {      if (functionSpaceType_target==FaceElements) {
1941        return true;              return true;
1942     } else if (functionSpaceType_target==ReducedFaceElements) {      } else if (functionSpaceType_target==ReducedFaceElements) {
1943        return true;              return true;
1944     } else {      } else {
1945        return false;              return false;
1946     }      }
1947     case(ReducedFaceElements):     case(ReducedFaceElements):
1948     if (functionSpaceType_target==ReducedFaceElements) {      if (functionSpaceType_target==ReducedFaceElements) {
1949        return true;              return true;
1950     } else {      } else {
1951        return false;          return false;
1952     }      }
1953     case(Points):     case(Points):
1954     if (functionSpaceType_target==Points) {      if (functionSpaceType_target==Points) {
1955        return true;              return true;
1956     } else {      } else {
1957        return false;              return false;
1958     }      }
1959     case(ContactElementsZero):     case(ContactElementsZero):
1960     case(ContactElementsOne):     case(ContactElementsOne):
1961     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1962        return true;              return true;
1963     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1964        return true;              return true;
1965     } else {      } else {
1966        return false;              return false;
1967     }      }
1968     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1969     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1970     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1971        return true;              return true;
1972     } else {      } else {
1973        return false;              return false;
1974     }      }
1975     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1976     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1977     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1978     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1979     case(Nodes):      case(Nodes):
1980     case(ReducedNodes):      case(ReducedNodes):
1981     case(Elements):      case(Elements):
1982     case(ReducedElements):      case(ReducedElements):
1983     case(Points):      case(Points):
1984     case(FaceElements):      case(FaceElements):
1985     case(ReducedFaceElements):      case(ReducedFaceElements):
1986     case(ContactElementsZero):      case(ContactElementsZero):
1987     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1988     case(ContactElementsOne):      case(ContactElementsOne):
1989     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1990     return true;      return true;
1991     default:      default:
1992        stringstream temp;          stringstream temp;
1993        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1994        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1995     }      }
1996     break;      break;
1997     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1998     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
1999     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
2000     case(ReducedNodes):      case(ReducedNodes):
2001     case(Elements):      case(Elements):
2002     case(ReducedElements):      case(ReducedElements):
2003     case(FaceElements):      case(FaceElements):
2004     case(ReducedFaceElements):      case(ReducedFaceElements):
2005     case(Points):      case(Points):
2006     case(ContactElementsZero):      case(ContactElementsZero):
2007     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
2008     case(ContactElementsOne):      case(ContactElementsOne):
2009     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
2010     return true;      return true;
2011     case(Nodes):      case(Nodes):
2012     case(DegreesOfFreedom):      case(DegreesOfFreedom):
2013     return false;      return false;
2014     default:      default:
2015        stringstream temp;          stringstream temp;
2016        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;
2017        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
2018     }      }
2019     break;      break;
2020     default:     default:
2021        stringstream temp;        stringstream temp;
2022        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 2049  bool MeshAdapter::operator!=(const Abstr
2049    
2050  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
2051  {  {
2052     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2053       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2054     checkPasoError();     checkPasoError();
2055     return out;     return out;
2056  }  }
2057  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
2058  {  {
2059     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2060       int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2061     checkPasoError();     checkPasoError();
2062     return out;     return out;
2063  }  }
2064    
2065  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2066  {  {
2067     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2068  }  }
2069    
2070  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2071  {  {
2072     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2073  }  }
2074    
2075  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2076  {  {
2077     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2078  }  }
2079    
2080  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2081  {  {
2082     int *out = NULL;     int *out = NULL;
2083     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2027  void MeshAdapter::setTags(const int func Line 2236  void MeshAdapter::setTags(const int func
2236     return;     return;
2237  }  }
2238    
2239  void MeshAdapter::setTagMap(const std::string& name,  int tag)  void MeshAdapter::setTagMap(const string& name,  int tag)
2240  {  {
2241     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2242     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 2244  void MeshAdapter::setTagMap(const std::s
2244     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2245  }  }
2246    
2247  int MeshAdapter::getTag(const std::string& name) const  int MeshAdapter::getTag(const string& name) const
2248  {  {
2249     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2250     int tag=0;     int tag=0;
# Line 2045  int MeshAdapter::getTag(const std::strin Line 2254  int MeshAdapter::getTag(const std::strin
2254     return tag;     return tag;
2255  }  }
2256    
2257  bool MeshAdapter::isValidTagName(const std::string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2258  {  {
2259     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2260     return Finley_Mesh_isValidTagName(mesh,name.c_str());     return Finley_Mesh_isValidTagName(mesh,name.c_str());
2261  }  }
2262    
2263  std::string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2264  {  {
2265     stringstream temp;     stringstream temp;
2266     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2105  int MeshAdapter::getNumberOfTagsInUse(in Line 2314  int MeshAdapter::getNumberOfTagsInUse(in
2314    }    }
2315    return numTags;    return numTags;
2316  }  }
2317  int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  
2318    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2319  {  {
2320    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2321    index_t* tags=NULL;    index_t* tags=NULL;
# Line 2171  bool MeshAdapter::canTag(int functionSpa Line 2381  bool MeshAdapter::canTag(int functionSpa
2381    }    }
2382  }  }
2383    
2384    AbstractDomain::StatusType MeshAdapter::getStatus() const
2385    {
2386      Finley_Mesh* mesh=m_finleyMesh.get();
2387      return Finley_Mesh_getStatus(mesh);
2388    }
2389    
2390    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2391    {
2392      
2393      Finley_Mesh* mesh=m_finleyMesh.get();
2394      int order =-1;
2395      switch(functionSpaceCode) {
2396       case(Nodes):
2397       case(DegreesOfFreedom):
2398              order=mesh->approximationOrder;
2399              break;
2400       case(ReducedNodes):
2401       case(ReducedDegreesOfFreedom):
2402              order=mesh->reducedApproximationOrder;
2403              break;
2404       case(Elements):
2405       case(FaceElements):
2406       case(Points):
2407       case(ContactElementsZero):
2408       case(ContactElementsOne):
2409              order=mesh->integrationOrder;
2410              break;
2411       case(ReducedElements):
2412       case(ReducedFaceElements):
2413       case(ReducedContactElementsZero):
2414       case(ReducedContactElementsOne):
2415              order=mesh->reducedIntegrationOrder;
2416              break;
2417       default:
2418          stringstream temp;
2419          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2420          throw FinleyAdapterException(temp.str());
2421      }
2422      return order;
2423    }
2424    
2425    bool MeshAdapter::supportsContactElements() const
2426    {
2427      return true;
2428    }
2429    
2430    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2431    {
2432      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2433    }
2434    
2435    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2436    {
2437      Finley_ReferenceElementSet_dealloc(m_refSet);
2438    }
2439    
2440    void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2441    {
2442          const int dim = getDim();
2443          int numPoints=boost::python::extract<int>(points.attr("__len__")());
2444          int numTags=boost::python::extract<int>(tags.attr("__len__")());
2445          Finley_Mesh* mesh=m_finleyMesh.get();
2446          
2447          if  ( (numTags > 0) and ( numPoints !=  numTags ) )
2448         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2449          
2450          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2451          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2452          
2453          for (int i=0;i<numPoints;++i) {
2454           int tag_id=-1;
2455           int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2456           if  ( numComps !=   dim ) {
2457                   stringstream temp;          
2458                   temp << "Error - illegal number of components " << numComps << " for point " << i;              
2459                   throw FinleyAdapterException(temp.str());
2460           }
2461           points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2462           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2463           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2464          
2465           if ( numTags > 0) {
2466                  boost::python::extract<string> ex_str(tags[i]);
2467              if  ( ex_str.check() ) {
2468                  tag_id=getTag( ex_str());
2469              } else {
2470                   boost::python::extract<int> ex_int(tags[i]);
2471                   if ( ex_int.check() ) {
2472                       tag_id=ex_int();
2473                   } else {
2474                    stringstream temp;          
2475                        temp << "Error - unable to extract tag for point " << i;
2476                    throw FinleyAdapterException(temp.str());
2477                  }
2478              }
2479           }      
2480               tags_ptr[i]=tag_id;
2481          }
2482          
2483          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2484          checkPasoError();
2485          
2486          TMPMEMFREE(points_ptr);
2487          TMPMEMFREE(tags_ptr);
2488    }
2489    
2490    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2491    {  
2492        boost::python::list points =  boost::python::list();
2493        boost::python::list tags =  boost::python::list();
2494        points.append(point);
2495        tags.append(tag);
2496        addDiracPoints(points, tags);
2497    }
2498    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2499    {
2500            boost::python::list points =   boost::python::list();
2501            boost::python::list tags =   boost::python::list();
2502            points.append(point);
2503            tags.append(tag);
2504            addDiracPoints(points, tags);
2505    }  
2506  }  // end of namespace  }  // end of namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26