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

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

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

revision 1346 by ksteube, Wed Nov 14 22:48:12 2007 UTC revision 3736 by jfenwick, Mon Dec 12 01:26:51 2011 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  }  }
25  #include <vector>  
26    #include <boost/python/import.hpp>
27    
28  using namespace std;  using namespace std;
29  using namespace escript;  using namespace escript;
30    using namespace paso;
31    
32  namespace finley {  namespace finley {
33    
# Line 48  const int MeshAdapter::ReducedContactEle Line 50  const int MeshAdapter::ReducedContactEle
50    
51  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
52  {  {
53    setFunctionSpaceTypeNames();     setFunctionSpaceTypeNames();
54    //     //
55    // need to use a null_deleter as Finley_Mesh_free deletes the pointer     // need to use a null_deleter as Finley_Mesh_free deletes the pointer
56    // for us.     // for us.
57    m_finleyMesh.reset(finleyMesh,null_deleter());     m_finleyMesh.reset(finleyMesh,null_deleter());
58  }  }
59    
60  //  //
# Line 60  MeshAdapter::MeshAdapter(Finley_Mesh* fi Line 62  MeshAdapter::MeshAdapter(Finley_Mesh* fi
62  MeshAdapter::MeshAdapter(const MeshAdapter& in):  MeshAdapter::MeshAdapter(const MeshAdapter& in):
63  m_finleyMesh(in.m_finleyMesh)  m_finleyMesh(in.m_finleyMesh)
64  {  {
65    setFunctionSpaceTypeNames();     setFunctionSpaceTypeNames();
66  }  }
67    
68  MeshAdapter::~MeshAdapter()  MeshAdapter::~MeshAdapter()
69  {  {
70    //     //
71    // I hope the case for the pointer being zero has been taken care of.     // I hope the case for the pointer being zero has been taken care of.
72    //  cout << "In MeshAdapter destructor." << endl;     //  cout << "In MeshAdapter destructor." << endl;
73    if (m_finleyMesh.unique()) {     if (m_finleyMesh.unique()) {
74      Finley_Mesh_free(m_finleyMesh.get());        Finley_Mesh_free(m_finleyMesh.get());
75    }     }
76  }  }
77    
78  int MeshAdapter::getMPISize() const  int MeshAdapter::getMPISize() const
# Line 81  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());
122    Finley_Mesh_write(m_finleyMesh.get(),fName);     Finley_Mesh_write(m_finleyMesh.get(),fName);
123    checkFinleyError();     checkFinleyError();
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[25]; /* ksteube how big should ncdims be? */     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
136     NcVar *ids, *data;     NcVar *ids;
137     int *int_ptr;     int *int_ptr;
138     Finley_Mesh *mesh = m_finleyMesh.get();     Finley_Mesh *mesh = m_finleyMesh.get();
139       Finley_TagMap* tag_map;
140       int num_Tags = 0;
141     int mpi_size             = mesh->MPIInfo->size;     int mpi_size             = mesh->MPIInfo->size;
142     int mpi_rank             = mesh->MPIInfo->rank;     int mpi_rank             = mesh->MPIInfo->rank;
143     int numDim               = mesh->Nodes->numDim;     int numDim               = mesh->Nodes->numDim;
# Line 119  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     char *newFileName = Paso_MPI_appendRankToFileName(strdup(fileName.c_str()), mpi_size, mpi_rank);  #ifdef ESYS_MPI
153       MPI_Status status;
154    #endif
155    
156    /* Incoming token indicates it's my turn to write */
157    #ifdef ESYS_MPI
158       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 = Esys_MPI_appendRankToFileName(fileName.c_str(),
162                                                         mpi_size, mpi_rank);
163    
164       /* Figure out how much storage is required for tags */
165       tag_map = mesh->TagMap;
166       num_Tags = 0;
167       while (tag_map) {
168          num_Tags++;
169          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)
211          if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
212             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->referenceElementSet->referenceElement->Type->TypeId) )
244          throw DataException(msgPrefix+"add_att(Elements_TypeId)");
245       if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
246          throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
247       if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
248          throw DataException(msgPrefix+"add_att(ContactElements_TypeId)");
249       if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
250          throw DataException(msgPrefix+"add_att(Points_TypeId)");
251       if (!dataFile.add_att("num_Tags", num_Tags) )
252          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 (numNodes>0) {     if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
258          throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
259     // Nodes Id     int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
260     if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )     if (! (ids->put(int_ptr, mpi_size+1)) )
261          throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
    int_ptr = &mesh->Nodes->Id[0];  
    if (! (ids->put(int_ptr, numNodes)) )  
         throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);  
   
    // Nodes Tag  
    if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )  
         throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);  
    int_ptr = &mesh->Nodes->Tag[0];  
    if (! (ids->put(int_ptr, numNodes)) )  
         throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);  
   
    // Nodes gDOF  
    if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )  
         throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);  
    int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];  
    if (! (ids->put(int_ptr, numNodes)) )  
         throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);  
   
    // Nodes global node index  
    if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )  
         throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);  
    int_ptr = &mesh->Nodes->globalNodesIndex[0];  
    if (! (ids->put(int_ptr, numNodes)) )  
         throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);  
   
    // Nodes grDof  
    if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )  
         throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);  
    int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];  
    if (! (ids->put(int_ptr, numNodes)) )  
         throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);  
   
    // Nodes grNI  
    if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )  
         throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);  
    int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];  
    if (! (ids->put(int_ptr, numNodes)) )  
         throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);  
   
    // Nodes Coordinates  
    if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )  
         throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);  
    if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )  
         throw DataException("Error - MeshAdapter::dump: copy Nodes_Coordinates to netCDF buffer failed: " + *newFileName);  
262    
263     // Nodes degreesOfFreedomDistribution     // Nodes degreesOfFreedomDistribution
264     if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )     if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
265          throw DataException("Error - MeshAdapter::dump: appending Nodes_DofDistribution to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
266     int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];     int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
267     if (! (ids->put(int_ptr, mpi_size+1)) )     if (! (ids->put(int_ptr, mpi_size+1)) )
268          throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);        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) {
273    
274          // Nodes Id
275          if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
276             throw DataException(msgPrefix+"add_var(Nodes_Id)");
277          int_ptr = &mesh->Nodes->Id[0];
278          if (! (ids->put(int_ptr, numNodes)) )
279             throw DataException(msgPrefix+"put(Nodes_Id)");
280    
281          // Nodes Tag
282          if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
283             throw DataException(msgPrefix+"add_var(Nodes_Tag)");
284          int_ptr = &mesh->Nodes->Tag[0];
285          if (! (ids->put(int_ptr, numNodes)) )
286             throw DataException(msgPrefix+"put(Nodes_Tag)");
287    
288          // Nodes gDOF
289          if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
290             throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
291          int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
292          if (! (ids->put(int_ptr, numNodes)) )
293             throw DataException(msgPrefix+"put(Nodes_gDOF)");
294    
295          // Nodes global node index
296          if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
297             throw DataException(msgPrefix+"add_var(Nodes_gNI)");
298          int_ptr = &mesh->Nodes->globalNodesIndex[0];
299          if (! (ids->put(int_ptr, numNodes)) )
300             throw DataException(msgPrefix+"put(Nodes_gNI)");
301    
302          // Nodes grDof
303          if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
304             throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
305          int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
306          if (! (ids->put(int_ptr, numNodes)) )
307             throw DataException(msgPrefix+"put(Nodes_grDfI)");
308    
309          // Nodes grNI
310          if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
311             throw DataException(msgPrefix+"add_var(Nodes_grNI)");
312          int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
313          if (! (ids->put(int_ptr, numNodes)) )
314             throw DataException(msgPrefix+"put(Nodes_grNI)");
315    
316          // Nodes Coordinates
317          if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
318             throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
319          if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
320             throw DataException(msgPrefix+"put(Nodes_Coordinates)");
321    
322     }     }
323    
# Line 254  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    
328       // Temp storage to gather node IDs        // Elements_Id
329       int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
330             throw DataException(msgPrefix+"add_var(Elements_Id)");
331       // Elements ReferenceElement->ElementTypeId        int_ptr = &mesh->Elements->Id[0];
332       if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )        if (! (ids->put(int_ptr, num_Elements)) )
333          throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Id)");
334    
335       // Elements_Id        // Elements_Tag
336       if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
337          throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Tag)");
338       int_ptr = &mesh->Elements->Id[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_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Tag)");
341    
342       // Elements_Tag        // Elements_Owner
343       if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
344          throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Owner)");
345       int_ptr = &mesh->Elements->Tag[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_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Owner)");
348    
349       // Elements_Owner        // Elements_Color
350       if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
351          throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Color)");
352       int_ptr = &mesh->Elements->Owner[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_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Color)");
355    
356       // Elements_Color        // Elements_Nodes
357       if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
358          throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Nodes)");
359       int_ptr = &mesh->Elements->Color[0];        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
360       if (! (ids->put(int_ptr, num_Elements)) )           throw DataException(msgPrefix+"put(Elements_Nodes)");
         throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);  
   
      // 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)]];  
      if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )  
         throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);  
      if (! (ids->put(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes)) )  
         throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);  
   
      TMPMEMFREE(Elements_Nodes);  
361    
362     }     }
363    
# Line 306  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    
368       // Temp storage to gather node IDs        // FaceElements_Id
369       int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
370             throw DataException(msgPrefix+"add_var(FaceElements_Id)");
371       // FaceElements ReferenceElement->ElementTypeId        int_ptr = &mesh->FaceElements->Id[0];
372       if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )        if (! (ids->put(int_ptr, num_FaceElements)) )
373          throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Id)");
374    
375       // FaceElements_Id        // FaceElements_Tag
376       if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
377          throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
378       int_ptr = &mesh->FaceElements->Id[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_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Tag)");
381    
382       // FaceElements_Tag        // FaceElements_Owner
383       if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
384          throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
385       int_ptr = &mesh->FaceElements->Tag[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_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Owner)");
388    
389       // FaceElements_Owner        // FaceElements_Color
390       if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
391          throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Color)");
392       int_ptr = &mesh->FaceElements->Owner[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_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Color)");
395    
396       // FaceElements_Color        // FaceElements_Nodes
397       if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
398          throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
399       int_ptr = &mesh->FaceElements->Color[0];        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
400       if (! (ids->put(int_ptr, num_FaceElements)) )           throw DataException(msgPrefix+"put(FaceElements_Nodes)");
         throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);  
   
      // 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)]];  
      if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )  
         throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);  
      if (! (ids->put(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )  
         throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);  
   
      TMPMEMFREE(FaceElements_Nodes);  
401    
402     }     }
403    
# Line 358  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    
408       // Temp storage to gather node IDs        // ContactElements_Id
409       int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
410             throw DataException(msgPrefix+"add_var(ContactElements_Id)");
411       // ContactElements ReferenceElement->ElementTypeId        int_ptr = &mesh->ContactElements->Id[0];
412       if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )        if (! (ids->put(int_ptr, num_ContactElements)) )
413          throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Id)");
414    
415       // ContactElements_Id        // ContactElements_Tag
416       if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
417          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Tag)");
418       int_ptr = &mesh->ContactElements->Id[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_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Tag)");
421    
422       // ContactElements_Tag        // ContactElements_Owner
423       if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
424          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Owner)");
425       int_ptr = &mesh->ContactElements->Tag[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_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Owner)");
428    
429       // ContactElements_Owner        // ContactElements_Color
430       if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
431          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Color)");
432       int_ptr = &mesh->ContactElements->Owner[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_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Color)");
435    
436       // ContactElements_Color        // ContactElements_Nodes
437       if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
438          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Nodes)");
439       int_ptr = &mesh->ContactElements->Color[0];        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
440       if (! (ids->put(int_ptr, num_ContactElements)) )           throw DataException(msgPrefix+"put(ContactElements_Nodes)");
         throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);  
   
      // 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)]];  
      if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )  
         throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);  
      if (! (ids->put(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )  
         throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);  
   
      TMPMEMFREE(ContactElements_Nodes);  
441    
442     }     }
443    
# Line 410  void MeshAdapter::dump(const std::string Line 445  void MeshAdapter::dump(const std::string
445    
446     if (num_Points>0) {     if (num_Points>0) {
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");
   
      // Temp storage to gather node IDs  
      int *Points_Nodes = TMPMEMALLOC(num_Points,int);  
   
      // Points ReferenceElement->ElementTypeId  
      if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )  
         throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);  
   
      // Points_Id  
      if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )  
         throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);  
      int_ptr = &mesh->Points->Id[0];  
      if (! (ids->put(int_ptr, num_Points)) )  
         throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);  
   
      // Points_Tag  
      if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )  
         throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);  
      int_ptr = &mesh->Points->Tag[0];  
      if (! (ids->put(int_ptr, num_Points)) )  
         throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);  
   
      // Points_Owner  
      if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )  
         throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);  
      int_ptr = &mesh->Points->Owner[0];  
      if (! (ids->put(int_ptr, num_Points)) )  
         throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);  
   
      // Points_Color  
      if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )  
         throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);  
      int_ptr = &mesh->Points->Color[0];  
      if (! (ids->put(int_ptr, num_Points)) )  
         throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);  
   
      // Points_Nodes  
      // 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)]];  
      if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )  
         throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);  
      if (! (ids->put(&(Points_Nodes[0]), num_Points)) )  
         throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);  
449    
450       TMPMEMFREE(Points_Nodes);        // Points_Id
451          if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
452             throw DataException(msgPrefix+"add_var(Points_Id)");
453          int_ptr = &mesh->Points->Id[0];
454          if (! (ids->put(int_ptr, num_Points)) )
455             throw DataException(msgPrefix+"put(Points_Id)");
456    
457          // Points_Tag
458          if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
459             throw DataException(msgPrefix+"add_var(Points_Tag)");
460          int_ptr = &mesh->Points->Tag[0];
461          if (! (ids->put(int_ptr, num_Points)) )
462             throw DataException(msgPrefix+"put(Points_Tag)");
463    
464          // Points_Owner
465          if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
466             throw DataException(msgPrefix+"add_var(Points_Owner)");
467          int_ptr = &mesh->Points->Owner[0];
468          if (! (ids->put(int_ptr, num_Points)) )
469             throw DataException(msgPrefix+"put(Points_Owner)");
470    
471          // Points_Color
472          if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
473             throw DataException(msgPrefix+"add_var(Points_Color)");
474          int_ptr = &mesh->Points->Color[0];
475          if (! (ids->put(int_ptr, num_Points)) )
476             throw DataException(msgPrefix+"put(Points_Color)");
477    
478          // Points_Nodes
479          // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
480          if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
481             throw DataException(msgPrefix+"add_var(Points_Nodes)");
482          if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
483             throw DataException(msgPrefix+"put(Points_Nodes)");
484    
485     }     }
486    
487     // // // // // TagMap // // // // //     // // // // // TagMap // // // // //
488    
489       fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump does not save the TagMap yet\n\n\n");     if (num_Tags>0) {
490    
491          // Temp storage to gather node IDs
492          int *Tags_keys = TMPMEMALLOC(num_Tags, int);
493          char name_temp[4096];
494    
495          /* Copy tag data into temp arrays */
496          tag_map = mesh->TagMap;
497          if (tag_map) {
498             int i = 0;
499             while (tag_map) {
500                Tags_keys[i++] = tag_map->tag_key;
501                tag_map=tag_map->next;
502             }
503          }
504    
505          // Tags_keys
506          if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
507             throw DataException(msgPrefix+"add_var(Tags_keys)");
508          int_ptr = &Tags_keys[0];
509          if (! (ids->put(int_ptr, num_Tags)) )
510             throw DataException(msgPrefix+"put(Tags_keys)");
511    
512          // Tags_names_*
513          // This is an array of strings, it should be stored as an array but
514          // 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;
517          if (tag_map) {
518             int i = 0;
519             while (tag_map) {
520                sprintf(name_temp, "Tags_name_%d", i);
521                if (!dataFile.add_att(name_temp, tag_map->name) )
522                   throw DataException(msgPrefix+"add_att(Tags_names_XX)");
523                tag_map=tag_map->next;
524                i++;
525             }
526          }
527    
528          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 473  void MeshAdapter::dump(const std::string Line 543  void MeshAdapter::dump(const std::string
543    
544  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
545  {  {
546    return "FinleyMesh";     return "FinleyMesh";
547  }  }
548    
549  string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const  string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
550  {  {
551    FunctionSpaceNamesMapType::iterator loc;     FunctionSpaceNamesMapType::iterator loc;
552    loc=m_functionSpaceTypeNames.find(functionSpaceType);     loc=m_functionSpaceTypeNames.find(functionSpaceType);
553    if (loc==m_functionSpaceTypeNames.end()) {     if (loc==m_functionSpaceTypeNames.end()) {
554      return "Invalid function space type code.";        return "Invalid function space type code.";
555    } else {     } else {
556      return loc->second;        return loc->second;
557    }     }
558  }  }
559    
560  bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const  bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
561  {  {
562    FunctionSpaceNamesMapType::iterator loc;     FunctionSpaceNamesMapType::iterator loc;
563    loc=m_functionSpaceTypeNames.find(functionSpaceType);     loc=m_functionSpaceTypeNames.find(functionSpaceType);
564    return (loc!=m_functionSpaceTypeNames.end());     return (loc!=m_functionSpaceTypeNames.end());
565  }  }
566    
567  void MeshAdapter::setFunctionSpaceTypeNames()  void MeshAdapter::setFunctionSpaceTypeNames()
568  {  {
569    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
570      (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));     (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));
571    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
572      (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));     (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));
573    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
574      (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));     (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
575    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
576      (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));     (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));
577    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
578      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));     (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
579    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
580      (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));     (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
581    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
582      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));     (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
583    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
584      (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));     (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
585    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
586      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));     (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
587    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
588      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));     (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
589    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
590      (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));     (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
591    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
592      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));     (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
593    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
594      (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));     (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
595  }  }
596    
597  int MeshAdapter::getContinuousFunctionCode() const  int MeshAdapter::getContinuousFunctionCode() const
598  {  {
599    return Nodes;     return Nodes;
600  }  }
601  int MeshAdapter::getReducedContinuousFunctionCode() const  int MeshAdapter::getReducedContinuousFunctionCode() const
602  {  {
603    return ReducedNodes;     return ReducedNodes;
604  }  }
605    
606  int MeshAdapter::getFunctionCode() const  int MeshAdapter::getFunctionCode() const
607  {  {
608    return Elements;     return Elements;
609  }  }
610  int MeshAdapter::getReducedFunctionCode() const  int MeshAdapter::getReducedFunctionCode() const
611  {  {
612    return ReducedElements;     return ReducedElements;
613  }  }
614    
615  int MeshAdapter::getFunctionOnBoundaryCode() const  int MeshAdapter::getFunctionOnBoundaryCode() const
616  {  {
617    return FaceElements;     return FaceElements;
618  }  }
619  int MeshAdapter::getReducedFunctionOnBoundaryCode() const  int MeshAdapter::getReducedFunctionOnBoundaryCode() const
620  {  {
621    return ReducedFaceElements;     return ReducedFaceElements;
622  }  }
623    
624  int MeshAdapter::getFunctionOnContactZeroCode() const  int MeshAdapter::getFunctionOnContactZeroCode() const
625  {  {
626    return ContactElementsZero;     return ContactElementsZero;
627  }  }
628  int MeshAdapter::getReducedFunctionOnContactZeroCode() const  int MeshAdapter::getReducedFunctionOnContactZeroCode() const
629  {  {
630    return ReducedContactElementsZero;     return ReducedContactElementsZero;
631  }  }
632    
633  int MeshAdapter::getFunctionOnContactOneCode() const  int MeshAdapter::getFunctionOnContactOneCode() const
634  {  {
635    return ContactElementsOne;     return ContactElementsOne;
636  }  }
637  int MeshAdapter::getReducedFunctionOnContactOneCode() const  int MeshAdapter::getReducedFunctionOnContactOneCode() const
638  {  {
639    return ReducedContactElementsOne;     return ReducedContactElementsOne;
640  }  }
641    
642  int MeshAdapter::getSolutionCode() const  int MeshAdapter::getSolutionCode() const
643  {  {
644    return DegreesOfFreedom;     return DegreesOfFreedom;
645  }  }
646    
647  int MeshAdapter::getReducedSolutionCode() const  int MeshAdapter::getReducedSolutionCode() const
648  {  {
649    return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
650  }  }
651    
652  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
653  {  {
654    return Points;     return Points;
655  }  }
656    
657  //  //
# Line 589  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    checkFinleyError();     int numDim=Finley_Mesh_getDim(mesh);
664    return numDim;     checkFinleyError();
665       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  //  //
# Line 604  pair<int,int> MeshAdapter::getDataShape( Line 683  pair<int,int> MeshAdapter::getDataShape(
683     int numSamples=0;     int numSamples=0;
684     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
685     switch (functionSpaceCode) {     switch (functionSpaceCode) {
686        case(Nodes):     case(Nodes):
687             numDataPointsPerSample=1;     numDataPointsPerSample=1;
688             numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);     numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);
689             break;     break;
690        case(ReducedNodes):     case(ReducedNodes):
691             numDataPointsPerSample=1;     numDataPointsPerSample=1;
692             numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);     numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
693             break;     break;
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;
718        case(Points):     case(Points):
719             if (mesh->Points!=NULL) {     if (mesh->Points!=NULL) {
720               numDataPointsPerSample=1;        numDataPointsPerSample=1;
721               numSamples=mesh->Points->numElements;        numSamples=mesh->Points->numElements;
722             }     }
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;
748        case(DegreesOfFreedom):     case(DegreesOfFreedom):
749             if (mesh->Nodes!=NULL) {     if (mesh->Nodes!=NULL) {
750               numDataPointsPerSample=1;        numDataPointsPerSample=1;
751               numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);        numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
752             }     }
753             break;     break;
754        case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
755             if (mesh->Nodes!=NULL) {     if (mesh->Nodes!=NULL) {
756               numDataPointsPerSample=1;        numDataPointsPerSample=1;
757               numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);        numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
758             }     }
759             break;     break;
760        default:     default:
761          stringstream temp;        stringstream temp;
762          temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
763          throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
764          break;        break;
765        }     }
766        return pair<int,int>(numDataPointsPerSample,numSamples);     return pair<int,int>(numDataPointsPerSample,numSamples);
767  }  }
768    
769  //  //
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 707  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 750  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 759  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
868    //
869    void MeshAdapter::addPDEToTransportProblem(
870                                               AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
871                                               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,
873                                               const escript::Data& d, const escript::Data& y,
874                                               const escript::Data& d_contact,const escript::Data& y_contact,
875                                               const escript::Data& d_dirac, const escript::Data& y_dirac) const
876    {
877       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
878       if (tpa==0)
879       {
880        throw FinleyAdapterException("finley only supports Paso transport problems.");
881       }
882    
883    
884       DataTypes::ShapeType shape;
885       source.expand();
886       escriptDataC _source=source.getDataC();
887       escriptDataC _M=M.getDataC();
888       escriptDataC _A=A.getDataC();
889       escriptDataC _B=B.getDataC();
890       escriptDataC _C=C.getDataC();
891       escriptDataC _D=D.getDataC();
892       escriptDataC _X=X.getDataC();
893       escriptDataC _Y=Y.getDataC();
894       escriptDataC _d=d.getDataC();
895       escriptDataC _y=y.getDataC();
896       escriptDataC _d_contact=d_contact.getDataC();
897       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();
903       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 );
906       checkFinleyError();
907    
908       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
909       checkFinleyError();
910    
911       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
912       checkFinleyError();
913    
914       Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
915       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 766  void MeshAdapter::addPDEToRHS( escript:: Line 924  void MeshAdapter::addPDEToRHS( escript::
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)
932      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");        throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
933    
934    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
935    escriptDataC _target=target.getDataC();     escriptDataC _target=target.getDataC();
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);
945                 break;        break;
946             case(Elements):        case(Elements):
947             case(ReducedElements):        case(ReducedElements):
948                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
949                 break;        break;
950             case(FaceElements):        case(FaceElements):
951             case(ReducedFaceElements):        case(ReducedFaceElements):
952                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
953                 break;        break;
954             case(Points):        case(Points):
955                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
956                 break;        break;
957             case(ContactElementsZero):        case(ContactElementsZero):
958             case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
959                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
960                 break;        break;
961             case(ContactElementsOne):        case(ContactElementsOne):
962             case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
963                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
964                 break;        break;
965             default:        default:
966                 stringstream temp;           stringstream temp;
967                 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
968                 throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
969                 break;           break;
970          }        }
971          break;        break;
972       case(ReducedNodes):     case(ReducedNodes):
973          switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
974             case(Nodes):        case(Nodes):
975             case(ReducedNodes):        case(ReducedNodes):
976             case(DegreesOfFreedom):        case(DegreesOfFreedom):
977             case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
978                 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
979                 break;        break;
980             case(Elements):        case(Elements):
981             case(ReducedElements):        case(ReducedElements):
982                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
983                 break;        break;
984             case(FaceElements):        case(FaceElements):
985             case(ReducedFaceElements):        case(ReducedFaceElements):
986                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
987                 break;        break;
988             case(Points):        case(Points):
989                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
990                 break;        break;
991             case(ContactElementsZero):        case(ContactElementsZero):
992             case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
993                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
994                 break;        break;
995             case(ContactElementsOne):        case(ContactElementsOne):
996             case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
997                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
998                 break;        break;
999             default:        default:
1000                 stringstream temp;           stringstream temp;
1001                 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1002                 throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
1003                 break;           break;
1004          }        }
1005          break;        break;
1006       case(Elements):     case(Elements):
1007          if (target.getFunctionSpace().getTypeCode()==Elements) {        if (target.getFunctionSpace().getTypeCode()==Elements) {
1008             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1009          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1010             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
1011          } else {        } else {
1012             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
1013          }        }
1014          break;        break;
1015       case(ReducedElements):     case(ReducedElements):
1016          if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1017             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1018          } else {        } else {
1019             throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
1020          }        }
1021          break;        break;
1022       case(FaceElements):     case(FaceElements):
1023          if (target.getFunctionSpace().getTypeCode()==FaceElements) {        if (target.getFunctionSpace().getTypeCode()==FaceElements) {
1024             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1025          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1026             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1027          } else {        } else {
1028             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1029         }        }
1030         break;        break;
1031       case(ReducedFaceElements):     case(ReducedFaceElements):
1032          if (target.getFunctionSpace().getTypeCode()==FaceElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1033             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1034          } else {        } else {
1035             throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
1036         }        }
1037         break;        break;
1038       case(Points):     case(Points):
1039          if (target.getFunctionSpace().getTypeCode()==Points) {        if (target.getFunctionSpace().getTypeCode()==Points) {
1040             Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1041          } else {        } else {
1042             throw FinleyAdapterException("Error - No interpolation with data on points possible.");           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1043          }        }
1044          break;        break;
1045       case(ContactElementsZero):     case(ContactElementsZero):
1046       case(ContactElementsOne):     case(ContactElementsOne):
1047          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {        if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1048             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1049          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1050             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1051          } else {        } else {
1052             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1053             break;        }
1054          }        break;
1055          break;     case(ReducedContactElementsZero):
1056       case(ReducedContactElementsZero):     case(ReducedContactElementsOne):
1057       case(ReducedContactElementsOne):        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1058          if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1059             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);        } else {
1060          } else {           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1061             throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");        }
1062             break;        break;
1063          }     case(DegreesOfFreedom):      
1064          break;        switch(target.getFunctionSpace().getTypeCode()) {
1065       case(DegreesOfFreedom):              case(ReducedDegreesOfFreedom):
1066          switch(target.getFunctionSpace().getTypeCode()) {        case(DegreesOfFreedom):
1067             case(ReducedDegreesOfFreedom):        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1068             case(DegreesOfFreedom):        break;
1069                Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);    
1070                break;        case(Nodes):
1071          case(ReducedNodes):
1072             case(Nodes):        if (getMPISize()>1) {
1073             case(ReducedNodes):           escript::Data temp=escript::Data(in);
1074                if (getMPISize()>1) {           temp.expand();
1075                    escript::Data temp=escript::Data(in);           escriptDataC _in2 = temp.getDataC();
1076                    temp.expand();           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1077                    escriptDataC _in2 = temp.getDataC();        } else {
1078                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1079                } else {        }
1080                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        break;
1081                }        case(Elements):
1082                break;        case(ReducedElements):
1083             case(Elements):        if (getMPISize()>1) {
1084             case(ReducedElements):           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1085                if (getMPISize()>1) {           escriptDataC _in2 = temp.getDataC();
1086                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1087                   escriptDataC _in2 = temp.getDataC();        } else {
1088                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1089                } else {        }
1090                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        break;
1091                }        case(FaceElements):
1092                break;        case(ReducedFaceElements):
1093             case(FaceElements):        if (getMPISize()>1) {
1094             case(ReducedFaceElements):           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1095                if (getMPISize()>1) {           escriptDataC _in2 = temp.getDataC();
1096                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1097                   escriptDataC _in2 = temp.getDataC();    
1098                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);        } else {
1099             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1100                } else {        }
1101                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
1102                }        case(Points):
1103                break;        if (getMPISize()>1) {
1104             case(Points):           //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1105                if (getMPISize()>1) {           //escriptDataC _in2 = temp.getDataC();
1106                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );        } else {
1107                   escriptDataC _in2 = temp.getDataC();           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1108                } else {        }
1109                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        break;
1110                }        case(ContactElementsZero):
1111                break;        case(ContactElementsOne):
1112             case(ContactElementsZero):        case(ReducedContactElementsZero):
1113             case(ContactElementsOne):        case(ReducedContactElementsOne):
1114             case(ReducedContactElementsZero):        if (getMPISize()>1) {
1115             case(ReducedContactElementsOne):           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1116                if (getMPISize()>1) {           escriptDataC _in2 = temp.getDataC();
1117                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1118                   escriptDataC _in2 = temp.getDataC();        } else {
1119                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1120                } else {        }
1121                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        break;
1122                }        default:
1123               break;           stringstream temp;
1124             default:           temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1125               stringstream temp;           throw FinleyAdapterException(temp.str());
1126               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           break;
1127               throw FinleyAdapterException(temp.str());        }
1128               break;        break;
1129          }     case(ReducedDegreesOfFreedom):
1130          break;        switch(target.getFunctionSpace().getTypeCode()) {
1131       case(ReducedDegreesOfFreedom):        case(Nodes):
1132         switch(target.getFunctionSpace().getTypeCode()) {        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1133            case(Nodes):        break;
1134               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");        case(ReducedNodes):
1135               break;        if (getMPISize()>1) {
1136            case(ReducedNodes):           escript::Data temp=escript::Data(in);
1137                if (getMPISize()>1) {           temp.expand();
1138                    escript::Data temp=escript::Data(in);           escriptDataC _in2 = temp.getDataC();
1139                    temp.expand();           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1140                    escriptDataC _in2 = temp.getDataC();        } else {
1141                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1142                } else {        }
1143                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        break;
1144                }        case(DegreesOfFreedom):
1145                break;        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1146            case(DegreesOfFreedom):        break;
1147               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");        case(ReducedDegreesOfFreedom):
1148               break;        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1149            case(ReducedDegreesOfFreedom):        break;
1150               Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(Elements):
1151               break;        case(ReducedElements):
1152            case(Elements):        if (getMPISize()>1) {
1153            case(ReducedElements):           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1154                if (getMPISize()>1) {           escriptDataC _in2 = temp.getDataC();
1155                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1156                   escriptDataC _in2 = temp.getDataC();        } else {
1157                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1158                } else {        }
1159                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        break;
1160               }        case(FaceElements):
1161               break;        case(ReducedFaceElements):
1162            case(FaceElements):        if (getMPISize()>1) {
1163            case(ReducedFaceElements):           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1164                if (getMPISize()>1) {           escriptDataC _in2 = temp.getDataC();
1165                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1166                   escriptDataC _in2 = temp.getDataC();        } else {
1167                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1168                } else {        }
1169                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
1170                }        case(Points):
1171               break;        if (getMPISize()>1) {
1172            case(Points):           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1173                if (getMPISize()>1) {           escriptDataC _in2 = temp.getDataC();
1174                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1175                   escriptDataC _in2 = temp.getDataC();        } else {
1176                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1177                } else {        }
1178                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        break;
1179               }        case(ContactElementsZero):
1180               break;        case(ContactElementsOne):
1181            case(ContactElementsZero):        case(ReducedContactElementsZero):
1182            case(ContactElementsOne):        case(ReducedContactElementsOne):
1183            case(ReducedContactElementsZero):        if (getMPISize()>1) {
1184            case(ReducedContactElementsOne):           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1185                if (getMPISize()>1) {           escriptDataC _in2 = temp.getDataC();
1186                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1187                   escriptDataC _in2 = temp.getDataC();        } else {
1188                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1189                } else {        }
1190                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        break;
1191               }        default:
1192               break;           stringstream temp;
1193            default:           temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1194               stringstream temp;           throw FinleyAdapterException(temp.str());
1195               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           break;
1196               throw FinleyAdapterException(temp.str());        }
1197               break;        break;
1198         }     default:
1199         break;        stringstream temp;
1200       default:        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1201          stringstream temp;        throw FinleyAdapterException(temp.str());
1202          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();        break;
1203          throw FinleyAdapterException(temp.str());     }
1204          break;     checkFinleyError();
   }  
   checkFinleyError();  
1205  }  }
1206    
1207  //  //
# Line 1053  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();
1216    // in case of values node coordinates we can do the job directly:     // in case of values node coordinates we can do the job directly:
1217    if (arg.getFunctionSpace().getTypeCode()==Nodes) {     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
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:
1225       interpolateOnDomain(arg,tmp_data);        interpolateOnDomain(arg,tmp_data);
1226    }     }
1227    checkFinleyError();     checkFinleyError();
1228  }  }
1229    
1230  //  //
# Line 1076  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    if (normalDomain!=*this)     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1237       throw FinleyAdapterException("Error - Illegal domain of normal locations");     if (normalDomain!=*this)
1238    Finley_Mesh* mesh=m_finleyMesh.get();        throw FinleyAdapterException("Error - Illegal domain of normal locations");
1239    escriptDataC _normal=normal.getDataC();     Finley_Mesh* mesh=m_finleyMesh.get();
1240    switch(normal.getFunctionSpace().getTypeCode()) {     escriptDataC _normal=normal.getDataC();
1241      case(Nodes):     switch(normal.getFunctionSpace().getTypeCode()) {
1242        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");     case(Nodes):
1243        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1244      case(ReducedNodes):     break;
1245        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");     case(ReducedNodes):
1246        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1247      case(Elements):     break;
1248        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");     case(Elements):
1249        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1250      case(ReducedElements):     break;
1251        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");     case(ReducedElements):
1252        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1253      case (FaceElements):     break;
1254        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     case (FaceElements):
1255        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1256      case (ReducedFaceElements):     break;
1257        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     case (ReducedFaceElements):
1258        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1259      case(Points):     break;
1260        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");     case(Points):
1261        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1262      case (ContactElementsOne):     break;
1263      case (ContactElementsZero):     case (ContactElementsOne):
1264        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);     case (ContactElementsZero):
1265        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1266      case (ReducedContactElementsOne):     break;
1267      case (ReducedContactElementsZero):     case (ReducedContactElementsOne):
1268        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);     case (ReducedContactElementsZero):
1269        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1270      case(DegreesOfFreedom):     break;
1271        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");     case(DegreesOfFreedom):
1272        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1273      case(ReducedDegreesOfFreedom):     break;
1274        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");     case(ReducedDegreesOfFreedom):
1275        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1276      default:     break;
1277       default:
1278        stringstream temp;        stringstream temp;
1279        temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();        temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1280        throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1281        break;        break;
1282    }     }
1283    checkFinleyError();     checkFinleyError();
1284  }  }
1285    
1286  //  //
# Line 1131  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       throw FinleyAdapterException("Error - Illegal domain of interpolation target");     if (targetDomain!=this)
1294          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.");
1297  }  }
1298    
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    
1308    double blocktimer_start = blocktimer_time();     double blocktimer_start = blocktimer_time();
1309    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1310    escriptDataC _temp;     escriptDataC _temp;
1311    escript::Data temp;     escript::Data temp;
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;
1324       case(Elements):     case(Elements):
1325          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1326          break;     break;
1327       case(ReducedElements):     case(ReducedElements):
1328          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1329          break;     break;
1330       case(FaceElements):     case(FaceElements):
1331          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1332          break;     break;
1333       case(ReducedFaceElements):     case(ReducedFaceElements):
1334          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1335          break;     break;
1336       case(Points):     case(Points):
1337          throw FinleyAdapterException("Error - Integral of data on points is not supported.");     throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1338          break;     break;
1339       case(ContactElementsZero):     case(ContactElementsZero):
1340          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1341          break;     break;
1342       case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1343          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1344          break;     break;
1345       case(ContactElementsOne):     case(ContactElementsOne):
1346          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1347          break;     break;
1348       case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
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;
1361       default:     default:
1362          stringstream temp;        stringstream temp;
1363          temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();        temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1364          throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1365          break;        break;
1366    }     }
1367    checkFinleyError();     checkFinleyError();
1368    blocktimer_increment("integrate()", blocktimer_start);     blocktimer_increment("integrate()", blocktimer_start);
1369  }  }
1370    
1371  //  //
# Line 1215  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    
1383    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1384    escriptDataC _grad=grad.getDataC();     escriptDataC _grad=grad.getDataC();
1385    escriptDataC nodeDataC;     escriptDataC nodeDataC;
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();
1396        }        }
1397    } else {     } else {
1398       nodeDataC = arg.getDataC();        nodeDataC = arg.getDataC();
1399    }     }
1400    switch(grad.getFunctionSpace().getTypeCode()) {     switch(grad.getFunctionSpace().getTypeCode()) {
1401         case(Nodes):     case(Nodes):
1402            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");     throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1403            break;     break;
1404         case(ReducedNodes):     case(ReducedNodes):
1405            throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1406            break;     break;
1407         case(Elements):     case(Elements):
1408            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1409            break;     break;
1410         case(ReducedElements):     case(ReducedElements):
1411            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1412            break;     break;
1413         case(FaceElements):     case(FaceElements):
1414            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1415            break;     break;
1416         case(ReducedFaceElements):     case(ReducedFaceElements):
1417            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1418            break;     break;
1419         case(Points):     case(Points):
1420            throw FinleyAdapterException("Error - Gradient at points is not supported.");     throw FinleyAdapterException("Error - Gradient at points is not supported.");
1421            break;     break;
1422         case(ContactElementsZero):     case(ContactElementsZero):
1423            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1424            break;     break;
1425         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1426            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1427            break;     break;
1428         case(ContactElementsOne):     case(ContactElementsOne):
1429            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1430            break;     break;
1431         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1432            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1433            break;     break;
1434         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1435            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1436            break;     break;
1437         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1438            throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1439            break;     break;
1440         default:     default:
1441            stringstream temp;        stringstream temp;
1442            temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();        temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1443            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1444            break;        break;
1445    }     }
1446    checkFinleyError();     checkFinleyError();
1447  }  }
1448    
1449  //  //
# Line 1293  void MeshAdapter::setToGradient(escript: Line 1451  void MeshAdapter::setToGradient(escript:
1451  //  //
1452  void MeshAdapter::setToSize(escript::Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
1453  {  {
1454    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1455    escriptDataC tmp=size.getDataC();     escriptDataC tmp=size.getDataC();
1456    switch(size.getFunctionSpace().getTypeCode()) {     switch(size.getFunctionSpace().getTypeCode()) {
1457         case(Nodes):     case(Nodes):
1458            throw FinleyAdapterException("Error - Size of nodes is not supported.");     throw FinleyAdapterException("Error - Size of nodes is not supported.");
1459            break;     break;
1460         case(ReducedNodes):     case(ReducedNodes):
1461            throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1462            break;     break;
1463         case(Elements):     case(Elements):
1464            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1465            break;     break;
1466         case(ReducedElements):     case(ReducedElements):
1467            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1468            break;     break;
1469         case(FaceElements):     case(FaceElements):
1470            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1471            break;     break;
1472         case(ReducedFaceElements):     case(ReducedFaceElements):
1473            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1474            break;     break;
1475         case(Points):     case(Points):
1476            throw FinleyAdapterException("Error - Size of point elements is not supported.");     throw FinleyAdapterException("Error - Size of point elements is not supported.");
1477            break;     break;
1478         case(ContactElementsZero):     case(ContactElementsZero):
1479         case(ContactElementsOne):     case(ContactElementsOne):
1480            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1481            break;     break;
1482         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1483         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1484            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1485            break;     break;
1486         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1487            throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1488            break;     break;
1489         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1490            throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1491            break;     break;
1492         default:     default:
1493            stringstream temp;        stringstream temp;
1494            temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();        temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1495            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1496            break;        break;
1497    }     }
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    checkFinleyError();         Finley_Mesh_setCoordinates(mesh,&tmp);
1514  }     } else {
1515           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1516  // saves a data array in openDX format:         tmp = new_x_inter.getDataC();
1517  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const         Finley_Mesh_setCoordinates(mesh,&tmp);
1518  {     }
1519      int MAX_namelength=256;     checkFinleyError();
1520      const int num_data=boost::python::extract<int>(arg.attr("__len__")());  }
   /* win32 refactor */  
   char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;  
   for(int i=0;i<num_data;i++)  
   {  
     names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
   }  
   
   char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;  
   escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  
   escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  
   
   boost::python::list keys=arg.keys();  
   for (int i=0;i<num_data;++i) {  
          std::string n=boost::python::extract<std::string>(keys[i]);  
          escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);  
          if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)  
              throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");  
          data[i]=d.getDataC();  
          ptr_data[i]=&(data[i]);  
          c_names[i]=&(names[i][0]);  
          if (n.length()>MAX_namelength-1) {  
             strncpy(c_names[i],n.c_str(),MAX_namelength-1);  
             c_names[i][MAX_namelength-1]='\0';  
          } else {  
             strcpy(c_names[i],n.c_str());  
          }  
     }  
     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);  
     checkFinleyError();  
       
       /* win32 refactor */  
   TMPMEMFREE(c_names);  
   TMPMEMFREE(data);  
   TMPMEMFREE(ptr_data);  
   for(int i=0;i<num_data;i++)  
   {  
     TMPMEMFREE(names[i]);  
   }  
   TMPMEMFREE(names);  
1521    
1522      return;  //
1523    // 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       numData = boost::python::extract<int>(arg.attr("__len__")());
1529       /* win32 refactor */
1530       names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1531       data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1532       dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1533    
1534       boost::python::list keys=arg.keys();
1535       for (int i=0; i<numData; ++i) {
1536          string n=boost::python::extract<string>(keys[i]);
1537          escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1538          if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1539             throw FinleyAdapterException("Error: Data must be defined on same Domain");
1540          data[i] = d.getDataC();
1541          dataPtr[i] = &(data[i]);
1542          names[i] = TMPMEMALLOC(n.length()+1, char);
1543          strcpy(names[i], n.c_str());
1544       }
1545  }  }
1546    
1547  // saves a data array in openVTK format:  //
1548  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  // 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 MAX_namelength=256;     int num_data;
1553      const int num_data=boost::python::extract<int>(arg.attr("__len__")());     char **names;
1554    /* win32 refactor */     escriptDataC *data;
1555    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     escriptDataC **ptr_data;
   for(int i=0;i<num_data;i++)  
   {  
     names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
   }  
1556    
1557    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1558    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;     Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1559    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;     checkFinleyError();
1560    
1561      boost::python::list keys=arg.keys();     /* win32 refactor */
1562      for (int i=0;i<num_data;++i) {     TMPMEMFREE(data);
1563           std::string n=boost::python::extract<std::string>(keys[i]);     TMPMEMFREE(ptr_data);
1564           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);     for(int i=0; i<num_data; i++)
1565           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)     {
1566               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");        TMPMEMFREE(names[i]);
1567           data[i]=d.getDataC();     }
1568           ptr_data[i]=&(data[i]);     TMPMEMFREE(names);
          c_names[i]=&(names[i][0]);  
          if (n.length()>MAX_namelength-1) {  
             strncpy(c_names[i],n.c_str(),MAX_namelength-1);  
             c_names[i][MAX_namelength-1]='\0';  
          } else {  
             strcpy(c_names[i],n.c_str());  
          }  
     }  
     Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);  
1569    
1570    checkFinleyError();     return;
1571    /* win32 refactor */  }
   TMPMEMFREE(c_names);  
   TMPMEMFREE(data);  
   TMPMEMFREE(ptr_data);  
   for(int i=0;i<num_data;i++)  
   {  
     TMPMEMFREE(names[i]);  
   }  
   TMPMEMFREE(names);  
1572    
1573      return;  //
1574    // 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        boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1579        pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1580                  metadata, metadata_schema, arg);
1581  }  }
1582                                                                                                                                                                            
1583                                                                                                                                                                            bool MeshAdapter::ownSample(int fs_code, index_t id) const
1584  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  {
1585  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  #ifdef ESYS_MPI
1586                        const int row_blocksize,      index_t myFirstNode=0, myLastNode=0, k=0;
1587                        const escript::FunctionSpace& row_functionspace,      index_t* globalNodeIndex=0;
1588                        const int column_blocksize,      Finley_Mesh* mesh_p=m_finleyMesh.get();
1589                        const escript::FunctionSpace& column_functionspace,      if (fs_code == FINLEY_REDUCED_NODES)
1590                        const int type) const      {
1591  {      myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1592      int reduceRowOrder=0;      myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1593      int reduceColOrder=0;      globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
     // is the domain right?  
     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());  
     if (row_domain!=*this)  
           throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");  
     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());  
     if (col_domain!=*this)  
           throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");  
     // is the function space type right  
     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {  
         reduceRowOrder=0;  
     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {  
         reduceRowOrder=1;  
     } else {  
         throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");  
1594      }      }
1595      if (column_functionspace.getTypeCode()==DegreesOfFreedom) {      else
1596          reduceColOrder=0;      {
1597      } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1598          reduceColOrder=1;      myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1599      } else {      globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
         throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");  
1600      }      }
1601      // generate matrix:      k=globalNodeIndex[id];
1602            return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1603      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);  #endif
1604      checkFinleyError();      return true;
1605      Paso_SystemMatrix* fsystemMatrix;  }
1606      int trilinos = 0;  
1607      if (trilinos) {  
1608    
1609    //
1610    // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1611    //
1612    ASM_ptr MeshAdapter::newSystemMatrix(
1613                                                     const int row_blocksize,
1614                                                     const escript::FunctionSpace& row_functionspace,
1615                                                     const int column_blocksize,
1616                                                     const escript::FunctionSpace& column_functionspace,
1617                                                     const int type) const
1618    {
1619       int reduceRowOrder=0;
1620       int reduceColOrder=0;
1621       // is the domain right?
1622       const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1623       if (row_domain!=*this)
1624          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()));
1626       if (col_domain!=*this)
1627          throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1628       // is the function space type right
1629       if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1630          reduceRowOrder=0;
1631       } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1632          reduceRowOrder=1;
1633       } else {
1634          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1635       }
1636       if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1637          reduceColOrder=0;
1638       } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1639          reduceColOrder=1;
1640       } else {
1641          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1642       }
1643       // generate matrix:
1644    
1645       Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1646       checkFinleyError();
1647       Paso_SystemMatrix* fsystemMatrix;
1648       int trilinos = 0;
1649       if (trilinos) {
1650  #ifdef TRILINOS  #ifdef TRILINOS
1651        /* Allocation Epetra_VrbMatrix here */        /* Allocation Epetra_VrbMatrix here */
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
1666    //
1667    ATP_ptr MeshAdapter::newTransportProblem(
1668                                                             const bool useBackwardEuler,
1669                                                             const int blocksize,
1670                                                             const escript::FunctionSpace& functionspace,
1671                                                             const int type) const
1672    {
1673       int reduceOrder=0;
1674       // is the domain right?
1675       const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1676       if (domain!=*this)
1677          throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1678       // is the function space type right
1679       if (functionspace.getTypeCode()==DegreesOfFreedom) {
1680          reduceOrder=0;
1681       } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1682          reduceOrder=1;
1683       } else {
1684          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1685       }
1686       // generate matrix:
1687    
1688       Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1689       checkFinleyError();
1690       Paso_TransportProblem* transportProblem;
1691       transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1692       checkPasoError();
1693       Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1694       TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1695       return ATP_ptr(tpa);
1696    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1697  }  }
1698    
1699  //  //
# Line 1510  SystemMatrixAdapter MeshAdapter::newSyst Line 1705  SystemMatrixAdapter MeshAdapter::newSyst
1705  // returns true if data at the atom_type is considered as being cell centered:  // returns true if data at the atom_type is considered as being cell centered:
1706  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1707  {  {
1708    switch(functionSpaceCode) {     switch(functionSpaceCode) {
1709         case(Nodes):     case(Nodes):
1710         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1711         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1712            return false;     return false;
1713            break;     break;
1714         case(Elements):     case(Elements):
1715         case(FaceElements):     case(FaceElements):
1716         case(Points):     case(Points):
1717         case(ContactElementsZero):     case(ContactElementsZero):
1718         case(ContactElementsOne):     case(ContactElementsOne):
1719         case(ReducedElements):     case(ReducedElements):
1720         case(ReducedFaceElements):     case(ReducedFaceElements):
1721         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1722         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1723            return true;     return true;
1724            break;     break;
1725         default:     default:
1726            stringstream temp;        stringstream temp;
1727            temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;        temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1728            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1729            break;        break;
1730    }     }
1731    checkFinleyError();     checkFinleyError();
1732    return false;     return false;
1733    }
1734    
1735    bool
1736    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1737    {
1738       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1739        class 1: DOF <-> Nodes
1740        class 2: ReducedDOF <-> ReducedNodes
1741        class 3: Points
1742        class 4: Elements
1743        class 5: ReducedElements
1744        class 6: FaceElements
1745        class 7: ReducedFaceElements
1746        class 8: ContactElementZero <-> ContactElementOne
1747        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1748    
1749       There is also a set of lines. Interpolation is possible down a line but not between lines.
1750       class 1 and 2 belong to all lines so aren't considered.
1751        line 0: class 3
1752        line 1: class 4,5
1753        line 2: class 6,7
1754        line 3: class 8,9
1755    
1756       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1757       eg hasnodes is true if we have at least one instance of Nodes.
1758       */
1759        if (fs.empty())
1760        {
1761            return false;
1762        }
1763        vector<int> hasclass(10);
1764        vector<int> hasline(4);
1765        bool hasnodes=false;
1766        bool hasrednodes=false;
1767        bool hascez=false;
1768        bool hasrcez=false;
1769        for (int i=0;i<fs.size();++i)
1770        {
1771        switch(fs[i])
1772        {
1773        case(Nodes):   hasnodes=true;   // no break is deliberate
1774        case(DegreesOfFreedom):
1775            hasclass[1]=1;
1776            break;
1777        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1778        case(ReducedDegreesOfFreedom):
1779            hasclass[2]=1;
1780            break;
1781        case(Points):
1782            hasline[0]=1;
1783            hasclass[3]=1;
1784            break;
1785        case(Elements):
1786            hasclass[4]=1;
1787            hasline[1]=1;
1788            break;
1789        case(ReducedElements):
1790            hasclass[5]=1;
1791            hasline[1]=1;
1792            break;
1793        case(FaceElements):
1794            hasclass[6]=1;
1795            hasline[2]=1;
1796            break;
1797        case(ReducedFaceElements):
1798            hasclass[7]=1;
1799            hasline[2]=1;
1800            break;
1801        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1802        case(ContactElementsOne):
1803            hasclass[8]=1;
1804            hasline[3]=1;
1805            break;
1806        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1807        case(ReducedContactElementsOne):
1808            hasclass[9]=1;
1809            hasline[3]=1;
1810            break;
1811        default:
1812            return false;
1813        }
1814        }
1815        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1816        // fail if we have more than one leaf group
1817    
1818        if (totlines>1)
1819        {
1820        return false;   // there are at least two branches we can't interpolate between
1821        }
1822        else if (totlines==1)
1823        {
1824        if (hasline[0]==1)      // we have points
1825        {
1826            resultcode=Points;
1827        }
1828        else if (hasline[1]==1)
1829        {
1830            if (hasclass[5]==1)
1831            {
1832            resultcode=ReducedElements;
1833            }
1834            else
1835            {
1836            resultcode=Elements;
1837            }
1838        }
1839        else if (hasline[2]==1)
1840        {
1841            if (hasclass[7]==1)
1842            {
1843            resultcode=ReducedFaceElements;
1844            }
1845            else
1846            {
1847            resultcode=FaceElements;
1848            }
1849        }
1850        else    // so we must be in line3
1851        {
1852            if (hasclass[9]==1)
1853            {
1854            // need something from class 9
1855            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1856            }
1857            else
1858            {
1859            // something from class 8
1860            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1861            }
1862        }
1863        }
1864        else    // totlines==0
1865        {
1866        if (hasclass[2]==1)
1867        {
1868            // something from class 2
1869            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1870        }
1871        else
1872        {   // something from class 1
1873            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1874        }
1875        }
1876        return true;
1877  }  }
1878    
1879  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1880  {  {
1881    switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1882       case(Nodes):     case(Nodes):
1883          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1884             case(Nodes):      case(Nodes):
1885             case(ReducedNodes):      case(ReducedNodes):
1886             case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1887             case(DegreesOfFreedom):      case(DegreesOfFreedom):
1888             case(Elements):      case(Elements):
1889             case(ReducedElements):      case(ReducedElements):
1890             case(FaceElements):      case(FaceElements):
1891             case(ReducedFaceElements):      case(ReducedFaceElements):
1892             case(Points):      case(Points):
1893             case(ContactElementsZero):      case(ContactElementsZero):
1894             case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1895             case(ContactElementsOne):      case(ContactElementsOne):
1896             case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1897                 return true;      return true;
1898             default:      default:
1899                 stringstream temp;            stringstream temp;
1900                 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1901                 throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1902          }     }
1903          break;     break;
1904       case(ReducedNodes):     case(ReducedNodes):
1905          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1906             case(ReducedNodes):      case(ReducedNodes):
1907             case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1908             case(Elements):      case(Elements):
1909             case(ReducedElements):      case(ReducedElements):
1910             case(FaceElements):      case(FaceElements):
1911             case(ReducedFaceElements):      case(ReducedFaceElements):
1912             case(Points):      case(Points):
1913             case(ContactElementsZero):      case(ContactElementsZero):
1914             case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1915             case(ContactElementsOne):      case(ContactElementsOne):
1916             case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1917                 return true;      return true;
1918            case(Nodes):      case(Nodes):
1919            case(DegreesOfFreedom):      case(DegreesOfFreedom):
1920               return false;      return false;
1921             default:      default:
1922                 stringstream temp;          stringstream temp;
1923                 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1924                 throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1925          }     }
1926          break;     break;
1927       case(Elements):     case(Elements):
1928          if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1929             return true;        return true;
1930          } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1931             return true;        return true;
         } else {  
            return false;  
         }  
      case(ReducedElements):  
         if (functionSpaceType_target==ReducedElements) {  
            return true;  
         } else {  
            return false;  
         }  
      case(FaceElements):  
         if (functionSpaceType_target==FaceElements) {  
            return true;  
         } else if (functionSpaceType_target==ReducedFaceElements) {  
            return true;  
         } else {  
            return false;  
         }  
      case(ReducedFaceElements):  
         if (functionSpaceType_target==ReducedFaceElements) {  
            return true;  
         } else {  
            return false;  
         }  
      case(Points):  
         if (functionSpaceType_target==Points) {  
            return true;  
         } else {  
            return false;  
         }  
      case(ContactElementsZero):  
      case(ContactElementsOne):  
         if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {  
            return true;  
         } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {  
            return true;  
         } else {  
            return false;  
         }  
      case(ReducedContactElementsZero):  
      case(ReducedContactElementsOne):  
         if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {  
            return true;  
1932          } else {          } else {
1933             return false;            return false;
         }  
      case(DegreesOfFreedom):  
         switch(functionSpaceType_target) {  
            case(ReducedDegreesOfFreedom):  
            case(DegreesOfFreedom):  
            case(Nodes):  
            case(ReducedNodes):  
            case(Elements):  
            case(ReducedElements):  
            case(Points):  
            case(FaceElements):  
            case(ReducedFaceElements):  
            case(ContactElementsZero):  
            case(ReducedContactElementsZero):  
            case(ContactElementsOne):  
            case(ReducedContactElementsOne):  
               return true;  
            default:  
              stringstream temp;  
              temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;  
              throw FinleyAdapterException(temp.str());  
1934          }          }
1935          break;     case(ReducedElements):
1936       case(ReducedDegreesOfFreedom):      if (functionSpaceType_target==ReducedElements) {
1937         switch(functionSpaceType_target) {        return true;
1938            case(ReducedDegreesOfFreedom):      } else {
1939            case(ReducedNodes):            return false;
1940            case(Elements):      }
1941            case(ReducedElements):     case(FaceElements):
1942            case(FaceElements):      if (functionSpaceType_target==FaceElements) {
1943            case(ReducedFaceElements):              return true;
1944            case(Points):      } else if (functionSpaceType_target==ReducedFaceElements) {
1945            case(ContactElementsZero):              return true;
1946            case(ReducedContactElementsZero):      } else {
1947            case(ContactElementsOne):              return false;
1948            case(ReducedContactElementsOne):      }
1949                return true;     case(ReducedFaceElements):
1950            case(Nodes):      if (functionSpaceType_target==ReducedFaceElements) {
1951            case(DegreesOfFreedom):              return true;
1952               return false;      } else {
1953            default:          return false;
1954               stringstream temp;      }
1955               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;     case(Points):
1956               throw FinleyAdapterException(temp.str());      if (functionSpaceType_target==Points) {
1957         }              return true;
1958         break;      } else {
1959       default:              return false;
1960          stringstream temp;      }
1961          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;     case(ContactElementsZero):
1962          throw FinleyAdapterException(temp.str());     case(ContactElementsOne):
1963          break;      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1964    }              return true;
1965    checkFinleyError();      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1966    return false;              return true;
1967        } else {
1968                return false;
1969        }
1970       case(ReducedContactElementsZero):
1971       case(ReducedContactElementsOne):
1972        if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1973                return true;
1974        } else {
1975                return false;
1976        }
1977       case(DegreesOfFreedom):
1978        switch(functionSpaceType_target) {
1979        case(ReducedDegreesOfFreedom):
1980        case(DegreesOfFreedom):
1981        case(Nodes):
1982        case(ReducedNodes):
1983        case(Elements):
1984        case(ReducedElements):
1985        case(Points):
1986        case(FaceElements):
1987        case(ReducedFaceElements):
1988        case(ContactElementsZero):
1989        case(ReducedContactElementsZero):
1990        case(ContactElementsOne):
1991        case(ReducedContactElementsOne):
1992        return true;
1993        default:
1994            stringstream temp;
1995            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1996            throw FinleyAdapterException(temp.str());
1997        }
1998        break;
1999       case(ReducedDegreesOfFreedom):
2000       switch(functionSpaceType_target) {
2001        case(ReducedDegreesOfFreedom):
2002        case(ReducedNodes):
2003        case(Elements):
2004        case(ReducedElements):
2005        case(FaceElements):
2006        case(ReducedFaceElements):
2007        case(Points):
2008        case(ContactElementsZero):
2009        case(ReducedContactElementsZero):
2010        case(ContactElementsOne):
2011        case(ReducedContactElementsOne):
2012        return true;
2013        case(Nodes):
2014        case(DegreesOfFreedom):
2015        return false;
2016        default:
2017            stringstream temp;
2018            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
2019            throw FinleyAdapterException(temp.str());
2020        }
2021        break;
2022       default:
2023          stringstream temp;
2024          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
2025          throw FinleyAdapterException(temp.str());
2026          break;
2027       }
2028       checkFinleyError();
2029       return false;
2030  }  }
2031    
2032  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
2033  {  {
2034      return false;     return false;
2035  }  }
2036    
2037  bool MeshAdapter::operator==(const AbstractDomain& other) const  bool MeshAdapter::operator==(const AbstractDomain& other) const
2038  {  {
2039    const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);     const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2040    if (temp!=0) {     if (temp!=0) {
2041      return (m_finleyMesh==temp->m_finleyMesh);        return (m_finleyMesh==temp->m_finleyMesh);
2042    } else {     } else {
2043      return false;        return false;
2044    }     }
2045  }  }
2046    
2047  bool MeshAdapter::operator!=(const AbstractDomain& other) const  bool MeshAdapter::operator!=(const AbstractDomain& other) const
2048  {  {
2049    return !(operator==(other));     return !(operator==(other));
2050  }  }
2051    
2052  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
2053  {  {
2054     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2055     checkPasoError();     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2056     return out;             package, symmetry, mesh->MPIInfo);
2057    }
2058    
2059    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2060    {
2061       Finley_Mesh* mesh=m_finleyMesh.get();
2062       return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2063               package, symmetry, mesh->MPIInfo);
2064  }  }
2065    
2066  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2067  {  {
2068    return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2069  }  }
2070    
2071  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2072  {  {
2073    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2074  }  }
2075    
2076  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2077  {  {
2078    return function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2079  }  }
2080    
2081  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2082  {  {
2083    int *out=0,i;     int *out = NULL;
2084    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2085    switch (functionSpaceType) {     switch (functionSpaceType) {
2086      case(Nodes):     case(Nodes):
2087        out=mesh->Nodes->Id;     out=mesh->Nodes->Id;
2088        break;     break;
2089      case(ReducedNodes):     case(ReducedNodes):
2090        out=mesh->Nodes->reducedNodesId;     out=mesh->Nodes->reducedNodesId;
2091        break;     break;
2092      case(Elements):     case(Elements):
2093        out=mesh->Elements->Id;     out=mesh->Elements->Id;
2094        break;     break;
2095      case(ReducedElements):     case(ReducedElements):
2096        out=mesh->Elements->Id;     out=mesh->Elements->Id;
2097        break;     break;
2098      case(FaceElements):     case(FaceElements):
2099        out=mesh->FaceElements->Id;     out=mesh->FaceElements->Id;
2100        break;     break;
2101      case(ReducedFaceElements):     case(ReducedFaceElements):
2102        out=mesh->FaceElements->Id;     out=mesh->FaceElements->Id;
2103        break;     break;
2104      case(Points):     case(Points):
2105        out=mesh->Points->Id;     out=mesh->Points->Id;
2106        break;     break;
2107      case(ContactElementsZero):     case(ContactElementsZero):
2108        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2109        break;     break;
2110      case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
2111        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2112        break;     break;
2113      case(ContactElementsOne):     case(ContactElementsOne):
2114        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2115        break;     break;
2116      case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
2117        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2118        break;     break;
2119      case(DegreesOfFreedom):     case(DegreesOfFreedom):
2120        out=mesh->Nodes->degreesOfFreedomId;     out=mesh->Nodes->degreesOfFreedomId;
2121        break;     break;
2122      case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2123        out=mesh->Nodes->reducedDegreesOfFreedomId;     out=mesh->Nodes->reducedDegreesOfFreedomId;
2124        break;     break;
2125      default:     default:
2126        stringstream temp;        stringstream temp;
2127        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2128        throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2129        break;        break;
2130    }     }
2131    return out;     return out;
2132  }  }
2133  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2134  {  {
2135    int out=0;     int out=0;
2136    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2137    switch (functionSpaceType) {     switch (functionSpaceType) {
2138    case(Nodes):     case(Nodes):
2139      out=mesh->Nodes->Tag[sampleNo];     out=mesh->Nodes->Tag[sampleNo];
2140      break;     break;
2141    case(ReducedNodes):     case(ReducedNodes):
2142      throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2143      break;     break;
2144    case(Elements):     case(Elements):
2145      out=mesh->Elements->Tag[sampleNo];     out=mesh->Elements->Tag[sampleNo];
2146      break;     break;
2147    case(ReducedElements):     case(ReducedElements):
2148      out=mesh->Elements->Tag[sampleNo];     out=mesh->Elements->Tag[sampleNo];
2149      break;     break;
2150    case(FaceElements):     case(FaceElements):
2151      out=mesh->FaceElements->Tag[sampleNo];     out=mesh->FaceElements->Tag[sampleNo];
2152      break;     break;
2153    case(ReducedFaceElements):     case(ReducedFaceElements):
2154      out=mesh->FaceElements->Tag[sampleNo];     out=mesh->FaceElements->Tag[sampleNo];
2155      break;     break;
2156    case(Points):     case(Points):
2157      out=mesh->Points->Tag[sampleNo];     out=mesh->Points->Tag[sampleNo];
2158      break;     break;
2159    case(ContactElementsZero):     case(ContactElementsZero):
2160      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2161      break;     break;
2162    case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
2163      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2164      break;     break;
2165    case(ContactElementsOne):     case(ContactElementsOne):
2166      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2167      break;     break;
2168    case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
2169      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2170      break;     break;
2171    case(DegreesOfFreedom):     case(DegreesOfFreedom):
2172      throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2173      break;     break;
2174    case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2175      throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");     throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2176      break;     break;
2177    default:     default:
2178      stringstream temp;        stringstream temp;
2179      temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2180      throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2181      break;        break;
2182    }     }
2183    return out;     return out;
2184  }  }
2185    
2186    
2187  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2188  {  {
2189       Finley_Mesh* mesh=m_finleyMesh.get();
2190       escriptDataC tmp=mask.getDataC();
2191       switch(functionSpaceType) {
2192       case(Nodes):
2193       Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2194       break;
2195       case(ReducedNodes):
2196       throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2197       break;
2198       case(DegreesOfFreedom):
2199       throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2200       break;
2201       case(ReducedDegreesOfFreedom):
2202       throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2203       break;
2204       case(Elements):
2205       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2206       break;
2207       case(ReducedElements):
2208       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2209       break;
2210       case(FaceElements):
2211       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2212       break;
2213       case(ReducedFaceElements):
2214       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2215       break;
2216       case(Points):
2217       Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2218       break;
2219       case(ContactElementsZero):
2220       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2221       break;
2222       case(ReducedContactElementsZero):
2223       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2224       break;
2225       case(ContactElementsOne):
2226       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2227       break;
2228       case(ReducedContactElementsOne):
2229       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2230       break;
2231       default:
2232          stringstream temp;
2233          temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2234          throw FinleyAdapterException(temp.str());
2235       }
2236       checkFinleyError();
2237       return;
2238    }
2239    
2240    void MeshAdapter::setTagMap(const string& name,  int tag)
2241    {
2242       Finley_Mesh* mesh=m_finleyMesh.get();
2243       Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2244       checkFinleyError();
2245       // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2246    }
2247    
2248    int MeshAdapter::getTag(const string& name) const
2249    {
2250       Finley_Mesh* mesh=m_finleyMesh.get();
2251       int tag=0;
2252       tag=Finley_Mesh_getTag(mesh, name.c_str());
2253       checkFinleyError();
2254       // throwStandardException("MeshAdapter::getTag is not implemented.");
2255       return tag;
2256    }
2257    
2258    bool MeshAdapter::isValidTagName(const string& name) const
2259    {
2260       Finley_Mesh* mesh=m_finleyMesh.get();
2261       return Finley_Mesh_isValidTagName(mesh,name.c_str());
2262    }
2263    
2264    string MeshAdapter::showTagNames() const
2265    {
2266       stringstream temp;
2267       Finley_Mesh* mesh=m_finleyMesh.get();
2268       Finley_TagMap* tag_map=mesh->TagMap;
2269       while (tag_map) {
2270          temp << tag_map->name;
2271          tag_map=tag_map->next;
2272          if (tag_map) temp << ", ";
2273       }
2274       return temp.str();
2275    }
2276    
2277    int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2278    {
2279    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2280    escriptDataC tmp=mask.getDataC();    dim_t numTags=0;
2281    switch(functionSpaceType) {    switch(functionSpaceCode) {
2282         case(Nodes):     case(Nodes):
2283            Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);            numTags=mesh->Nodes->numTagsInUse;
2284            break;            break;
2285         case(ReducedNodes):     case(ReducedNodes):
2286            throw FinleyAdapterException("Error - ReducedNodes does not support tags");            throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2287            break;            break;
2288         case(DegreesOfFreedom):     case(DegreesOfFreedom):
2289            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2290            break;            break;
2291         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2292            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2293            break;            break;
2294         case(Elements):     case(Elements):
2295            Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     case(ReducedElements):
2296            break;            numTags=mesh->Elements->numTagsInUse;
2297         case(ReducedElements):            break;
2298            Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     case(FaceElements):
2299            break;     case(ReducedFaceElements):
2300         case(FaceElements):            numTags=mesh->FaceElements->numTagsInUse;
2301            Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);            break;
2302            break;     case(Points):
2303         case(ReducedFaceElements):            numTags=mesh->Points->numTagsInUse;
2304            Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);            break;
2305       case(ContactElementsZero):
2306       case(ReducedContactElementsZero):
2307       case(ContactElementsOne):
2308       case(ReducedContactElementsOne):
2309              numTags=mesh->ContactElements->numTagsInUse;
2310            break;            break;
2311         case(Points):     default:
2312            Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);        stringstream temp;
2313          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2314          throw FinleyAdapterException(temp.str());
2315      }
2316      return numTags;
2317    }
2318    
2319    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2320    {
2321      Finley_Mesh* mesh=m_finleyMesh.get();
2322      index_t* tags=NULL;
2323      switch(functionSpaceCode) {
2324       case(Nodes):
2325              tags=mesh->Nodes->tagsInUse;
2326            break;            break;
2327         case(ContactElementsZero):     case(ReducedNodes):
2328            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2329            break;            break;
2330         case(ReducedContactElementsZero):     case(DegreesOfFreedom):
2331            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2332            break;            break;
2333         case(ContactElementsOne):     case(ReducedDegreesOfFreedom):
2334            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2335            break;            break;
2336         case(ReducedContactElementsOne):     case(Elements):
2337            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);     case(ReducedElements):
2338              tags=mesh->Elements->tagsInUse;
2339              break;
2340       case(FaceElements):
2341       case(ReducedFaceElements):
2342              tags=mesh->FaceElements->tagsInUse;
2343              break;
2344       case(Points):
2345              tags=mesh->Points->tagsInUse;
2346              break;
2347       case(ContactElementsZero):
2348       case(ReducedContactElementsZero):
2349       case(ContactElementsOne):
2350       case(ReducedContactElementsOne):
2351              tags=mesh->ContactElements->tagsInUse;
2352            break;            break;
2353         default:     default:
2354            stringstream temp;        stringstream temp;
2355            temp << "Error - Finley does not know anything about function space type " << functionSpaceType;        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2356            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2357    }    }
2358    checkFinleyError();    return tags;
   return;  
2359  }  }
2360    
2361  void MeshAdapter::setTagMap(const std::string& name,  int tag)  
2362    bool MeshAdapter::canTag(int functionSpaceCode) const
2363  {  {
2364    Finley_Mesh* mesh=m_finleyMesh.get();    switch(functionSpaceCode) {
2365    Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     case(Nodes):
2366    checkPasoError();     case(Elements):
2367    // throwStandardException("MeshAdapter::set TagMap is not implemented.");     case(ReducedElements):
2368       case(FaceElements):
2369       case(ReducedFaceElements):
2370       case(Points):
2371       case(ContactElementsZero):
2372       case(ReducedContactElementsZero):
2373       case(ContactElementsOne):
2374       case(ReducedContactElementsOne):
2375              return true;
2376       case(ReducedNodes):
2377       case(DegreesOfFreedom):
2378       case(ReducedDegreesOfFreedom):
2379          return false;
2380       default:
2381        return false;
2382      }
2383  }  }
2384    
2385  int MeshAdapter::getTag(const std::string& name) const  AbstractDomain::StatusType MeshAdapter::getStatus() const
2386  {  {
2387    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2388    int tag=0;    return Finley_Mesh_getStatus(mesh);
   tag=Finley_Mesh_getTag(mesh, name.c_str());  
   checkPasoError();  
   // throwStandardException("MeshAdapter::getTag is not implemented.");  
   return tag;  
2389  }  }
2390    
2391  bool MeshAdapter::isValidTagName(const std::string& name) const  int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2392  {  {
2393      
2394    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2395    return Finley_Mesh_isValidTagName(mesh,name.c_str());    int order =-1;
2396      switch(functionSpaceCode) {
2397       case(Nodes):
2398       case(DegreesOfFreedom):
2399              order=mesh->approximationOrder;
2400              break;
2401       case(ReducedNodes):
2402       case(ReducedDegreesOfFreedom):
2403              order=mesh->reducedApproximationOrder;
2404              break;
2405       case(Elements):
2406       case(FaceElements):
2407       case(Points):
2408       case(ContactElementsZero):
2409       case(ContactElementsOne):
2410              order=mesh->integrationOrder;
2411              break;
2412       case(ReducedElements):
2413       case(ReducedFaceElements):
2414       case(ReducedContactElementsZero):
2415       case(ReducedContactElementsOne):
2416              order=mesh->reducedIntegrationOrder;
2417              break;
2418       default:
2419          stringstream temp;
2420          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2421          throw FinleyAdapterException(temp.str());
2422      }
2423      return order;
2424  }  }
2425    
2426  std::string MeshAdapter::showTagNames() const  bool MeshAdapter::supportsContactElements() const
2427  {  {
2428    stringstream temp;    return true;
   Finley_Mesh* mesh=m_finleyMesh.get();  
   Finley_TagMap* tag_map=mesh->TagMap;  
   while (tag_map) {  
      temp << tag_map->name;  
      tag_map=tag_map->next;  
      if (tag_map) temp << ", ";  
   }  
   return temp.str();  
2429  }  }
2430    
2431    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2432    {
2433      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2434    }
2435    
2436    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2437    {
2438      Finley_ReferenceElementSet_dealloc(m_refSet);
2439    }
2440    
2441    // points will be flattened
2442    void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2443    {
2444          const int dim = getDim();
2445          int numPoints=points.size()/dim;
2446          int numTags=tags.size();
2447          Finley_Mesh* mesh=m_finleyMesh.get();
2448          
2449          if ( points.size() % dim != 0 )
2450          {
2451        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2452          }
2453          
2454          if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2455         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2456          
2457          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2458          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2459          
2460          for (int i=0;i<numPoints;++i) {
2461           points_ptr[ i * dim     ] = points[i * dim ];
2462           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2463           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2464               tags_ptr[i]=tags[i];
2465          }
2466          
2467          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2468          checkFinleyError();
2469          
2470          TMPMEMFREE(points_ptr);
2471          TMPMEMFREE(tags_ptr);
2472    }
2473    
2474    
2475    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2476    // {
2477    //       const int dim = getDim();
2478    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2479    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2480    //       Finley_Mesh* mesh=m_finleyMesh.get();
2481    //      
2482    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2483    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2484    //      
2485    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2486    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2487    //      
2488    //       for (int i=0;i<numPoints;++i) {
2489    //     int tag_id=-1;
2490    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2491    //     if  ( numComps !=   dim ) {
2492    //                stringstream temp;            
2493    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2494    //                throw FinleyAdapterException(temp.str());
2495    //     }
2496    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2497    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2498    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2499    //    
2500    //     if ( numTags > 0) {
2501    //            boost::python::extract<string> ex_str(tags[i]);
2502    //        if  ( ex_str.check() ) {
2503    //            tag_id=getTag( ex_str());
2504    //        } else {
2505    //             boost::python::extract<int> ex_int(tags[i]);
2506    //             if ( ex_int.check() ) {
2507    //                 tag_id=ex_int();
2508    //             } else {
2509    //              stringstream temp;          
2510    //                  temp << "Error - unable to extract tag for point " << i;
2511    //              throw FinleyAdapterException(temp.str());
2512    //            }
2513    //        }
2514    //     }      
2515    //            tags_ptr[i]=tag_id;
2516    //       }
2517    //      
2518    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2519    //       checkPasoError();
2520    //      
2521    //       TMPMEMFREE(points_ptr);
2522    //       TMPMEMFREE(tags_ptr);
2523    // }
2524    
2525    /*
2526    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2527    {  
2528        boost::python::list points =  boost::python::list();
2529        boost::python::list tags =  boost::python::list();
2530        points.append(point);
2531        tags.append(tag);
2532        addDiracPoints(points, tags);
2533    }
2534    */
2535    
2536    /*
2537    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2538    {
2539            boost::python::list points =   boost::python::list();
2540            boost::python::list tags =   boost::python::list();
2541            points.append(point);
2542            tags.append(tag);
2543            addDiracPoints(points, tags);
2544    }
2545    */
2546  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1346  
changed lines
  Added in v.3736

  ViewVC Help
Powered by ViewVC 1.1.26