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

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

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

revision 1748 by ksteube, Wed Sep 3 06:10:39 2008 UTC revision 3998 by caltinay, Thu Sep 27 01:17:28 2012 UTC
# Line 1  Line 1 
1    
2  /* $Id$ */  /*****************************************************************************
3    *
4  /*******************************************************  * Copyright (c) 2003-2012 by University of Queensland
5   *  * http://www.uq.edu.au
6   *           Copyright 2003-2007 by ACceSS MNRF  *
7   *       Copyright 2007 by University of Queensland  * Primary Business: Queensland, Australia
8   *  * Licensed under the Open Software License version 3.0
9   *                http://esscc.uq.edu.au  * http://www.opensource.org/licenses/osl-3.0.php
10   *        Primary Business: Queensland, Australia  *
11   *  Licensed under the Open Software License version 3.0  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12   *     http://www.opensource.org/licenses/osl-3.0.php  * 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"
# Line 20  Line 22 
22  #include <netcdfcpp.h>  #include <netcdfcpp.h>
23  #endif  #endif
24  extern "C" {  extern "C" {
25  #include "escript/blocktimer.h"  #include "esysUtils/blocktimer.h"
26  }  }
 #include <vector>  
27    
28  #define IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH  256  #include <boost/python/import.hpp>
29    
30  using namespace std;  using namespace std;
31  using namespace escript;  using namespace escript;
32    using namespace paso;
33    
34  namespace finley {  namespace finley {
35    
# Line 83  int MeshAdapter::getMPIRank() const Line 85  int MeshAdapter::getMPIRank() const
85  {  {
86     return m_finleyMesh.get()->MPIInfo->rank;     return m_finleyMesh.get()->MPIInfo->rank;
87  }  }
88    void MeshAdapter::MPIBarrier() const
89    {
90    #ifdef ESYS_MPI
91       MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
92    #endif
93       return;
94    }
95    bool MeshAdapter::onMasterProcessor() const
96    {
97       return m_finleyMesh.get()->MPIInfo->rank == 0;
98    }
99    
100    
101    #ifdef ESYS_MPI
102      MPI_Comm
103    #else
104      unsigned int
105    #endif
106    MeshAdapter::getMPIComm() const
107    {
108    #ifdef ESYS_MPI
109        return m_finleyMesh->MPIInfo->comm;
110    #else
111        return 0;
112    #endif
113    }
114    
115    
116  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
117     return m_finleyMesh.get();     return m_finleyMesh.get();
118  }  }
119    
120  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const string& fileName) const
121  {  {
122     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;
123     strcpy(fName,fileName.c_str());     strcpy(fName,fileName.c_str());
# Line 98  void MeshAdapter::write(const std::strin Line 126  void MeshAdapter::write(const std::strin
126     TMPMEMFREE(fName);     TMPMEMFREE(fName);
127  }  }
128    
129  void MeshAdapter::Print_Mesh_Info(const bool full=false) const  void MeshAdapter::Print_Mesh_Info(const bool full) const
130  {  {
131     Finley_PrintMesh_Info(m_finleyMesh.get(), full);     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
132  }  }
133    
134  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const string& fileName) const
135  {  {
136  #ifdef USE_NETCDF  #ifdef USE_NETCDF
137     const NcDim* ncdims[12];     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
138     NcVar *ids;     NcVar *ids;
139     int *int_ptr;     int *int_ptr;
140     Finley_Mesh *mesh = m_finleyMesh.get();     Finley_Mesh *mesh = m_finleyMesh.get();
# Line 123  void MeshAdapter::dump(const std::string Line 151  void MeshAdapter::dump(const std::string
151     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
152     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
153     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
154    #ifdef ESYS_MPI
155       MPI_Status status;
156    #endif
157    
158     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),  /* Incoming token indicates it's my turn to write */
159    #ifdef ESYS_MPI
160       if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
161    #endif
162    
163       char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
164                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
165    
166     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
167     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
168     if (tag_map) {     num_Tags = 0;
169        while (tag_map) {     while (tag_map) {
170           num_Tags++;        num_Tags++;
171           tag_map=tag_map->next;        tag_map=tag_map->next;
       }  
172     }     }
173    
174     // NetCDF error handler     // NetCDF error handler
175     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
176     // Create the file.     // Create the file.
177     NcFile dataFile(newFileName, NcFile::Replace);     NcFile dataFile(newFileName, NcFile::Replace);
178       string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
179     // check if writing was successful     // check if writing was successful
180     if (!dataFile.is_valid())     if (!dataFile.is_valid())
181        throw DataException("Error - MeshAdapter::dump: opening of NetCDF file for output failed: " + *newFileName);        throw DataException(msgPrefix+"Open file for output");
182    
183     // Define dimensions (num_Elements and dim_Elements are identical, dim_Elements only appears if > 0)     // Define dimensions (num_Elements and dim_Elements are identical,
184       // dim_Elements only appears if > 0)
185     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
186        throw DataException("Error - MeshAdapter::dump: appending dimension numNodes to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numNodes)");
187     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
188        throw DataException("Error - MeshAdapter::dump: appending dimension numDim to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numDim)");
189     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)) )
190        throw DataException("Error - MeshAdapter::dump: appending dimension mpi_size to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(mpi_size)");
191     if (num_Elements>0)     if (num_Elements>0)
192        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
193           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements)");
194     if (num_FaceElements>0)     if (num_FaceElements>0)
195        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
196           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
197     if (num_ContactElements>0)     if (num_ContactElements>0)
198        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
199           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements)");
200     if (num_Points>0)     if (num_Points>0)
201        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
202           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Points to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Points)");
203     if (num_Elements>0)     if (num_Elements>0)
204        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
205           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
206     if (num_FaceElements>0)     if (num_FaceElements>0)
207        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
208           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
209     if (num_ContactElements>0)     if (num_ContactElements>0)
210        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
211           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements_numNodes)");
212     if (num_Tags>0)     if (num_Tags>0)
213        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
214           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Tags)");
215    
216     // Attributes: MPI size, MPI rank, Name, order, reduced_order     // Attributes: MPI size, MPI rank, Name, order, reduced_order
217     if (!dataFile.add_att("mpi_size", mpi_size) )     if (!dataFile.add_att("mpi_size", mpi_size) )
218        throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_size)");
219     if (!dataFile.add_att("mpi_rank", mpi_rank) )     if (!dataFile.add_att("mpi_rank", mpi_rank) )
220        throw DataException("Error - MeshAdapter::dump: appending mpi_rank to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_rank)");
221     if (!dataFile.add_att("Name",mesh->Name) )     if (!dataFile.add_att("Name",mesh->Name) )
222        throw DataException("Error - MeshAdapter::dump: appending Name to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Name)");
223     if (!dataFile.add_att("numDim",numDim) )     if (!dataFile.add_att("numDim",numDim) )
224        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
225     if (!dataFile.add_att("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
226        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
227     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
228        throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(reduced_order)");
229     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
230        throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(numNodes)");
231     if (!dataFile.add_att("num_Elements",num_Elements) )     if (!dataFile.add_att("num_Elements",num_Elements) )
232        throw DataException("Error - MeshAdapter::dump: appending num_Elements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements)");
233     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
234        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements)");
235     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
236        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements)");
237     if (!dataFile.add_att("num_Points",num_Points) )     if (!dataFile.add_att("num_Points",num_Points) )
238        throw DataException("Error - MeshAdapter::dump: appending num_Points to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Points)");
239     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
240        throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
241     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
242        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
243     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
244        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements_numNodes)");
245     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
246        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Elements_TypeId)");
247     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
248        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
249     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
250        throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(ContactElements_TypeId)");
251     if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
252        throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Points_TypeId)");
253     if (!dataFile.add_att("num_Tags", num_Tags) )     if (!dataFile.add_att("num_Tags", num_Tags) )
254        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Tags)");
255    
256     // // // // // Nodes // // // // //     // // // // // Nodes // // // // //
257    
258     // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)     // Nodes nodeDistribution
259       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
260          throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
261       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
262       if (! (ids->put(int_ptr, mpi_size+1)) )
263          throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
264    
265       // Nodes degreesOfFreedomDistribution
266       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
267          throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
268       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
269       if (! (ids->put(int_ptr, mpi_size+1)) )
270          throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
271    
272       // Only write nodes if non-empty because NetCDF doesn't like empty arrays
273       // (it treats them as NC_UNLIMITED)
274     if (numNodes>0) {     if (numNodes>0) {
275    
276        // Nodes Id        // Nodes Id
277        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
278           throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Id)");
279        int_ptr = &mesh->Nodes->Id[0];        int_ptr = &mesh->Nodes->Id[0];
280        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
281           throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Id)");
282    
283        // Nodes Tag        // Nodes Tag
284        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
285           throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Tag)");
286        int_ptr = &mesh->Nodes->Tag[0];        int_ptr = &mesh->Nodes->Tag[0];
287        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
288           throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Tag)");
289    
290        // Nodes gDOF        // Nodes gDOF
291        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
292           throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
293        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
294        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
295           throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gDOF)");
296    
297        // Nodes global node index        // Nodes global node index
298        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
299           throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gNI)");
300        int_ptr = &mesh->Nodes->globalNodesIndex[0];        int_ptr = &mesh->Nodes->globalNodesIndex[0];
301        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
302           throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gNI)");
303    
304        // Nodes grDof        // Nodes grDof
305        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
306           throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
307        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
308        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
309           throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grDfI)");
310    
311        // Nodes grNI        // Nodes grNI
312        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
313           throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grNI)");
314        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
315        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
316           throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grNI)");
317    
318        // Nodes Coordinates        // Nodes Coordinates
319        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
320           throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
321        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
322           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);  
323    
324     }     }
325    
# Line 282  void MeshAdapter::dump(const std::string Line 327  void MeshAdapter::dump(const std::string
327    
328     if (num_Elements>0) {     if (num_Elements>0) {
329    
       // Temp storage to gather node IDs  
       int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);  
   
330        // Elements_Id        // Elements_Id
331        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
332           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Id)");
333        int_ptr = &mesh->Elements->Id[0];        int_ptr = &mesh->Elements->Id[0];
334        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
335           throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Id)");
336    
337        // Elements_Tag        // Elements_Tag
338        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
339           throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Tag)");
340        int_ptr = &mesh->Elements->Tag[0];        int_ptr = &mesh->Elements->Tag[0];
341        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
342           throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Tag)");
343    
344        // Elements_Owner        // Elements_Owner
345        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
346           throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Owner)");
347        int_ptr = &mesh->Elements->Owner[0];        int_ptr = &mesh->Elements->Owner[0];
348        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
349           throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Owner)");
350    
351        // Elements_Color        // Elements_Color
352        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
353           throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Color)");
354        int_ptr = &mesh->Elements->Color[0];        int_ptr = &mesh->Elements->Color[0];
355        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
356           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Color)");
357    
358        // Elements_Nodes        // Elements_Nodes
       for (int i=0; i<num_Elements; i++)  
          for (int j=0; j<num_Elements_numNodes; j++)  
             Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)] = mesh->Nodes->Id[mesh->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]];  
359        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
360           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Nodes)");
361        if (! (ids->put(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes)) )        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
362           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Nodes)");
   
       TMPMEMFREE(Elements_Nodes);  
363    
364     }     }
365    
# Line 330  void MeshAdapter::dump(const std::string Line 367  void MeshAdapter::dump(const std::string
367    
368     if (num_FaceElements>0) {     if (num_FaceElements>0) {
369    
       // Temp storage to gather node IDs  
       int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);  
   
370        // FaceElements_Id        // FaceElements_Id
371        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
372           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Id)");
373        int_ptr = &mesh->FaceElements->Id[0];        int_ptr = &mesh->FaceElements->Id[0];
374        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
375           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Id)");
376    
377        // FaceElements_Tag        // FaceElements_Tag
378        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
379           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
380        int_ptr = &mesh->FaceElements->Tag[0];        int_ptr = &mesh->FaceElements->Tag[0];
381        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
382           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Tag)");
383    
384        // FaceElements_Owner        // FaceElements_Owner
385        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
386           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
387        int_ptr = &mesh->FaceElements->Owner[0];        int_ptr = &mesh->FaceElements->Owner[0];
388        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
389           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Owner)");
390    
391        // FaceElements_Color        // FaceElements_Color
392        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
393           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Color)");
394        int_ptr = &mesh->FaceElements->Color[0];        int_ptr = &mesh->FaceElements->Color[0];
395        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
396           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Color)");
397    
398        // FaceElements_Nodes        // FaceElements_Nodes
       for (int i=0; i<num_FaceElements; i++)  
          for (int j=0; j<num_FaceElements_numNodes; j++)  
             FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = mesh->Nodes->Id[mesh->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]];  
399        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
400           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
401        if (! (ids->put(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
402           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Nodes)");
   
       TMPMEMFREE(FaceElements_Nodes);  
403    
404     }     }
405    
# Line 378  void MeshAdapter::dump(const std::string Line 407  void MeshAdapter::dump(const std::string
407    
408     if (num_ContactElements>0) {     if (num_ContactElements>0) {
409    
       // Temp storage to gather node IDs  
       int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);  
   
410        // ContactElements_Id        // ContactElements_Id
411        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
412           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Id)");
413        int_ptr = &mesh->ContactElements->Id[0];        int_ptr = &mesh->ContactElements->Id[0];
414        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
415           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Id)");
416    
417        // ContactElements_Tag        // ContactElements_Tag
418        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
419           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Tag)");
420        int_ptr = &mesh->ContactElements->Tag[0];        int_ptr = &mesh->ContactElements->Tag[0];
421        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
422           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Tag)");
423    
424        // ContactElements_Owner        // ContactElements_Owner
425        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
426           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Owner)");
427        int_ptr = &mesh->ContactElements->Owner[0];        int_ptr = &mesh->ContactElements->Owner[0];
428        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
429           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Owner)");
430    
431        // ContactElements_Color        // ContactElements_Color
432        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
433           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Color)");
434        int_ptr = &mesh->ContactElements->Color[0];        int_ptr = &mesh->ContactElements->Color[0];
435        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
436           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Color)");
437    
438        // ContactElements_Nodes        // ContactElements_Nodes
       for (int i=0; i<num_ContactElements; i++)  
          for (int j=0; j<num_ContactElements_numNodes; j++)  
             ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)] = mesh->Nodes->Id[mesh->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]];  
439        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
440           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Nodes)");
441        if (! (ids->put(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
442           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Nodes)");
   
       TMPMEMFREE(ContactElements_Nodes);  
443    
444     }     }
445    
# Line 428  void MeshAdapter::dump(const std::string Line 449  void MeshAdapter::dump(const std::string
449    
450        fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");        fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");
451    
       // Temp storage to gather node IDs  
       int *Points_Nodes = TMPMEMALLOC(num_Points,int);  
   
452        // Points_Id        // Points_Id
453        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
454           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Id)");
455        int_ptr = &mesh->Points->Id[0];        int_ptr = &mesh->Points->Id[0];
456        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
457           throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Id)");
458    
459        // Points_Tag        // Points_Tag
460        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
461           throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Tag)");
462        int_ptr = &mesh->Points->Tag[0];        int_ptr = &mesh->Points->Tag[0];
463        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
464           throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Tag)");
465    
466        // Points_Owner        // Points_Owner
467        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
468           throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Owner)");
469        int_ptr = &mesh->Points->Owner[0];        int_ptr = &mesh->Points->Owner[0];
470        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
471           throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Owner)");
472    
473        // Points_Color        // Points_Color
474        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
475           throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Color)");
476        int_ptr = &mesh->Points->Color[0];        int_ptr = &mesh->Points->Color[0];
477        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
478           throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Color)");
479    
480        // Points_Nodes        // Points_Nodes
481        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
       for (int i=0; i<num_Points; i++)  
          Points_Nodes[i] = mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]];  
482        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
483           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Nodes)");
484        if (! (ids->put(&(Points_Nodes[0]), num_Points)) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
485           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Nodes)");
   
       TMPMEMFREE(Points_Nodes);  
486    
487     }     }
488    
# Line 492  void MeshAdapter::dump(const std::string Line 506  void MeshAdapter::dump(const std::string
506    
507        // Tags_keys        // Tags_keys
508        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
509           throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Tags_keys)");
510        int_ptr = &Tags_keys[0];        int_ptr = &Tags_keys[0];
511        if (! (ids->put(int_ptr, num_Tags)) )        if (! (ids->put(int_ptr, num_Tags)) )
512           throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Tags_keys)");
513    
514        // Tags_names_*        // Tags_names_*
515        // 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
516        // 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
517          // manual doesn't tell how to do an array of strings
518        tag_map = mesh->TagMap;        tag_map = mesh->TagMap;
519        if (tag_map) {        if (tag_map) {
520           int i = 0;           int i = 0;
521           while (tag_map) {           while (tag_map) {
522              sprintf(name_temp, "Tags_name_%d", i);              sprintf(name_temp, "Tags_name_%d", i);
523              if (!dataFile.add_att(name_temp, tag_map->name) )              if (!dataFile.add_att(name_temp, tag_map->name) )
524                 throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);                 throw DataException(msgPrefix+"add_att(Tags_names_XX)");
525              tag_map=tag_map->next;              tag_map=tag_map->next;
526              i++;              i++;
527           }           }
528        }        }
529    
530        TMPMEMFREE(Tags_keys);        TMPMEMFREE(Tags_keys);
   
531     }     }
532    
533    /* Send token to next MPI process so he can take his turn */
534    #ifdef ESYS_MPI
535       if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
536    #endif
537    
538     // NetCDF file is closed by destructor of NcFile object     // NetCDF file is closed by destructor of NcFile object
539    
540  #else  #else
541     Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");     Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
542  #endif  /* USE_NETCDF */  #endif  /* USE_NETCDF */
# Line 632  int MeshAdapter::getReducedSolutionCode( Line 651  int MeshAdapter::getReducedSolutionCode(
651     return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
652  }  }
653    
654  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
655  {  {
656     return Points;     return Points;
657  }  }
# Line 642  int MeshAdapter::getDiracDeltaFunctionCo Line 661  int MeshAdapter::getDiracDeltaFunctionCo
661  //  //
662  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
663  {  {
664     int numDim=Finley_Mesh_getDim(m_finleyMesh.get());     Finley_Mesh* mesh=m_finleyMesh.get();
665       int numDim=Finley_Mesh_getDim(mesh);
666     checkFinleyError();     checkFinleyError();
667     return numDim;     return numDim;
668  }  }
669    
670  //  //
671    // Return the number of data points summed across all MPI processes
672    //
673    int MeshAdapter::getNumDataPointsGlobal() const
674    {
675       return Finley_NodeFile_getGlobalNumNodes(m_finleyMesh.get()->Nodes);
676    }
677    
678    //
679  // return the number of data points per sample and the number of samples  // return the number of data points per sample and the number of samples
680  // needed to represent data on a parts of the mesh.  // needed to represent data on a parts of the mesh.
681  //  //
# Line 668  pair<int,int> MeshAdapter::getDataShape( Line 696  pair<int,int> MeshAdapter::getDataShape(
696     case(Elements):     case(Elements):
697     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
698        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
699        numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
700     }     }
701     break;     break;
702     case(ReducedElements):     case(ReducedElements):
703     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
704        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
705        numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
706     }     }
707     break;     break;
708     case(FaceElements):     case(FaceElements):
709     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
710        numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
711        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
712     }     }
713     break;     break;
714     case(ReducedFaceElements):     case(ReducedFaceElements):
715     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
716        numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
717        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
718     }     }
719     break;     break;
# Line 697  pair<int,int> MeshAdapter::getDataShape( Line 725  pair<int,int> MeshAdapter::getDataShape(
725     break;     break;
726     case(ContactElementsZero):     case(ContactElementsZero):
727     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
728        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
729        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
730     }     }
731     break;     break;
732     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
733     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
734        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
735        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
736     }     }
737     break;     break;
738     case(ContactElementsOne):     case(ContactElementsOne):
739     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
740        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
741        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
742     }     }
743     break;     break;
744     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
745     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
746        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
747        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
748     }     }
749     break;     break;
# Line 744  pair<int,int> MeshAdapter::getDataShape( Line 772  pair<int,int> MeshAdapter::getDataShape(
772  // 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:
773  //  //
774  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
775                                   SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
776                                   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,
777                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
778                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact,
779                                     const escript::Data& d_dirac,const escript::Data& y_dirac) const
780  {  {
781       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
782       if (smat==0)
783       {
784        throw FinleyAdapterException("finley only supports Paso system matrices.");
785       }
786     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
787     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
788     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 760  void MeshAdapter::addPDEToSystem( Line 794  void MeshAdapter::addPDEToSystem(
794     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
795     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
796     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
797       escriptDataC _d_dirac=d_dirac.getDataC();
798       escriptDataC _y_dirac=y_dirac.getDataC();
799    
800     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
801    
802     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 );
803     checkFinleyError();     checkFinleyError();
804    
805     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 );
806     checkFinleyError();     checkFinleyError();
807    
808     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 );
809       checkFinleyError();
810    
811        Finley_Assemble_PDE(mesh->Nodes,mesh->Points, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
812     checkFinleyError();     checkFinleyError();
813  }  }
814    
815  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
816                                          escript::Data& mat,                                          escript::Data& mat,
817                                          const escript::Data& D,                                          const escript::Data& D,
818                                          const escript::Data& d) const                                          const escript::Data& d,
819                                            const escript::Data& d_dirac,
820                                            const bool useHRZ) const
821  {  {
822     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
823     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
824     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
825       escriptDataC _d_dirac=d_dirac.getDataC();
826    
827     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
828    
829     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
830     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     checkFinleyError();
831      
832       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
833       checkFinleyError();
834    
835       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Points,&_mat, &_d_dirac, useHRZ);
836     checkFinleyError();     checkFinleyError();
837    
838  }  }
839    
840    
841  //  //
842  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
843  //  //
844  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
845  {  {
846     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
847    
# Line 803  void MeshAdapter::addPDEToRHS( escript:: Line 850  void MeshAdapter::addPDEToRHS( escript::
850     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
851     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
852     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
853       escriptDataC _y_dirac=y_dirac.getDataC();
854    
855     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 );
856     checkFinleyError();     checkFinleyError();
# Line 812  void MeshAdapter::addPDEToRHS( escript:: Line 860  void MeshAdapter::addPDEToRHS( escript::
860    
861     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 );
862     checkFinleyError();     checkFinleyError();
863    
864       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs , 0, 0, 0, 0, 0, &_y_dirac );
865       checkFinleyError();
866    
867  }  }
868  //  //
869  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
870  //  //
871  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
872                                             TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
873                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
874                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
875                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
876                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact,
877                                               const escript::Data& d_dirac, const escript::Data& y_dirac) const
878  {  {
879     DataArrayView::ShapeType shape;     TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
880       if (tpa==0)
881       {
882        throw FinleyAdapterException("finley only supports Paso transport problems.");
883       }
884    
885    
886       DataTypes::ShapeType shape;
887     source.expand();     source.expand();
888     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
889     escriptDataC _M=M.getDataC();     escriptDataC _M=M.getDataC();
# Line 837  void MeshAdapter::addPDEToTransportProbl Line 897  void MeshAdapter::addPDEToTransportProbl
897     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
898     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
899     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
900       escriptDataC _d_dirac=d_dirac.getDataC();
901       escriptDataC _y_dirac=y_dirac.getDataC();
902    
903    
904     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
905     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
906    
907     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 );
908     checkFinleyError();     checkFinleyError();
# Line 852  void MeshAdapter::addPDEToTransportProbl Line 915  void MeshAdapter::addPDEToTransportProbl
915    
916     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 );
917     checkFinleyError();     checkFinleyError();
918    
919       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
920       checkFinleyError();
921    
922  }  }
923    
924  //  //
# Line 859  void MeshAdapter::addPDEToTransportProbl Line 926  void MeshAdapter::addPDEToTransportProbl
926  //  //
927  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
928  {  {
929     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
930     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
931     if (inDomain!=*this)       if (inDomain!=*this)  
932        throw FinleyAdapterException("Error - Illegal domain of interpolant.");        throw FinleyAdapterException("Error - Illegal domain of interpolant.");
933     if (targetDomain!=*this)     if (targetDomain!=*this)
# Line 871  void MeshAdapter::interpolateOnDomain(es Line 938  void MeshAdapter::interpolateOnDomain(es
938     escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
939     switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
940     case(Nodes):     case(Nodes):
941     switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
942     case(Nodes):        case(Nodes):
943     case(ReducedNodes):        case(ReducedNodes):
944     case(DegreesOfFreedom):        case(DegreesOfFreedom):
945     case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
946     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());  
947        break;        break;
948     }        case(Elements):
949     break;        case(ReducedElements):
950     case(ReducedNodes):        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
951     switch(target.getFunctionSpace().getTypeCode()) {        break;
952     case(Nodes):        case(FaceElements):
953     case(ReducedNodes):        case(ReducedFaceElements):
954     case(DegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
955     case(ReducedDegreesOfFreedom):        break;
956     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(Points):
957     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
958     case(Elements):        break;
959     case(ReducedElements):        case(ContactElementsZero):
960     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        case(ReducedContactElementsZero):
961     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
962     case(FaceElements):        break;
963     case(ReducedFaceElements):        case(ContactElementsOne):
964     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        case(ReducedContactElementsOne):
965     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
966     case(Points):        break;
967     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        default:
968     break;           stringstream temp;
969     case(ContactElementsZero):           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
970     case(ReducedContactElementsZero):           throw FinleyAdapterException(temp.str());
971     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           break;
972     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());  
973        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):  
974     case(ReducedNodes):     case(ReducedNodes):
975     if (getMPISize()>1) {        switch(target.getFunctionSpace().getTypeCode()) {
976        escript::Data temp=escript::Data(in);        case(Nodes):
977        temp.expand();        case(ReducedNodes):
978        escriptDataC _in2 = temp.getDataC();        case(DegreesOfFreedom):
979        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);        case(ReducedDegreesOfFreedom):
    } else {  
980        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
981     }        break;
982     break;        case(Elements):
983     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 {  
984        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
985     }        break;
986     break;        case(FaceElements):
987     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 {  
988        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
989     }        break;
990     break;        case(Points):
    case(Points):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
    } else {  
991        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
992     }        break;
993     break;        case(ContactElementsZero):
994     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 {  
995        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());  
996        break;        break;
997     }        case(ContactElementsOne):
998     break;        case(ReducedContactElementsOne):
999     case(ReducedDegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1000     switch(target.getFunctionSpace().getTypeCode()) {        break;
1001     case(Nodes):        default:
1002     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           stringstream temp;
1003     break;           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1004     case(ReducedNodes):           throw FinleyAdapterException(temp.str());
1005     if (getMPISize()>1) {           break;
1006        escript::Data temp=escript::Data(in);        }
1007        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;  
1008     case(Elements):     case(Elements):
1009          if (target.getFunctionSpace().getTypeCode()==Elements) {
1010             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1011          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1012             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
1013          } else {
1014             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
1015          }
1016          break;
1017     case(ReducedElements):     case(ReducedElements):
1018     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1019        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1020        escriptDataC _in2 = temp.getDataC();        } else {
1021        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
1022     } else {        }
1023        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        break;
    }  
    break;  
1024     case(FaceElements):     case(FaceElements):
1025          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
1026             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1027          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1028             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1029          } else {
1030             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1031          }
1032          break;
1033     case(ReducedFaceElements):     case(ReducedFaceElements):
1034     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1035        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1036        escriptDataC _in2 = temp.getDataC();        } else {
1037        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
1038     } else {        }
1039        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
    }  
    break;  
1040     case(Points):     case(Points):
1041     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==Points) {
1042        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1043        escriptDataC _in2 = temp.getDataC();        } else {
1044        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1045     } else {        }
1046        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        break;
    }  
    break;  
1047     case(ContactElementsZero):     case(ContactElementsZero):
1048     case(ContactElementsOne):     case(ContactElementsOne):
1049          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1050             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1051          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1052             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1053          } else {
1054             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1055          }
1056          break;
1057     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1058     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1059     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1060        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1061        escriptDataC _in2 = temp.getDataC();        } else {
1062        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1063     } else {        }
1064        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        break;
1065     }     case(DegreesOfFreedom):      
1066     break;        switch(target.getFunctionSpace().getTypeCode()) {
1067     default:        case(ReducedDegreesOfFreedom):
1068        stringstream temp;        case(DegreesOfFreedom):
1069        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1070        throw FinleyAdapterException(temp.str());        break;
1071      
1072          case(Nodes):
1073          case(ReducedNodes):
1074          if (getMPISize()>1) {
1075             escript::Data temp=escript::Data(in);
1076             temp.expand();
1077             escriptDataC _in2 = temp.getDataC();
1078             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1079          } else {
1080             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1081          }
1082          break;
1083          case(Elements):
1084          case(ReducedElements):
1085          if (getMPISize()>1) {
1086             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1087             escriptDataC _in2 = temp.getDataC();
1088             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1089          } else {
1090             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1091          }
1092          break;
1093          case(FaceElements):
1094          case(ReducedFaceElements):
1095          if (getMPISize()>1) {
1096             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1097             escriptDataC _in2 = temp.getDataC();
1098             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1099      
1100          } else {
1101             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1102          }
1103          break;
1104          case(Points):
1105          if (getMPISize()>1) {
1106             //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1107             //escriptDataC _in2 = temp.getDataC();
1108          } else {
1109             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1110          }
1111          break;
1112          case(ContactElementsZero):
1113          case(ContactElementsOne):
1114          case(ReducedContactElementsZero):
1115          case(ReducedContactElementsOne):
1116          if (getMPISize()>1) {
1117             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1118             escriptDataC _in2 = temp.getDataC();
1119             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1120          } else {
1121             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1122          }
1123          break;
1124          default:
1125             stringstream temp;
1126             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1127             throw FinleyAdapterException(temp.str());
1128             break;
1129          }
1130          break;
1131       case(ReducedDegreesOfFreedom):
1132          switch(target.getFunctionSpace().getTypeCode()) {
1133          case(Nodes):
1134          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1135          break;
1136          case(ReducedNodes):
1137          if (getMPISize()>1) {
1138             escript::Data temp=escript::Data(in);
1139             temp.expand();
1140             escriptDataC _in2 = temp.getDataC();
1141             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1142          } else {
1143             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1144          }
1145          break;
1146          case(DegreesOfFreedom):
1147          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1148          break;
1149          case(ReducedDegreesOfFreedom):
1150          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1151          break;
1152          case(Elements):
1153          case(ReducedElements):
1154          if (getMPISize()>1) {
1155             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1156             escriptDataC _in2 = temp.getDataC();
1157             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1158          } else {
1159             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1160          }
1161          break;
1162          case(FaceElements):
1163          case(ReducedFaceElements):
1164          if (getMPISize()>1) {
1165             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1166             escriptDataC _in2 = temp.getDataC();
1167             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1168          } else {
1169             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1170          }
1171          break;
1172          case(Points):
1173          if (getMPISize()>1) {
1174             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1175             escriptDataC _in2 = temp.getDataC();
1176             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1177          } else {
1178             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1179          }
1180          break;
1181          case(ContactElementsZero):
1182          case(ContactElementsOne):
1183          case(ReducedContactElementsZero):
1184          case(ReducedContactElementsOne):
1185          if (getMPISize()>1) {
1186             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1187             escriptDataC _in2 = temp.getDataC();
1188             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1189          } else {
1190             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1191          }
1192          break;
1193          default:
1194             stringstream temp;
1195             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1196             throw FinleyAdapterException(temp.str());
1197             break;
1198          }
1199        break;        break;
    }  
    break;  
1200     default:     default:
1201        stringstream temp;        stringstream temp;
1202        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 1144  void MeshAdapter::interpolateOnDomain(es Line 1211  void MeshAdapter::interpolateOnDomain(es
1211  //  //
1212  void MeshAdapter::setToX(escript::Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
1213  {  {
1214     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1215     if (argDomain!=*this)     if (argDomain!=*this)
1216        throw FinleyAdapterException("Error - Illegal domain of data point locations");        throw FinleyAdapterException("Error - Illegal domain of data point locations");
1217     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1153  void MeshAdapter::setToX(escript::Data& Line 1220  void MeshAdapter::setToX(escript::Data&
1220        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1221        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1222     } else {     } else {
1223        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1224        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1225        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1226        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1167  void MeshAdapter::setToX(escript::Data& Line 1234  void MeshAdapter::setToX(escript::Data&
1234  //  //
1235  void MeshAdapter::setToNormal(escript::Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
1236  {  {
1237     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1238       const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1239     if (normalDomain!=*this)     if (normalDomain!=*this)
1240        throw FinleyAdapterException("Error - Illegal domain of normal locations");        throw FinleyAdapterException("Error - Illegal domain of normal locations");
1241     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1222  void MeshAdapter::setToNormal(escript::D Line 1290  void MeshAdapter::setToNormal(escript::D
1290  //  //
1291  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1292  {  {
1293     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1294     if (targetDomain!=*this)     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1295       if (targetDomain!=this)
1296        throw FinleyAdapterException("Error - Illegal domain of interpolation target");        throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1297    
1298     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
# Line 1232  void MeshAdapter::interpolateACross(escr Line 1301  void MeshAdapter::interpolateACross(escr
1301  //  //
1302  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
1303  //  //
1304  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1305  {  {
1306     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1307     if (argDomain!=*this)     if (argDomain!=*this)
1308        throw FinleyAdapterException("Error - Illegal domain of integration kernel");        throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1309    
# Line 1245  void MeshAdapter::setToIntegrals(std::ve Line 1314  void MeshAdapter::setToIntegrals(std::ve
1314     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1315     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1316     case(Nodes):     case(Nodes):
1317     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1318     _temp=temp.getDataC();     _temp=temp.getDataC();
1319     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1320     break;     break;
1321     case(ReducedNodes):     case(ReducedNodes):
1322     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1323     _temp=temp.getDataC();     _temp=temp.getDataC();
1324     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1325     break;     break;
# Line 1282  void MeshAdapter::setToIntegrals(std::ve Line 1351  void MeshAdapter::setToIntegrals(std::ve
1351     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1352     break;     break;
1353     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1354     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1355     _temp=temp.getDataC();     _temp=temp.getDataC();
1356     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1357     break;     break;
1358     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1359     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1360     _temp=temp.getDataC();     _temp=temp.getDataC();
1361     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1362     break;     break;
# Line 1306  void MeshAdapter::setToIntegrals(std::ve Line 1375  void MeshAdapter::setToIntegrals(std::ve
1375  //  //
1376  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1377  {  {
1378     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1379     if (argDomain!=*this)     if (argDomain!=*this)
1380        throw FinleyAdapterException("Error - Illegal domain of gradient argument");        throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1381     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1382     if (gradDomain!=*this)     if (gradDomain!=*this)
1383        throw FinleyAdapterException("Error - Illegal domain of gradient");        throw FinleyAdapterException("Error - Illegal domain of gradient");
1384    
# Line 1319  void MeshAdapter::setToGradient(escript: Line 1388  void MeshAdapter::setToGradient(escript:
1388     escript::Data temp;     escript::Data temp;
1389     if (getMPISize()>1) {     if (getMPISize()>1) {
1390        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1391           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1392           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1393        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1394           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1395           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1396        } else {        } else {
1397           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1431  void MeshAdapter::setToSize(escript::Dat Line 1500  void MeshAdapter::setToSize(escript::Dat
1500     checkFinleyError();     checkFinleyError();
1501  }  }
1502    
1503  // sets the location of nodes:  //
1504    // sets the location of nodes
1505    //
1506  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1507  {  {
1508     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1509     escriptDataC tmp;     escriptDataC tmp;
1510     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1511     if (newDomain!=*this)     if (newDomain!=*this)
1512        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1513     tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1514     Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1515           Finley_Mesh_setCoordinates(mesh,&tmp);
1516       } else {
1517           throw FinleyAdapterException("As of version escript3.3 - SetNewX only accepts ContinuousFunction arguments please interpolate.");    
1518    /*       escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1519           tmp = new_x_inter.getDataC();
1520           Finley_Mesh_setCoordinates(mesh,&tmp);*/
1521       }
1522     checkFinleyError();     checkFinleyError();
1523  }  }
1524    
1525  // saves a data array in openDX format:  //
1526  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  // Helper for the save* methods. Extracts optional data variable names and the
1527    // corresponding pointers from python dictionary. Caller must free arrays.
1528    //
1529    void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1530  {  {
1531     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;     numData = boost::python::extract<int>(arg.attr("__len__")());
    const int num_data=boost::python::extract<int>(arg.attr("__len__")());  
1532     /* win32 refactor */     /* win32 refactor */
1533     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1534     for(int i=0;i<num_data;i++)     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1535     {     dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
       names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
    }  
   
    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  
    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  
1536    
1537     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1538     for (int i=0;i<num_data;++i) {     for (int i=0; i<numData; ++i) {
1539        std::string n=boost::python::extract<std::string>(keys[i]);        string n=boost::python::extract<string>(keys[i]);
1540        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1541        if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1542           throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");           throw FinleyAdapterException("Error: Data must be defined on same Domain");
1543        data[i]=d.getDataC();        data[i] = d.getDataC();
1544        ptr_data[i]= &(data[i]);        dataPtr[i] = &(data[i]);
1545        if (n.length()>MAX_namelength-1) {        names[i] = TMPMEMALLOC(n.length()+1, char);
1546           strncpy(names[i],n.c_str(),MAX_namelength-1);        strcpy(names[i], n.c_str());
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
1547     }     }
1548     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  }
1549    
1550    //
1551    // saves mesh and optionally data arrays in openDX format
1552    //
1553    void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1554    {
1555       int num_data;
1556       char **names;
1557       escriptDataC *data;
1558       escriptDataC **ptr_data;
1559    
1560       extractArgsFromDict(arg, num_data, names, data, ptr_data);
1561       Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1562     checkFinleyError();     checkFinleyError();
1563    
1564     /* win32 refactor */     /* win32 refactor */
1565     TMPMEMFREE(data);     TMPMEMFREE(data);
1566     TMPMEMFREE(ptr_data);     TMPMEMFREE(ptr_data);
1567     for(int i=0;i<num_data;i++)     for(int i=0; i<num_data; i++)
1568     {     {
1569        TMPMEMFREE(names[i]);        TMPMEMFREE(names[i]);
1570     }     }
1571     TMPMEMFREE(names);     TMPMEMFREE(names);
   
    return;  
1572  }  }
1573    
1574  // saves a data array in openVTK format:  //
1575  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  // saves mesh and optionally data arrays in VTK format
1576    //
1577    void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1578  {  {
1579     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1580     const int num_data=boost::python::extract<int>(arg.attr("__len__")());      pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1581     /* win32 refactor */                metadata, metadata_schema, arg);
1582     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;  }
1583     for(int i=0;i<num_data;i++)  
1584     {  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1585        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  {
1586     }      if (getMPISize() > 1) {
1587    #ifdef ESYS_MPI
1588     escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;          index_t myFirstNode=0, myLastNode=0, k=0;
1589     escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;          index_t* globalNodeIndex=0;
1590            Finley_Mesh* mesh_p=m_finleyMesh.get();
1591            /*
1592             * this method is only used by saveDataCSV which would use the returned
1593             * values for reduced nodes wrongly so this case is disabled for now
1594            if (fs_code == FINLEY_REDUCED_NODES)
1595            {
1596                myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1597                myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1598                globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1599            }
1600            else
1601            */
1602            if (fs_code == FINLEY_NODES)
1603            {
1604                myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1605                myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1606                globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1607            }
1608            else
1609            {
1610                throw FinleyAdapterException("Unsupported function space type for ownSample()");
1611            }
1612    
1613     boost::python::list keys=arg.keys();          k=globalNodeIndex[id];
1614     for (int i=0;i<num_data;++i) {          return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1615        std::string n=boost::python::extract<std::string>(keys[i]);  #endif
1616        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);      }
1617        if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)      return true;
1618           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);  
1619    
    checkFinleyError();  
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0;i<num_data;i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
1620    
1621     return;  //
1622  }  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1623                                                                                                                                                                        //
1624                                                                                                                                                                        ASM_ptr MeshAdapter::newSystemMatrix(
 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  
 SystemMatrixAdapter MeshAdapter::newSystemMatrix(  
1625                                                   const int row_blocksize,                                                   const int row_blocksize,
1626                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1627                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1546  SystemMatrixAdapter MeshAdapter::newSyst Line 1631  SystemMatrixAdapter MeshAdapter::newSyst
1631     int reduceRowOrder=0;     int reduceRowOrder=0;
1632     int reduceColOrder=0;     int reduceColOrder=0;
1633     // is the domain right?     // is the domain right?
1634     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1635     if (row_domain!=*this)     if (row_domain!=*this)
1636        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.");
1637     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1638     if (col_domain!=*this)     if (col_domain!=*this)
1639        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1640     // is the function space type right     // is the function space type right
# Line 1579  SystemMatrixAdapter MeshAdapter::newSyst Line 1664  SystemMatrixAdapter MeshAdapter::newSyst
1664  #endif  #endif
1665     }     }
1666     else {     else {
1667        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1668     }     }
1669     checkPasoError();     checkPasoError();
1670     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1671     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1672       return ASM_ptr(sma);
1673    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1674  }  }
1675    
1676    //
1677  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1678  TransportProblemAdapter MeshAdapter::newTransportProblem(  //
1679                                                           const double theta,  ATP_ptr MeshAdapter::newTransportProblem(
1680                                                           const int blocksize,                                                           const int blocksize,
1681                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1682                                                           const int type) const                                                           const int type) const
1683  {  {
1684     int reduceOrder=0;     int reduceOrder=0;
1685     // is the domain right?     // is the domain right?
1686     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(functionspace.getDomain());     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1687     if (domain!=*this)     if (domain!=*this)
1688        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1689     // is the function space type right     // is the function space type right
# Line 1609  TransportProblemAdapter MeshAdapter::new Line 1698  TransportProblemAdapter MeshAdapter::new
1698    
1699     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1700     checkFinleyError();     checkFinleyError();
1701     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1702     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1703     checkPasoError();     checkPasoError();
1704     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1705     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1706       return ATP_ptr(tpa);
1707    //   return TransportProblemAdapter(transportProblem, blocksize,functionspace);
1708  }  }
1709    
1710  //  //
# Line 1652  bool MeshAdapter::isCellOriented(int fun Line 1743  bool MeshAdapter::isCellOriented(int fun
1743     return false;     return false;
1744  }  }
1745    
1746    bool
1747    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1748    {
1749       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1750        class 1: DOF <-> Nodes
1751        class 2: ReducedDOF <-> ReducedNodes
1752        class 3: Points
1753        class 4: Elements
1754        class 5: ReducedElements
1755        class 6: FaceElements
1756        class 7: ReducedFaceElements
1757        class 8: ContactElementZero <-> ContactElementOne
1758        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1759    
1760       There is also a set of lines. Interpolation is possible down a line but not between lines.
1761       class 1 and 2 belong to all lines so aren't considered.
1762        line 0: class 3
1763        line 1: class 4,5
1764        line 2: class 6,7
1765        line 3: class 8,9
1766    
1767       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1768       eg hasnodes is true if we have at least one instance of Nodes.
1769       */
1770        if (fs.empty())
1771        {
1772            return false;
1773        }
1774        vector<int> hasclass(10);
1775        vector<int> hasline(4);
1776        bool hasnodes=false;
1777        bool hasrednodes=false;
1778        bool hascez=false;
1779        bool hasrcez=false;
1780        for (int i=0;i<fs.size();++i)
1781        {
1782        switch(fs[i])
1783        {
1784        case(Nodes):   hasnodes=true;   // no break is deliberate
1785        case(DegreesOfFreedom):
1786            hasclass[1]=1;
1787            break;
1788        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1789        case(ReducedDegreesOfFreedom):
1790            hasclass[2]=1;
1791            break;
1792        case(Points):
1793            hasline[0]=1;
1794            hasclass[3]=1;
1795            break;
1796        case(Elements):
1797            hasclass[4]=1;
1798            hasline[1]=1;
1799            break;
1800        case(ReducedElements):
1801            hasclass[5]=1;
1802            hasline[1]=1;
1803            break;
1804        case(FaceElements):
1805            hasclass[6]=1;
1806            hasline[2]=1;
1807            break;
1808        case(ReducedFaceElements):
1809            hasclass[7]=1;
1810            hasline[2]=1;
1811            break;
1812        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1813        case(ContactElementsOne):
1814            hasclass[8]=1;
1815            hasline[3]=1;
1816            break;
1817        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1818        case(ReducedContactElementsOne):
1819            hasclass[9]=1;
1820            hasline[3]=1;
1821            break;
1822        default:
1823            return false;
1824        }
1825        }
1826        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1827        // fail if we have more than one leaf group
1828    
1829        if (totlines>1)
1830        {
1831        return false;   // there are at least two branches we can't interpolate between
1832        }
1833        else if (totlines==1)
1834        {
1835        if (hasline[0]==1)      // we have points
1836        {
1837            resultcode=Points;
1838        }
1839        else if (hasline[1]==1)
1840        {
1841            if (hasclass[5]==1)
1842            {
1843            resultcode=ReducedElements;
1844            }
1845            else
1846            {
1847            resultcode=Elements;
1848            }
1849        }
1850        else if (hasline[2]==1)
1851        {
1852            if (hasclass[7]==1)
1853            {
1854            resultcode=ReducedFaceElements;
1855            }
1856            else
1857            {
1858            resultcode=FaceElements;
1859            }
1860        }
1861        else    // so we must be in line3
1862        {
1863            if (hasclass[9]==1)
1864            {
1865            // need something from class 9
1866            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1867            }
1868            else
1869            {
1870            // something from class 8
1871            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1872            }
1873        }
1874        }
1875        else    // totlines==0
1876        {
1877        if (hasclass[2]==1)
1878        {
1879            // something from class 2
1880            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1881        }
1882        else
1883        {   // something from class 1
1884            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1885        }
1886        }
1887        return true;
1888    }
1889    
1890  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1891  {  {
1892     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1893     case(Nodes):     case(Nodes):
1894     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1895     case(Nodes):      case(Nodes):
1896     case(ReducedNodes):      case(ReducedNodes):
1897     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1898     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1899     case(Elements):      case(Elements):
1900     case(ReducedElements):      case(ReducedElements):
1901     case(FaceElements):      case(FaceElements):
1902     case(ReducedFaceElements):      case(ReducedFaceElements):
1903     case(Points):      case(Points):
1904     case(ContactElementsZero):      case(ContactElementsZero):
1905     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1906     case(ContactElementsOne):      case(ContactElementsOne):
1907     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1908     return true;      return true;
1909     default:      default:
1910        stringstream temp;            stringstream temp;
1911        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;
1912        throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1913     }     }
1914     break;     break;
1915     case(ReducedNodes):     case(ReducedNodes):
1916     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1917     case(ReducedNodes):      case(ReducedNodes):
1918     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1919     case(Elements):      case(Elements):
1920     case(ReducedElements):      case(ReducedElements):
1921     case(FaceElements):      case(FaceElements):
1922     case(ReducedFaceElements):      case(ReducedFaceElements):
1923     case(Points):      case(Points):
1924     case(ContactElementsZero):      case(ContactElementsZero):
1925     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1926     case(ContactElementsOne):      case(ContactElementsOne):
1927     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1928     return true;      return true;
1929     case(Nodes):      case(Nodes):
1930     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1931     return false;      return false;
1932     default:      default:
1933        stringstream temp;          stringstream temp;
1934        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;
1935        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1936     }     }
1937     break;     break;
1938     case(Elements):     case(Elements):
1939     if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1940        return true;        return true;
1941     } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1942        return true;        return true;
1943     } else {          } else {
1944        return false;            return false;
1945     }          }
1946     case(ReducedElements):     case(ReducedElements):
1947     if (functionSpaceType_target==ReducedElements) {      if (functionSpaceType_target==ReducedElements) {
1948        return true;        return true;
1949     } else {      } else {
1950        return false;            return false;
1951     }      }
1952     case(FaceElements):     case(FaceElements):
1953     if (functionSpaceType_target==FaceElements) {      if (functionSpaceType_target==FaceElements) {
1954        return true;              return true;
1955     } else if (functionSpaceType_target==ReducedFaceElements) {      } else if (functionSpaceType_target==ReducedFaceElements) {
1956        return true;              return true;
1957     } else {      } else {
1958        return false;              return false;
1959     }      }
1960     case(ReducedFaceElements):     case(ReducedFaceElements):
1961     if (functionSpaceType_target==ReducedFaceElements) {      if (functionSpaceType_target==ReducedFaceElements) {
1962        return true;              return true;
1963     } else {      } else {
1964        return false;          return false;
1965     }      }
1966     case(Points):     case(Points):
1967     if (functionSpaceType_target==Points) {      if (functionSpaceType_target==Points) {
1968        return true;              return true;
1969     } else {      } else {
1970        return false;              return false;
1971     }      }
1972     case(ContactElementsZero):     case(ContactElementsZero):
1973     case(ContactElementsOne):     case(ContactElementsOne):
1974     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1975        return true;              return true;
1976     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1977        return true;              return true;
1978     } else {      } else {
1979        return false;              return false;
1980     }      }
1981     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1982     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1983     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1984        return true;              return true;
1985     } else {      } else {
1986        return false;              return false;
1987     }      }
1988     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1989     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1990     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1991     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1992     case(Nodes):      case(Nodes):
1993     case(ReducedNodes):      case(ReducedNodes):
1994     case(Elements):      case(Elements):
1995     case(ReducedElements):      case(ReducedElements):
1996     case(Points):      case(Points):
1997     case(FaceElements):      case(FaceElements):
1998     case(ReducedFaceElements):      case(ReducedFaceElements):
1999     case(ContactElementsZero):      case(ContactElementsZero):
2000     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
2001     case(ContactElementsOne):      case(ContactElementsOne):
2002     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
2003     return true;      return true;
2004     default:      default:
2005        stringstream temp;          stringstream temp;
2006        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;
2007        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
2008     }      }
2009     break;      break;
2010     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2011     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
2012     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
2013     case(ReducedNodes):      case(ReducedNodes):
2014     case(Elements):      case(Elements):
2015     case(ReducedElements):      case(ReducedElements):
2016     case(FaceElements):      case(FaceElements):
2017     case(ReducedFaceElements):      case(ReducedFaceElements):
2018     case(Points):      case(Points):
2019     case(ContactElementsZero):      case(ContactElementsZero):
2020     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
2021     case(ContactElementsOne):      case(ContactElementsOne):
2022     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
2023     return true;      return true;
2024     case(Nodes):      case(Nodes):
2025     case(DegreesOfFreedom):      case(DegreesOfFreedom):
2026     return false;      return false;
2027     default:      default:
2028        stringstream temp;          stringstream temp;
2029        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;
2030        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
2031     }      }
2032     break;      break;
2033     default:     default:
2034        stringstream temp;        stringstream temp;
2035        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 1825  bool MeshAdapter::operator!=(const Abstr Line 2060  bool MeshAdapter::operator!=(const Abstr
2060     return !(operator==(other));     return !(operator==(other));
2061  }  }
2062    
2063  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2064  {  {
2065     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2066     checkPasoError();     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2067     return out;             package, symmetry, mesh->MPIInfo);
2068    }
2069    
2070    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2071    {
2072       Finley_Mesh* mesh=m_finleyMesh.get();
2073       return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2074               package, symmetry, mesh->MPIInfo);
2075  }  }
2076    
2077  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2078  {  {
2079     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2080  }  }
2081    
2082  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2083  {  {
2084     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2085  }  }
2086    
2087  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2088  {  {
2089     return function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2090  }  }
2091    
2092  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2093  {  {
2094     int *out = NULL;     int *out = NULL;
2095     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2006  void MeshAdapter::setTags(const int func Line 2248  void MeshAdapter::setTags(const int func
2248     return;     return;
2249  }  }
2250    
2251  void MeshAdapter::setTagMap(const std::string& name,  int tag)  void MeshAdapter::setTagMap(const string& name,  int tag)
2252  {  {
2253     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2254     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2255     checkPasoError();     checkFinleyError();
2256     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2257  }  }
2258    
2259  int MeshAdapter::getTag(const std::string& name) const  int MeshAdapter::getTag(const string& name) const
2260  {  {
2261     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2262     int tag=0;     int tag=0;
2263     tag=Finley_Mesh_getTag(mesh, name.c_str());     tag=Finley_Mesh_getTag(mesh, name.c_str());
2264     checkPasoError();     checkFinleyError();
2265     // throwStandardException("MeshAdapter::getTag is not implemented.");     // throwStandardException("MeshAdapter::getTag is not implemented.");
2266     return tag;     return tag;
2267  }  }
2268    
2269  bool MeshAdapter::isValidTagName(const std::string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2270  {  {
2271     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2272     return Finley_Mesh_isValidTagName(mesh,name.c_str());     return Finley_Mesh_isValidTagName(mesh,name.c_str());
2273  }  }
2274    
2275  std::string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2276  {  {
2277     stringstream temp;     stringstream temp;
2278     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2084  int MeshAdapter::getNumberOfTagsInUse(in Line 2326  int MeshAdapter::getNumberOfTagsInUse(in
2326    }    }
2327    return numTags;    return numTags;
2328  }  }
2329  int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  
2330    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2331  {  {
2332    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2333    index_t* tags=NULL;    index_t* tags=NULL;
# Line 2127  int* MeshAdapter::borrowListOfTagsInUse( Line 2370  int* MeshAdapter::borrowListOfTagsInUse(
2370  }  }
2371    
2372    
2373    bool MeshAdapter::canTag(int functionSpaceCode) const
2374    {
2375      switch(functionSpaceCode) {
2376       case(Nodes):
2377       case(Elements):
2378       case(ReducedElements):
2379       case(FaceElements):
2380       case(ReducedFaceElements):
2381       case(Points):
2382       case(ContactElementsZero):
2383       case(ReducedContactElementsZero):
2384       case(ContactElementsOne):
2385       case(ReducedContactElementsOne):
2386              return true;
2387       case(ReducedNodes):
2388       case(DegreesOfFreedom):
2389       case(ReducedDegreesOfFreedom):
2390          return false;
2391       default:
2392        return false;
2393      }
2394    }
2395    
2396    AbstractDomain::StatusType MeshAdapter::getStatus() const
2397    {
2398      Finley_Mesh* mesh=m_finleyMesh.get();
2399      return Finley_Mesh_getStatus(mesh);
2400    }
2401    
2402    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2403    {
2404      
2405      Finley_Mesh* mesh=m_finleyMesh.get();
2406      int order =-1;
2407      switch(functionSpaceCode) {
2408       case(Nodes):
2409       case(DegreesOfFreedom):
2410              order=mesh->approximationOrder;
2411              break;
2412       case(ReducedNodes):
2413       case(ReducedDegreesOfFreedom):
2414              order=mesh->reducedApproximationOrder;
2415              break;
2416       case(Elements):
2417       case(FaceElements):
2418       case(Points):
2419       case(ContactElementsZero):
2420       case(ContactElementsOne):
2421              order=mesh->integrationOrder;
2422              break;
2423       case(ReducedElements):
2424       case(ReducedFaceElements):
2425       case(ReducedContactElementsZero):
2426       case(ReducedContactElementsOne):
2427              order=mesh->reducedIntegrationOrder;
2428              break;
2429       default:
2430          stringstream temp;
2431          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2432          throw FinleyAdapterException(temp.str());
2433      }
2434      return order;
2435    }
2436    
2437    bool MeshAdapter::supportsContactElements() const
2438    {
2439      return true;
2440    }
2441    
2442    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2443    {
2444      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2445    }
2446    
2447    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2448    {
2449      Finley_ReferenceElementSet_dealloc(m_refSet);
2450    }
2451    
2452    // points will be flattened
2453    void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2454    {
2455          const int dim = getDim();
2456          int numPoints=points.size()/dim;
2457          int numTags=tags.size();
2458          Finley_Mesh* mesh=m_finleyMesh.get();
2459          
2460          if ( points.size() % dim != 0 )
2461          {
2462        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2463          }
2464          
2465          if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2466         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2467          
2468          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2469          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2470          
2471          for (int i=0;i<numPoints;++i) {
2472           points_ptr[ i * dim     ] = points[i * dim ];
2473           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2474           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2475               tags_ptr[i]=tags[i];
2476          }
2477          
2478          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2479          checkFinleyError();
2480          
2481          TMPMEMFREE(points_ptr);
2482          TMPMEMFREE(tags_ptr);
2483    }
2484    
2485    
2486    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2487    // {
2488    //       const int dim = getDim();
2489    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2490    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2491    //       Finley_Mesh* mesh=m_finleyMesh.get();
2492    //      
2493    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2494    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2495    //      
2496    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2497    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2498    //      
2499    //       for (int i=0;i<numPoints;++i) {
2500    //     int tag_id=-1;
2501    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2502    //     if  ( numComps !=   dim ) {
2503    //                stringstream temp;            
2504    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2505    //                throw FinleyAdapterException(temp.str());
2506    //     }
2507    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2508    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2509    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2510    //    
2511    //     if ( numTags > 0) {
2512    //            boost::python::extract<string> ex_str(tags[i]);
2513    //        if  ( ex_str.check() ) {
2514    //            tag_id=getTag( ex_str());
2515    //        } else {
2516    //             boost::python::extract<int> ex_int(tags[i]);
2517    //             if ( ex_int.check() ) {
2518    //                 tag_id=ex_int();
2519    //             } else {
2520    //              stringstream temp;          
2521    //                  temp << "Error - unable to extract tag for point " << i;
2522    //              throw FinleyAdapterException(temp.str());
2523    //            }
2524    //        }
2525    //     }      
2526    //            tags_ptr[i]=tag_id;
2527    //       }
2528    //      
2529    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2530    //       checkPasoError();
2531    //      
2532    //       TMPMEMFREE(points_ptr);
2533    //       TMPMEMFREE(tags_ptr);
2534    // }
2535    
2536    /*
2537    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2538    {  
2539        boost::python::list points =  boost::python::list();
2540        boost::python::list tags =  boost::python::list();
2541        points.append(point);
2542        tags.append(tag);
2543        addDiracPoints(points, tags);
2544    }
2545    */
2546    
2547    /*
2548    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2549    {
2550            boost::python::list points =   boost::python::list();
2551            boost::python::list tags =   boost::python::list();
2552            points.append(point);
2553            tags.append(tag);
2554            addDiracPoints(points, tags);
2555    }
2556    */
2557  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1748  
changed lines
  Added in v.3998

  ViewVC Help
Powered by ViewVC 1.1.26