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

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

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

revision 2149 by caltinay, Wed Dec 10 05:08:17 2008 UTC revision 4346 by jfenwick, Tue Apr 2 04:46:45 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2013 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16    #include <pasowrap/PasoException.h>
17    #include <pasowrap/TransportProblemAdapter.h>
18  #include "MeshAdapter.h"  #include "MeshAdapter.h"
19  #include "escript/Data.h"  #include "escript/Data.h"
20  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
21  #ifdef USE_NETCDF  #ifdef USE_NETCDF
22  #include <netcdfcpp.h>  #include <netcdfcpp.h>
23  #endif  #endif
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #include "paso/Paso_MPI.h"  
 #endif  
 extern "C" {  
24  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
 }  
 #include <vector>  
25    
26  #define IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH  256  #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 88  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 99  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 112  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 137  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 314  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 354  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 394  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 436  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 491  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 636  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 681  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 710  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 757  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 773  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 816  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 825  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 850  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 865  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 884  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 1166  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 1247  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 1260  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 1297  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 1334  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 1446  void MeshAdapter::setToSize(escript::Dat Line 1498  void MeshAdapter::setToSize(escript::Dat
1498     checkFinleyError();     checkFinleyError();
1499  }  }
1500    
1501  // sets the location of nodes:  //
1502    // sets the location of nodes
1503    //
1504  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1505  {  {
1506     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1454  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     checkFinleyError();         Finley_Mesh_setCoordinates(mesh,&tmp);
1514  }     } else {
1515           throw FinleyAdapterException("As of version escript3.3 - SetNewX only accepts ContinuousFunction arguments please interpolate.");    
1516  // saves a data array in openDX format:  /*       escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1517  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const         tmp = new_x_inter.getDataC();
1518  {         Finley_Mesh_setCoordinates(mesh,&tmp);*/
    unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;  
    const int num_data=boost::python::extract<int>(arg.attr("__len__")());  
    /* win32 refactor */  
    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;  
    for(int i=0;i<num_data;i++)  
    {  
       names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
    }  
   
    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  
    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  
   
    boost::python::list keys=arg.keys();  
    for (int i=0;i<num_data;++i) {  
       std::string n=boost::python::extract<std::string>(keys[i]);  
       escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);  
       if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)  
          throw FinleyAdapterException("Error in saveDX: Data must be defined on same Domain");  
       data[i]=d.getDataC();  
       ptr_data[i]= &(data[i]);  
       if (n.length()>MAX_namelength-1) {  
          strncpy(names[i],n.c_str(),MAX_namelength-1);  
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
1519     }     }
    Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  
1520     checkFinleyError();     checkFinleyError();
   
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0;i<num_data;i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
   
    return;  
1521  }  }
1522    
1523  // saves a data array in openVTK format:  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1524  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  {
1525  {      if (getMPISize() > 1) {
1526     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;  #ifdef ESYS_MPI
1527     const int num_data=boost::python::extract<int>(arg.attr("__len__")());          index_t myFirstNode=0, myLastNode=0, k=0;
1528     /* win32 refactor */          index_t* globalNodeIndex=0;
1529     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;          Finley_Mesh* mesh_p=m_finleyMesh.get();
1530     for(int i=0;i<num_data;i++)          /*
1531     {           * this method is only used by saveDataCSV which would use the returned
1532        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;           * values for reduced nodes wrongly so this case is disabled for now
1533     }          if (fs_code == FINLEY_REDUCED_NODES)
1534            {
1535                myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1536                myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1537                globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1538            }
1539            else
1540            */
1541            if (fs_code == FINLEY_NODES)
1542            {
1543                myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1544                myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1545                globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1546            }
1547            else
1548            {
1549                throw FinleyAdapterException("Unsupported function space type for ownSample()");
1550            }
1551    
1552     escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;          k=globalNodeIndex[id];
1553     escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;          return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1554    #endif
1555        }
1556        return true;
1557    }
1558    
    boost::python::list keys=arg.keys();  
    for (int i=0;i<num_data;++i) {  
       std::string n=boost::python::extract<std::string>(keys[i]);  
       escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);  
       if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)  
          throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");  
       data[i]=d.getDataC();  
       ptr_data[i]=&(data[i]);  
       if (n.length()>MAX_namelength-1) {  
          strncpy(names[i],n.c_str(),MAX_namelength-1);  
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
    }  
    Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  
1559    
1560     checkFinleyError();  //
1561     /* win32 refactor */  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1562     TMPMEMFREE(data);  //
1563     TMPMEMFREE(ptr_data);  ASM_ptr MeshAdapter::newSystemMatrix(
    for(int i=0;i<num_data;i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
   
    return;  
 }  
                                                                                                                                                                       
                                                                                                                                                                       
 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  
 SystemMatrixAdapter MeshAdapter::newSystemMatrix(  
1564                                                   const int row_blocksize,                                                   const int row_blocksize,
1565                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1566                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1566  SystemMatrixAdapter MeshAdapter::newSyst Line 1575  SystemMatrixAdapter MeshAdapter::newSyst
1575        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1576     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1577     if (col_domain!=*this)     if (col_domain!=*this)
1578        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of column function space does not match the domain of matrix generator.");
1579     // is the function space type right     // is the function space type right
1580     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1581        reduceRowOrder=0;        reduceRowOrder=0;
# Line 1594  SystemMatrixAdapter MeshAdapter::newSyst Line 1603  SystemMatrixAdapter MeshAdapter::newSyst
1603  #endif  #endif
1604     }     }
1605     else {     else {
1606        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1607     }     }
1608     checkPasoError();     checkPasoError();
1609     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1610     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1611       return ASM_ptr(sma);
1612    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1613  }  }
1614    
1615    //
1616  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1617  TransportProblemAdapter MeshAdapter::newTransportProblem(  //
1618                                                           const double theta,  ATP_ptr MeshAdapter::newTransportProblem(
1619                                                           const int blocksize,                                                           const int blocksize,
1620                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1621                                                           const int type) const                                                           const int type) const
# Line 1624  TransportProblemAdapter MeshAdapter::new Line 1637  TransportProblemAdapter MeshAdapter::new
1637    
1638     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1639     checkFinleyError();     checkFinleyError();
1640     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1641     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1642     checkPasoError();     checkPasoError();
1643     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1644     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1645       return ATP_ptr(tpa);
1646    //   return TransportProblemAdapter(transportProblem, blocksize,functionspace);
1647  }  }
1648    
1649  //  //
# Line 1667  bool MeshAdapter::isCellOriented(int fun Line 1682  bool MeshAdapter::isCellOriented(int fun
1682     return false;     return false;
1683  }  }
1684    
1685    bool
1686    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1687    {
1688       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1689        class 1: DOF <-> Nodes
1690        class 2: ReducedDOF <-> ReducedNodes
1691        class 3: Points
1692        class 4: Elements
1693        class 5: ReducedElements
1694        class 6: FaceElements
1695        class 7: ReducedFaceElements
1696        class 8: ContactElementZero <-> ContactElementOne
1697        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1698    
1699       There is also a set of lines. Interpolation is possible down a line but not between lines.
1700       class 1 and 2 belong to all lines so aren't considered.
1701        line 0: class 3
1702        line 1: class 4,5
1703        line 2: class 6,7
1704        line 3: class 8,9
1705    
1706       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1707       eg hasnodes is true if we have at least one instance of Nodes.
1708       */
1709        if (fs.empty())
1710        {
1711            return false;
1712        }
1713        vector<int> hasclass(10);
1714        vector<int> hasline(4);
1715        bool hasnodes=false;
1716        bool hasrednodes=false;
1717        bool hascez=false;
1718        bool hasrcez=false;
1719        for (int i=0;i<fs.size();++i)
1720        {
1721        switch(fs[i])
1722        {
1723        case(Nodes):   hasnodes=true;   // no break is deliberate
1724        case(DegreesOfFreedom):
1725            hasclass[1]=1;
1726            break;
1727        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1728        case(ReducedDegreesOfFreedom):
1729            hasclass[2]=1;
1730            break;
1731        case(Points):
1732            hasline[0]=1;
1733            hasclass[3]=1;
1734            break;
1735        case(Elements):
1736            hasclass[4]=1;
1737            hasline[1]=1;
1738            break;
1739        case(ReducedElements):
1740            hasclass[5]=1;
1741            hasline[1]=1;
1742            break;
1743        case(FaceElements):
1744            hasclass[6]=1;
1745            hasline[2]=1;
1746            break;
1747        case(ReducedFaceElements):
1748            hasclass[7]=1;
1749            hasline[2]=1;
1750            break;
1751        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1752        case(ContactElementsOne):
1753            hasclass[8]=1;
1754            hasline[3]=1;
1755            break;
1756        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1757        case(ReducedContactElementsOne):
1758            hasclass[9]=1;
1759            hasline[3]=1;
1760            break;
1761        default:
1762            return false;
1763        }
1764        }
1765        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1766        // fail if we have more than one leaf group
1767    
1768        if (totlines>1)
1769        {
1770        return false;   // there are at least two branches we can't interpolate between
1771        }
1772        else if (totlines==1)
1773        {
1774        if (hasline[0]==1)      // we have points
1775        {
1776            resultcode=Points;
1777        }
1778        else if (hasline[1]==1)
1779        {
1780            if (hasclass[5]==1)
1781            {
1782            resultcode=ReducedElements;
1783            }
1784            else
1785            {
1786            resultcode=Elements;
1787            }
1788        }
1789        else if (hasline[2]==1)
1790        {
1791            if (hasclass[7]==1)
1792            {
1793            resultcode=ReducedFaceElements;
1794            }
1795            else
1796            {
1797            resultcode=FaceElements;
1798            }
1799        }
1800        else    // so we must be in line3
1801        {
1802            if (hasclass[9]==1)
1803            {
1804            // need something from class 9
1805            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1806            }
1807            else
1808            {
1809            // something from class 8
1810            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1811            }
1812        }
1813        }
1814        else    // totlines==0
1815        {
1816        if (hasclass[2]==1)
1817        {
1818            // something from class 2
1819            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1820        }
1821        else
1822        {   // something from class 1
1823            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1824        }
1825        }
1826        return true;
1827    }
1828    
1829  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1830  {  {
1831     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1832     case(Nodes):     case(Nodes):
1833     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1834     case(Nodes):      case(Nodes):
1835     case(ReducedNodes):      case(ReducedNodes):
1836     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1837     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1838     case(Elements):      case(Elements):
1839     case(ReducedElements):      case(ReducedElements):
1840     case(FaceElements):      case(FaceElements):
1841     case(ReducedFaceElements):      case(ReducedFaceElements):
1842     case(Points):      case(Points):
1843     case(ContactElementsZero):      case(ContactElementsZero):
1844     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1845     case(ContactElementsOne):      case(ContactElementsOne):
1846     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1847     return true;      return true;
1848     default:      default:
1849        stringstream temp;            stringstream temp;
1850        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;
1851        throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1852     }     }
1853     break;     break;
1854     case(ReducedNodes):     case(ReducedNodes):
1855     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1856     case(ReducedNodes):      case(ReducedNodes):
1857     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1858     case(Elements):      case(Elements):
1859     case(ReducedElements):      case(ReducedElements):
1860     case(FaceElements):      case(FaceElements):
1861     case(ReducedFaceElements):      case(ReducedFaceElements):
1862     case(Points):      case(Points):
1863     case(ContactElementsZero):      case(ContactElementsZero):
1864     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1865     case(ContactElementsOne):      case(ContactElementsOne):
1866     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1867     return true;      return true;
1868     case(Nodes):      case(Nodes):
1869     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1870     return false;      return false;
1871     default:      default:
1872        stringstream temp;          stringstream temp;
1873        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;
1874        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1875     }     }
1876     break;     break;
1877     case(Elements):     case(Elements):
1878     if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1879        return true;        return true;
1880     } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1881        return true;        return true;
1882     } else {          } else {
1883        return false;            return false;
1884     }          }
1885     case(ReducedElements):     case(ReducedElements):
1886     if (functionSpaceType_target==ReducedElements) {      if (functionSpaceType_target==ReducedElements) {
1887        return true;        return true;
1888     } else {      } else {
1889        return false;            return false;
1890     }      }
1891     case(FaceElements):     case(FaceElements):
1892     if (functionSpaceType_target==FaceElements) {      if (functionSpaceType_target==FaceElements) {
1893        return true;              return true;
1894     } else if (functionSpaceType_target==ReducedFaceElements) {      } else if (functionSpaceType_target==ReducedFaceElements) {
1895        return true;              return true;
1896     } else {      } else {
1897        return false;              return false;
1898     }      }
1899     case(ReducedFaceElements):     case(ReducedFaceElements):
1900     if (functionSpaceType_target==ReducedFaceElements) {      if (functionSpaceType_target==ReducedFaceElements) {
1901        return true;              return true;
1902     } else {      } else {
1903        return false;          return false;
1904     }      }
1905     case(Points):     case(Points):
1906     if (functionSpaceType_target==Points) {      if (functionSpaceType_target==Points) {
1907        return true;              return true;
1908     } else {      } else {
1909        return false;              return false;
1910     }      }
1911     case(ContactElementsZero):     case(ContactElementsZero):
1912     case(ContactElementsOne):     case(ContactElementsOne):
1913     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1914        return true;              return true;
1915     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1916        return true;              return true;
1917     } else {      } else {
1918        return false;              return false;
1919     }      }
1920     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1921     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1922     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1923        return true;              return true;
1924     } else {      } else {
1925        return false;              return false;
1926     }      }
1927     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1928     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1929     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1930     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1931     case(Nodes):      case(Nodes):
1932     case(ReducedNodes):      case(ReducedNodes):
1933     case(Elements):      case(Elements):
1934     case(ReducedElements):      case(ReducedElements):
1935     case(Points):      case(Points):
1936     case(FaceElements):      case(FaceElements):
1937     case(ReducedFaceElements):      case(ReducedFaceElements):
1938     case(ContactElementsZero):      case(ContactElementsZero):
1939     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1940     case(ContactElementsOne):      case(ContactElementsOne):
1941     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1942     return true;      return true;
1943     default:      default:
1944        stringstream temp;          stringstream temp;
1945        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;
1946        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1947     }      }
1948     break;      break;
1949     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1950     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
1951     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1952     case(ReducedNodes):      case(ReducedNodes):
1953     case(Elements):      case(Elements):
1954     case(ReducedElements):      case(ReducedElements):
1955     case(FaceElements):      case(FaceElements):
1956     case(ReducedFaceElements):      case(ReducedFaceElements):
1957     case(Points):      case(Points):
1958     case(ContactElementsZero):      case(ContactElementsZero):
1959     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1960     case(ContactElementsOne):      case(ContactElementsOne):
1961     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1962     return true;      return true;
1963     case(Nodes):      case(Nodes):
1964     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1965     return false;      return false;
1966     default:      default:
1967        stringstream temp;          stringstream temp;
1968        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;
1969        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1970     }      }
1971     break;      break;
1972     default:     default:
1973        stringstream temp;        stringstream temp;
1974        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 1820  bool MeshAdapter::probeInterpolationOnDo Line 1979  bool MeshAdapter::probeInterpolationOnDo
1979     return false;     return false;
1980  }  }
1981    
1982    signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1983    {
1984        if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1985        {
1986            return 1;
1987        }
1988        else if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1989        {
1990            return -1;
1991        }
1992        return 0;
1993    }  
1994      
1995  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1996  {  {
1997     return false;     return false;
# Line 1842  bool MeshAdapter::operator!=(const Abstr Line 2014  bool MeshAdapter::operator!=(const Abstr
2014    
2015  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
2016  {  {
2017     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2018     checkPasoError();     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2019     return out;             package, symmetry, mesh->MPIInfo);
2020  }  }
2021    
2022  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
2023  {  {
2024     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2025     checkPasoError();     return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2026     return out;             package, symmetry, mesh->MPIInfo);
2027  }  }
2028    
2029  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2030  {  {
2031     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2032  }  }
2033    
2034  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2035  {  {
2036     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2037  }  }
2038    
2039  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2040  {  {
2041     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2042  }  }
2043    
2044  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2045  {  {
2046     int *out = NULL;     int *out = NULL;
2047     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2027  void MeshAdapter::setTags(const int func Line 2200  void MeshAdapter::setTags(const int func
2200     return;     return;
2201  }  }
2202    
2203  void MeshAdapter::setTagMap(const std::string& name,  int tag)  void MeshAdapter::setTagMap(const string& name,  int tag)
2204  {  {
2205     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2206     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2207     checkPasoError();     checkFinleyError();
2208     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2209  }  }
2210    
2211  int MeshAdapter::getTag(const std::string& name) const  int MeshAdapter::getTag(const string& name) const
2212  {  {
2213     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2214     int tag=0;     int tag=0;
2215     tag=Finley_Mesh_getTag(mesh, name.c_str());     tag=Finley_Mesh_getTag(mesh, name.c_str());
2216     checkPasoError();     checkFinleyError();
2217     // throwStandardException("MeshAdapter::getTag is not implemented.");     // throwStandardException("MeshAdapter::getTag is not implemented.");
2218     return tag;     return tag;
2219  }  }
2220    
2221  bool MeshAdapter::isValidTagName(const std::string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2222  {  {
2223     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2224     return Finley_Mesh_isValidTagName(mesh,name.c_str());     return Finley_Mesh_isValidTagName(mesh,name.c_str());
2225  }  }
2226    
2227  std::string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2228  {  {
2229     stringstream temp;     stringstream temp;
2230     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2105  int MeshAdapter::getNumberOfTagsInUse(in Line 2278  int MeshAdapter::getNumberOfTagsInUse(in
2278    }    }
2279    return numTags;    return numTags;
2280  }  }
2281  int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  
2282    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2283  {  {
2284    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2285    index_t* tags=NULL;    index_t* tags=NULL;
# Line 2171  bool MeshAdapter::canTag(int functionSpa Line 2345  bool MeshAdapter::canTag(int functionSpa
2345    }    }
2346  }  }
2347    
2348    AbstractDomain::StatusType MeshAdapter::getStatus() const
2349    {
2350      Finley_Mesh* mesh=m_finleyMesh.get();
2351      return Finley_Mesh_getStatus(mesh);
2352    }
2353    
2354    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2355    {
2356      
2357      Finley_Mesh* mesh=m_finleyMesh.get();
2358      int order =-1;
2359      switch(functionSpaceCode) {
2360       case(Nodes):
2361       case(DegreesOfFreedom):
2362              order=mesh->approximationOrder;
2363              break;
2364       case(ReducedNodes):
2365       case(ReducedDegreesOfFreedom):
2366              order=mesh->reducedApproximationOrder;
2367              break;
2368       case(Elements):
2369       case(FaceElements):
2370       case(Points):
2371       case(ContactElementsZero):
2372       case(ContactElementsOne):
2373              order=mesh->integrationOrder;
2374              break;
2375       case(ReducedElements):
2376       case(ReducedFaceElements):
2377       case(ReducedContactElementsZero):
2378       case(ReducedContactElementsOne):
2379              order=mesh->reducedIntegrationOrder;
2380              break;
2381       default:
2382          stringstream temp;
2383          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2384          throw FinleyAdapterException(temp.str());
2385      }
2386      return order;
2387    }
2388    
2389    bool MeshAdapter::supportsContactElements() const
2390    {
2391      return true;
2392    }
2393    
2394    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2395    {
2396      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2397    }
2398    
2399    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2400    {
2401      Finley_ReferenceElementSet_dealloc(m_refSet);
2402    }
2403    
2404    // points will be flattened
2405    void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2406    {
2407          const int dim = getDim();
2408          int numPoints=points.size()/dim;
2409          int numTags=tags.size();
2410          Finley_Mesh* mesh=m_finleyMesh.get();
2411          
2412          if ( points.size() % dim != 0 )
2413          {
2414        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2415          }
2416          
2417          if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2418         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2419          
2420          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2421          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2422          
2423          for (int i=0;i<numPoints;++i) {
2424           points_ptr[ i * dim     ] = points[i * dim ];
2425           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2426           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2427               tags_ptr[i]=tags[i];
2428          }
2429          
2430          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2431          checkFinleyError();
2432          
2433          TMPMEMFREE(points_ptr);
2434          TMPMEMFREE(tags_ptr);
2435    }
2436    
2437    
2438    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2439    // {
2440    //       const int dim = getDim();
2441    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2442    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2443    //       Finley_Mesh* mesh=m_finleyMesh.get();
2444    //      
2445    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2446    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2447    //      
2448    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2449    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2450    //      
2451    //       for (int i=0;i<numPoints;++i) {
2452    //     int tag_id=-1;
2453    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2454    //     if  ( numComps !=   dim ) {
2455    //                stringstream temp;            
2456    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2457    //                throw FinleyAdapterException(temp.str());
2458    //     }
2459    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2460    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2461    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2462    //    
2463    //     if ( numTags > 0) {
2464    //            boost::python::extract<string> ex_str(tags[i]);
2465    //        if  ( ex_str.check() ) {
2466    //            tag_id=getTag( ex_str());
2467    //        } else {
2468    //             boost::python::extract<int> ex_int(tags[i]);
2469    //             if ( ex_int.check() ) {
2470    //                 tag_id=ex_int();
2471    //             } else {
2472    //              stringstream temp;          
2473    //                  temp << "Error - unable to extract tag for point " << i;
2474    //              throw FinleyAdapterException(temp.str());
2475    //            }
2476    //        }
2477    //     }      
2478    //            tags_ptr[i]=tag_id;
2479    //       }
2480    //      
2481    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2482    //       checkPasoError();
2483    //      
2484    //       TMPMEMFREE(points_ptr);
2485    //       TMPMEMFREE(tags_ptr);
2486    // }
2487    
2488    /*
2489    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2490    {  
2491        boost::python::list points =  boost::python::list();
2492        boost::python::list tags =  boost::python::list();
2493        points.append(point);
2494        tags.append(tag);
2495        addDiracPoints(points, tags);
2496    }
2497    */
2498    
2499    /*
2500    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2501    {
2502            boost::python::list points =   boost::python::list();
2503            boost::python::list tags =   boost::python::list();
2504            points.append(point);
2505            tags.append(tag);
2506            addDiracPoints(points, tags);
2507    }
2508    */
2509  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.2149  
changed lines
  Added in v.4346

  ViewVC Help
Powered by ViewVC 1.1.26