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

Legend:
Removed from v.2315  
changed lines
  Added in v.3914

  ViewVC Help
Powered by ViewVC 1.1.26