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

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

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

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

Legend:
Removed from v.1632  
changed lines
  Added in v.3793

  ViewVC Help
Powered by ViewVC 1.1.26