/[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 1375 by gross, Wed Jan 9 00:15:05 2008 UTC revision 3522 by gross, Tue May 24 00:57:58 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    
15  #include "MeshAdapter.h"  #include "MeshAdapter.h"
16  #include "escript/Data.h"  #include "escript/Data.h"
# Line 20  Line 19 
19  #include <netcdfcpp.h>  #include <netcdfcpp.h>
20  #endif  #endif
21  extern "C" {  extern "C" {
22  #include "escript/blocktimer.h"  #include "esysUtils/blocktimer.h"
23  }  }
24  #include <vector>  
25    #include <boost/python/import.hpp>
26    
27  using namespace std;  using namespace std;
28  using namespace escript;  using namespace escript;
# Line 48  const int MeshAdapter::ReducedContactEle Line 48  const int MeshAdapter::ReducedContactEle
48    
49  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
50  {  {
51    setFunctionSpaceTypeNames();     setFunctionSpaceTypeNames();
52    //     //
53    // 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
54    // for us.     // for us.
55    m_finleyMesh.reset(finleyMesh,null_deleter());     m_finleyMesh.reset(finleyMesh,null_deleter());
56  }  }
57    
58  //  //
# Line 60  MeshAdapter::MeshAdapter(Finley_Mesh* fi Line 60  MeshAdapter::MeshAdapter(Finley_Mesh* fi
60  MeshAdapter::MeshAdapter(const MeshAdapter& in):  MeshAdapter::MeshAdapter(const MeshAdapter& in):
61  m_finleyMesh(in.m_finleyMesh)  m_finleyMesh(in.m_finleyMesh)
62  {  {
63    setFunctionSpaceTypeNames();     setFunctionSpaceTypeNames();
64  }  }
65    
66  MeshAdapter::~MeshAdapter()  MeshAdapter::~MeshAdapter()
67  {  {
68    //     //
69    // 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.
70    //  cout << "In MeshAdapter destructor." << endl;     //  cout << "In MeshAdapter destructor." << endl;
71    if (m_finleyMesh.unique()) {     if (m_finleyMesh.unique()) {
72      Finley_Mesh_free(m_finleyMesh.get());        Finley_Mesh_free(m_finleyMesh.get());
73    }     }
74  }  }
75    
76  int MeshAdapter::getMPISize() const  int MeshAdapter::getMPISize() const
# Line 81  int MeshAdapter::getMPIRank() const Line 81  int MeshAdapter::getMPIRank() const
81  {  {
82     return m_finleyMesh.get()->MPIInfo->rank;     return m_finleyMesh.get()->MPIInfo->rank;
83  }  }
84    void MeshAdapter::MPIBarrier() const
85    {
86    #ifdef ESYS_MPI
87       MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
88    #endif
89       return;
90    }
91    bool MeshAdapter::onMasterProcessor() const
92    {
93       return m_finleyMesh.get()->MPIInfo->rank == 0;
94    }
95    
96    
97    #ifdef ESYS_MPI
98      MPI_Comm
99    #else
100      unsigned int
101    #endif
102    MeshAdapter::getMPIComm() const
103    {
104    #ifdef ESYS_MPI
105        return m_finleyMesh->MPIInfo->comm;
106    #else
107        return 0;
108    #endif
109    }
110    
111    
112  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
113     return m_finleyMesh.get();     return m_finleyMesh.get();
114  }  }
115    
116  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const string& fileName) const
117  {  {
118    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;
119    strcpy(fName,fileName.c_str());     strcpy(fName,fileName.c_str());
120    Finley_Mesh_write(m_finleyMesh.get(),fName);     Finley_Mesh_write(m_finleyMesh.get(),fName);
121    checkFinleyError();     checkFinleyError();
122    TMPMEMFREE(fName);     TMPMEMFREE(fName);
123  }  }
124    
125  void MeshAdapter::Print_Mesh_Info(const bool full=false) const  void MeshAdapter::Print_Mesh_Info(const bool full) const
126  {  {
127    Finley_PrintMesh_Info(m_finleyMesh.get(), full);     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
128  }  }
129    
130  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const string& fileName) const
131  {  {
132  #ifdef USE_NETCDF  #ifdef USE_NETCDF
133     const NcDim* ncdims[12];     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
134     NcVar *ids, *data;     NcVar *ids;
135     int *int_ptr;     int *int_ptr;
136     Finley_Mesh *mesh = m_finleyMesh.get();     Finley_Mesh *mesh = m_finleyMesh.get();
137     Finley_TagMap* tag_map;     Finley_TagMap* tag_map;
# Line 121  void MeshAdapter::dump(const std::string Line 147  void MeshAdapter::dump(const std::string
147     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
148     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
149     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
150     char *newFileName = Paso_MPI_appendRankToFileName(strdup(fileName.c_str()), mpi_size, mpi_rank);  #ifdef ESYS_MPI
151       MPI_Status status;
152    #endif
153    
154    /* Incoming token indicates it's my turn to write */
155    #ifdef ESYS_MPI
156       if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
157    #endif
158    
159       char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
160                                                         mpi_size, mpi_rank);
161    
162     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
163     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
164     if (tag_map) {     num_Tags = 0;
165       while (tag_map) {     while (tag_map) {
166      num_Tags++;        num_Tags++;
167          tag_map=tag_map->next;        tag_map=tag_map->next;
      }  
168     }     }
169    
170     // NetCDF error handler     // NetCDF error handler
171     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
172     // Create the file.     // Create the file.
173     NcFile dataFile(newFileName, NcFile::Replace);     NcFile dataFile(newFileName, NcFile::Replace);
174       string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
175     // check if writing was successful     // check if writing was successful
176     if (!dataFile.is_valid())     if (!dataFile.is_valid())
177          throw DataException("Error - MeshAdapter::dump: opening of NetCDF file for output failed: " + *newFileName);        throw DataException(msgPrefix+"Open file for output");
178    
179     // Define dimensions (num_Elements and dim_Elements are identical, dim_Elements only appears if > 0)     // Define dimensions (num_Elements and dim_Elements are identical,
180       // dim_Elements only appears if > 0)
181     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
182          throw DataException("Error - MeshAdapter::dump: appending dimension numNodes to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numNodes)");
183     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
184          throw DataException("Error - MeshAdapter::dump: appending dimension numDim to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numDim)");
185     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)) )
186          throw DataException("Error - MeshAdapter::dump: appending dimension mpi_size to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(mpi_size)");
187     if (num_Elements>0)     if (num_Elements>0)
188        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
189          throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements)");
190     if (num_FaceElements>0)     if (num_FaceElements>0)
191        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
192          throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
193     if (num_ContactElements>0)     if (num_ContactElements>0)
194        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
195          throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements)");
196     if (num_Points>0)     if (num_Points>0)
197        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
198          throw DataException("Error - MeshAdapter::dump: appending dimension dim_Points to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Points)");
199     if (num_Elements>0)     if (num_Elements>0)
200        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
201          throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
202     if (num_FaceElements>0)     if (num_FaceElements>0)
203        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
204          throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
205     if (num_ContactElements>0)     if (num_ContactElements>0)
206        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
207          throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements_numNodes)");
208     if (num_Tags>0)     if (num_Tags>0)
209        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
210          throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Tags)");
211    
212     // Attributes: MPI size, MPI rank, Name, order, reduced_order     // Attributes: MPI size, MPI rank, Name, order, reduced_order
213     if (!dataFile.add_att("mpi_size", mpi_size) )     if (!dataFile.add_att("mpi_size", mpi_size) )
214          throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_size)");
215     if (!dataFile.add_att("mpi_rank", mpi_rank) )     if (!dataFile.add_att("mpi_rank", mpi_rank) )
216          throw DataException("Error - MeshAdapter::dump: appending mpi_rank to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_rank)");
217     if (!dataFile.add_att("Name",mesh->Name) )     if (!dataFile.add_att("Name",mesh->Name) )
218          throw DataException("Error - MeshAdapter::dump: appending Name to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Name)");
219     if (!dataFile.add_att("numDim",numDim) )     if (!dataFile.add_att("numDim",numDim) )
220          throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
221     if (!dataFile.add_att("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
222          throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
223     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
224          throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(reduced_order)");
225     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
226          throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(numNodes)");
227     if (!dataFile.add_att("num_Elements",num_Elements) )     if (!dataFile.add_att("num_Elements",num_Elements) )
228          throw DataException("Error - MeshAdapter::dump: appending num_Elements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements)");
229     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
230          throw DataException("Error - MeshAdapter::dump: appending num_FaceElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements)");
231     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
232          throw DataException("Error - MeshAdapter::dump: appending num_ContactElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements)");
233     if (!dataFile.add_att("num_Points",num_Points) )     if (!dataFile.add_att("num_Points",num_Points) )
234          throw DataException("Error - MeshAdapter::dump: appending num_Points to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Points)");
235     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
236          throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
237     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
238          throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
239     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
240          throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements_numNodes)");
241     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
242        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Elements_TypeId)");
243     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
244        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
245     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
246        throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(ContactElements_TypeId)");
247     if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
248        throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Points_TypeId)");
249     if (!dataFile.add_att("num_Tags", num_Tags) )     if (!dataFile.add_att("num_Tags", num_Tags) )
250        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Tags)");
251    
252     // // // // // Nodes // // // // //     // // // // // Nodes // // // // //
253    
254     // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)     // Nodes nodeDistribution
255     if (numNodes>0) {     if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
256          throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
257     // Nodes Id     int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
258     if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )     if (! (ids->put(int_ptr, mpi_size+1)) )
259          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);  
260    
261     // Nodes degreesOfFreedomDistribution     // Nodes degreesOfFreedomDistribution
262     if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )     if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
263          throw DataException("Error - MeshAdapter::dump: appending Nodes_DofDistribution to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
264     int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];     int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
265     if (! (ids->put(int_ptr, mpi_size+1)) )     if (! (ids->put(int_ptr, mpi_size+1)) )
266          throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);        throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
267    
268       // Only write nodes if non-empty because NetCDF doesn't like empty arrays
269       // (it treats them as NC_UNLIMITED)
270       if (numNodes>0) {
271    
272          // Nodes Id
273          if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
274             throw DataException(msgPrefix+"add_var(Nodes_Id)");
275          int_ptr = &mesh->Nodes->Id[0];
276          if (! (ids->put(int_ptr, numNodes)) )
277             throw DataException(msgPrefix+"put(Nodes_Id)");
278    
279          // Nodes Tag
280          if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
281             throw DataException(msgPrefix+"add_var(Nodes_Tag)");
282          int_ptr = &mesh->Nodes->Tag[0];
283          if (! (ids->put(int_ptr, numNodes)) )
284             throw DataException(msgPrefix+"put(Nodes_Tag)");
285    
286          // Nodes gDOF
287          if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
288             throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
289          int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
290          if (! (ids->put(int_ptr, numNodes)) )
291             throw DataException(msgPrefix+"put(Nodes_gDOF)");
292    
293          // Nodes global node index
294          if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
295             throw DataException(msgPrefix+"add_var(Nodes_gNI)");
296          int_ptr = &mesh->Nodes->globalNodesIndex[0];
297          if (! (ids->put(int_ptr, numNodes)) )
298             throw DataException(msgPrefix+"put(Nodes_gNI)");
299    
300          // Nodes grDof
301          if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
302             throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
303          int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
304          if (! (ids->put(int_ptr, numNodes)) )
305             throw DataException(msgPrefix+"put(Nodes_grDfI)");
306    
307          // Nodes grNI
308          if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
309             throw DataException(msgPrefix+"add_var(Nodes_grNI)");
310          int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
311          if (! (ids->put(int_ptr, numNodes)) )
312             throw DataException(msgPrefix+"put(Nodes_grNI)");
313    
314          // Nodes Coordinates
315          if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
316             throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
317          if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
318             throw DataException(msgPrefix+"put(Nodes_Coordinates)");
319    
320     }     }
321    
# Line 278  void MeshAdapter::dump(const std::string Line 323  void MeshAdapter::dump(const std::string
323    
324     if (num_Elements>0) {     if (num_Elements>0) {
325    
326       // Temp storage to gather node IDs        // Elements_Id
327       int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
328             throw DataException(msgPrefix+"add_var(Elements_Id)");
329       // Elements_Id        int_ptr = &mesh->Elements->Id[0];
330       if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! (ids->put(int_ptr, num_Elements)) )
331          throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Id)");
332       int_ptr = &mesh->Elements->Id[0];  
333       if (! (ids->put(int_ptr, num_Elements)) )        // Elements_Tag
334          throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
335             throw DataException(msgPrefix+"add_var(Elements_Tag)");
336       // Elements_Tag        int_ptr = &mesh->Elements->Tag[0];
337       if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )        if (! (ids->put(int_ptr, num_Elements)) )
338          throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Tag)");
339       int_ptr = &mesh->Elements->Tag[0];  
340       if (! (ids->put(int_ptr, num_Elements)) )        // Elements_Owner
341          throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
342             throw DataException(msgPrefix+"add_var(Elements_Owner)");
343       // Elements_Owner        int_ptr = &mesh->Elements->Owner[0];
344       if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )        if (! (ids->put(int_ptr, num_Elements)) )
345          throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Owner)");
346       int_ptr = &mesh->Elements->Owner[0];  
347       if (! (ids->put(int_ptr, num_Elements)) )        // Elements_Color
348          throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
349             throw DataException(msgPrefix+"add_var(Elements_Color)");
350       // Elements_Color        int_ptr = &mesh->Elements->Color[0];
351       if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )        if (! (ids->put(int_ptr, num_Elements)) )
352          throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Color)");
353       int_ptr = &mesh->Elements->Color[0];  
354       if (! (ids->put(int_ptr, num_Elements)) )        // Elements_Nodes
355          throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
356             throw DataException(msgPrefix+"add_var(Elements_Nodes)");
357       // Elements_Nodes        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
358       for (int i=0; i<num_Elements; i++)           throw DataException(msgPrefix+"put(Elements_Nodes)");
        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);  
359    
360     }     }
361    
# Line 326  void MeshAdapter::dump(const std::string Line 363  void MeshAdapter::dump(const std::string
363    
364     if (num_FaceElements>0) {     if (num_FaceElements>0) {
365    
366       // Temp storage to gather node IDs        // FaceElements_Id
367       int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
368             throw DataException(msgPrefix+"add_var(FaceElements_Id)");
369       // FaceElements_Id        int_ptr = &mesh->FaceElements->Id[0];
370       if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! (ids->put(int_ptr, num_FaceElements)) )
371          throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Id)");
372       int_ptr = &mesh->FaceElements->Id[0];  
373       if (! (ids->put(int_ptr, num_FaceElements)) )        // FaceElements_Tag
374          throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
375             throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
376       // FaceElements_Tag        int_ptr = &mesh->FaceElements->Tag[0];
377       if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )        if (! (ids->put(int_ptr, num_FaceElements)) )
378          throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Tag)");
379       int_ptr = &mesh->FaceElements->Tag[0];  
380       if (! (ids->put(int_ptr, num_FaceElements)) )        // FaceElements_Owner
381          throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
382             throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
383       // FaceElements_Owner        int_ptr = &mesh->FaceElements->Owner[0];
384       if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )        if (! (ids->put(int_ptr, num_FaceElements)) )
385          throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Owner)");
386       int_ptr = &mesh->FaceElements->Owner[0];  
387       if (! (ids->put(int_ptr, num_FaceElements)) )        // FaceElements_Color
388          throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
389             throw DataException(msgPrefix+"add_var(FaceElements_Color)");
390       // FaceElements_Color        int_ptr = &mesh->FaceElements->Color[0];
391       if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )        if (! (ids->put(int_ptr, num_FaceElements)) )
392          throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Color)");
393       int_ptr = &mesh->FaceElements->Color[0];  
394       if (! (ids->put(int_ptr, num_FaceElements)) )        // FaceElements_Nodes
395          throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
396             throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
397       // FaceElements_Nodes        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
398       for (int i=0; i<num_FaceElements; i++)           throw DataException(msgPrefix+"put(FaceElements_Nodes)");
        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);  
399    
400     }     }
401    
# Line 374  void MeshAdapter::dump(const std::string Line 403  void MeshAdapter::dump(const std::string
403    
404     if (num_ContactElements>0) {     if (num_ContactElements>0) {
405    
406       // Temp storage to gather node IDs        // ContactElements_Id
407       int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
408             throw DataException(msgPrefix+"add_var(ContactElements_Id)");
409       // ContactElements_Id        int_ptr = &mesh->ContactElements->Id[0];
410       if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )        if (! (ids->put(int_ptr, num_ContactElements)) )
411          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Id)");
412       int_ptr = &mesh->ContactElements->Id[0];  
413       if (! (ids->put(int_ptr, num_ContactElements)) )        // ContactElements_Tag
414          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
415             throw DataException(msgPrefix+"add_var(ContactElements_Tag)");
416       // ContactElements_Tag        int_ptr = &mesh->ContactElements->Tag[0];
417       if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )        if (! (ids->put(int_ptr, num_ContactElements)) )
418          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Tag)");
419       int_ptr = &mesh->ContactElements->Tag[0];  
420       if (! (ids->put(int_ptr, num_ContactElements)) )        // ContactElements_Owner
421          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
422             throw DataException(msgPrefix+"add_var(ContactElements_Owner)");
423       // ContactElements_Owner        int_ptr = &mesh->ContactElements->Owner[0];
424       if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )        if (! (ids->put(int_ptr, num_ContactElements)) )
425          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Owner)");
426       int_ptr = &mesh->ContactElements->Owner[0];  
427       if (! (ids->put(int_ptr, num_ContactElements)) )        // ContactElements_Color
428          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
429             throw DataException(msgPrefix+"add_var(ContactElements_Color)");
430       // ContactElements_Color        int_ptr = &mesh->ContactElements->Color[0];
431       if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )        if (! (ids->put(int_ptr, num_ContactElements)) )
432          throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Color)");
433       int_ptr = &mesh->ContactElements->Color[0];  
434       if (! (ids->put(int_ptr, num_ContactElements)) )        // ContactElements_Nodes
435          throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
436             throw DataException(msgPrefix+"add_var(ContactElements_Nodes)");
437       // ContactElements_Nodes        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
438       for (int i=0; i<num_ContactElements; i++)           throw DataException(msgPrefix+"put(ContactElements_Nodes)");
        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);  
439    
440     }     }
441    
# Line 422  void MeshAdapter::dump(const std::string Line 443  void MeshAdapter::dump(const std::string
443    
444     if (num_Points>0) {     if (num_Points>0) {
445    
446       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_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);  
447    
448       TMPMEMFREE(Points_Nodes);        // Points_Id
449          if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
450             throw DataException(msgPrefix+"add_var(Points_Id)");
451          int_ptr = &mesh->Points->Id[0];
452          if (! (ids->put(int_ptr, num_Points)) )
453             throw DataException(msgPrefix+"put(Points_Id)");
454    
455          // Points_Tag
456          if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
457             throw DataException(msgPrefix+"add_var(Points_Tag)");
458          int_ptr = &mesh->Points->Tag[0];
459          if (! (ids->put(int_ptr, num_Points)) )
460             throw DataException(msgPrefix+"put(Points_Tag)");
461    
462          // Points_Owner
463          if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
464             throw DataException(msgPrefix+"add_var(Points_Owner)");
465          int_ptr = &mesh->Points->Owner[0];
466          if (! (ids->put(int_ptr, num_Points)) )
467             throw DataException(msgPrefix+"put(Points_Owner)");
468    
469          // Points_Color
470          if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
471             throw DataException(msgPrefix+"add_var(Points_Color)");
472          int_ptr = &mesh->Points->Color[0];
473          if (! (ids->put(int_ptr, num_Points)) )
474             throw DataException(msgPrefix+"put(Points_Color)");
475    
476          // Points_Nodes
477          // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
478          if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
479             throw DataException(msgPrefix+"add_var(Points_Nodes)");
480          if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
481             throw DataException(msgPrefix+"put(Points_Nodes)");
482    
483     }     }
484    
# Line 472  void MeshAdapter::dump(const std::string Line 486  void MeshAdapter::dump(const std::string
486    
487     if (num_Tags>0) {     if (num_Tags>0) {
488    
489       // Temp storage to gather node IDs        // Temp storage to gather node IDs
490       int *Tags_keys = TMPMEMALLOC(num_Tags, int);        int *Tags_keys = TMPMEMALLOC(num_Tags, int);
491       char name_temp[4096];        char name_temp[4096];
492    
493       /* Copy tag data into temp arrays */        /* Copy tag data into temp arrays */
494       tag_map = mesh->TagMap;        tag_map = mesh->TagMap;
495       if (tag_map) {        if (tag_map) {
496         int i = 0;           int i = 0;
497         while (tag_map) {           while (tag_map) {
498      Tags_keys[i++] = tag_map->tag_key;              Tags_keys[i++] = tag_map->tag_key;
499          tag_map=tag_map->next;              tag_map=tag_map->next;
500         }           }
501       }        }
   
      // Tags_keys  
      if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )  
         throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);  
      int_ptr = &Tags_keys[0];  
      if (! (ids->put(int_ptr, num_Tags)) )  
         throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);  
   
      // Tags_names_*  
      // This is an array of strings, it should be stored as an array but instead I have hacked in one attribute per string  
      // because the NetCDF manual doesn't tell how to do an array of strings  
      tag_map = mesh->TagMap;  
      if (tag_map) {  
        int i = 0;  
        while (tag_map) {  
          sprintf(name_temp, "Tags_name_%d", i);  
          if (!dataFile.add_att(name_temp, tag_map->name) )  
            throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);  
          tag_map=tag_map->next;  
      i++;  
        }  
      }  
502    
503       TMPMEMFREE(Tags_keys);        // Tags_keys
504          if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
505             throw DataException(msgPrefix+"add_var(Tags_keys)");
506          int_ptr = &Tags_keys[0];
507          if (! (ids->put(int_ptr, num_Tags)) )
508             throw DataException(msgPrefix+"put(Tags_keys)");
509    
510          // Tags_names_*
511          // This is an array of strings, it should be stored as an array but
512          // instead I have hacked in one attribute per string because the NetCDF
513          // manual doesn't tell how to do an array of strings
514          tag_map = mesh->TagMap;
515          if (tag_map) {
516             int i = 0;
517             while (tag_map) {
518                sprintf(name_temp, "Tags_name_%d", i);
519                if (!dataFile.add_att(name_temp, tag_map->name) )
520                   throw DataException(msgPrefix+"add_att(Tags_names_XX)");
521                tag_map=tag_map->next;
522                i++;
523             }
524          }
525    
526          TMPMEMFREE(Tags_keys);
527     }     }
528    
529    /* Send token to next MPI process so he can take his turn */
530    #ifdef ESYS_MPI
531       if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
532    #endif
533    
534     // NetCDF file is closed by destructor of NcFile object     // NetCDF file is closed by destructor of NcFile object
535    
536  #else  #else
537     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.");
538  #endif  /* USE_NETCDF */  #endif  /* USE_NETCDF */
# Line 522  void MeshAdapter::dump(const std::string Line 541  void MeshAdapter::dump(const std::string
541    
542  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
543  {  {
544    return "FinleyMesh";     return "FinleyMesh";
545  }  }
546    
547  string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const  string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
548  {  {
549    FunctionSpaceNamesMapType::iterator loc;     FunctionSpaceNamesMapType::iterator loc;
550    loc=m_functionSpaceTypeNames.find(functionSpaceType);     loc=m_functionSpaceTypeNames.find(functionSpaceType);
551    if (loc==m_functionSpaceTypeNames.end()) {     if (loc==m_functionSpaceTypeNames.end()) {
552      return "Invalid function space type code.";        return "Invalid function space type code.";
553    } else {     } else {
554      return loc->second;        return loc->second;
555    }     }
556  }  }
557    
558  bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const  bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
559  {  {
560    FunctionSpaceNamesMapType::iterator loc;     FunctionSpaceNamesMapType::iterator loc;
561    loc=m_functionSpaceTypeNames.find(functionSpaceType);     loc=m_functionSpaceTypeNames.find(functionSpaceType);
562    return (loc!=m_functionSpaceTypeNames.end());     return (loc!=m_functionSpaceTypeNames.end());
563  }  }
564    
565  void MeshAdapter::setFunctionSpaceTypeNames()  void MeshAdapter::setFunctionSpaceTypeNames()
566  {  {
567    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
568      (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));     (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));
569    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
570      (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));     (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));
571    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
572      (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));     (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
573    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
574      (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));     (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));
575    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
576      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));     (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
577    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
578      (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));     (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
579    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
580      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));     (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
581    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
582      (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));     (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
583    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
584      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));     (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
585    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
586      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));     (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
587    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
588      (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));     (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
589    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
590      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));     (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
591    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
592      (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));     (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
593  }  }
594    
595  int MeshAdapter::getContinuousFunctionCode() const  int MeshAdapter::getContinuousFunctionCode() const
596  {  {
597    return Nodes;     return Nodes;
598  }  }
599  int MeshAdapter::getReducedContinuousFunctionCode() const  int MeshAdapter::getReducedContinuousFunctionCode() const
600  {  {
601    return ReducedNodes;     return ReducedNodes;
602  }  }
603    
604  int MeshAdapter::getFunctionCode() const  int MeshAdapter::getFunctionCode() const
605  {  {
606    return Elements;     return Elements;
607  }  }
608  int MeshAdapter::getReducedFunctionCode() const  int MeshAdapter::getReducedFunctionCode() const
609  {  {
610    return ReducedElements;     return ReducedElements;
611  }  }
612    
613  int MeshAdapter::getFunctionOnBoundaryCode() const  int MeshAdapter::getFunctionOnBoundaryCode() const
614  {  {
615    return FaceElements;     return FaceElements;
616  }  }
617  int MeshAdapter::getReducedFunctionOnBoundaryCode() const  int MeshAdapter::getReducedFunctionOnBoundaryCode() const
618  {  {
619    return ReducedFaceElements;     return ReducedFaceElements;
620  }  }
621    
622  int MeshAdapter::getFunctionOnContactZeroCode() const  int MeshAdapter::getFunctionOnContactZeroCode() const
623  {  {
624    return ContactElementsZero;     return ContactElementsZero;
625  }  }
626  int MeshAdapter::getReducedFunctionOnContactZeroCode() const  int MeshAdapter::getReducedFunctionOnContactZeroCode() const
627  {  {
628    return ReducedContactElementsZero;     return ReducedContactElementsZero;
629  }  }
630    
631  int MeshAdapter::getFunctionOnContactOneCode() const  int MeshAdapter::getFunctionOnContactOneCode() const
632  {  {
633    return ContactElementsOne;     return ContactElementsOne;
634  }  }
635  int MeshAdapter::getReducedFunctionOnContactOneCode() const  int MeshAdapter::getReducedFunctionOnContactOneCode() const
636  {  {
637    return ReducedContactElementsOne;     return ReducedContactElementsOne;
638  }  }
639    
640  int MeshAdapter::getSolutionCode() const  int MeshAdapter::getSolutionCode() const
641  {  {
642    return DegreesOfFreedom;     return DegreesOfFreedom;
643  }  }
644    
645  int MeshAdapter::getReducedSolutionCode() const  int MeshAdapter::getReducedSolutionCode() const
646  {  {
647    return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
648  }  }
649    
650  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
651  {  {
652    return Points;     return Points;
653  }  }
654    
655  //  //
# Line 638  int MeshAdapter::getDiracDeltaFunctionCo Line 657  int MeshAdapter::getDiracDeltaFunctionCo
657  //  //
658  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
659  {  {
660    int numDim=Finley_Mesh_getDim(m_finleyMesh.get());     Finley_Mesh* mesh=m_finleyMesh.get();
661    checkFinleyError();     int numDim=Finley_Mesh_getDim(mesh);
662    return numDim;     checkFinleyError();
663       return numDim;
664    }
665    
666    //
667    // Return the number of data points summed across all MPI processes
668    //
669    int MeshAdapter::getNumDataPointsGlobal() const
670    {
671       return Finley_NodeFile_getGlobalNumNodes(m_finleyMesh.get()->Nodes);
672  }  }
673    
674  //  //
# Line 653  pair<int,int> MeshAdapter::getDataShape( Line 681  pair<int,int> MeshAdapter::getDataShape(
681     int numSamples=0;     int numSamples=0;
682     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
683     switch (functionSpaceCode) {     switch (functionSpaceCode) {
684        case(Nodes):     case(Nodes):
685             numDataPointsPerSample=1;     numDataPointsPerSample=1;
686             numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);     numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);
687             break;     break;
688        case(ReducedNodes):     case(ReducedNodes):
689             numDataPointsPerSample=1;     numDataPointsPerSample=1;
690             numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);     numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
691             break;     break;
692        case(Elements):     case(Elements):
693             if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
694               numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
695               numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
696             }     }
697             break;     break;
698        case(ReducedElements):     case(ReducedElements):
699             if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
700               numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
701               numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
702             }     }
703             break;     break;
704        case(FaceElements):     case(FaceElements):
705             if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
706                  numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
707                  numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
708             }     }
709             break;     break;
710        case(ReducedFaceElements):     case(ReducedFaceElements):
711             if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
712                  numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
713                  numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
714             }     }
715             break;     break;
716        case(Points):     case(Points):
717             if (mesh->Points!=NULL) {     if (mesh->Points!=NULL) {
718               numDataPointsPerSample=1;        numDataPointsPerSample=1;
719               numSamples=mesh->Points->numElements;        numSamples=mesh->Points->numElements;
720             }     }
721             break;     break;
722        case(ContactElementsZero):     case(ContactElementsZero):
723             if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
724               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
725               numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
726             }     }
727             break;     break;
728        case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
729             if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
730               numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
731               numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
732             }     }
733             break;     break;
734        case(ContactElementsOne):     case(ContactElementsOne):
735             if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
736               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
737               numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
738             }     }
739             break;     break;
740        case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
741             if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
742               numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
743               numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
744             }     }
745             break;     break;
746        case(DegreesOfFreedom):     case(DegreesOfFreedom):
747             if (mesh->Nodes!=NULL) {     if (mesh->Nodes!=NULL) {
748               numDataPointsPerSample=1;        numDataPointsPerSample=1;
749               numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);        numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
750             }     }
751             break;     break;
752        case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
753             if (mesh->Nodes!=NULL) {     if (mesh->Nodes!=NULL) {
754               numDataPointsPerSample=1;        numDataPointsPerSample=1;
755               numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);        numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
756             }     }
757             break;     break;
758        default:     default:
759          stringstream temp;        stringstream temp;
760          temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
761          throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
762          break;        break;
763        }     }
764        return pair<int,int>(numDataPointsPerSample,numSamples);     return pair<int,int>(numDataPointsPerSample,numSamples);
765  }  }
766    
767  //  //
768  // 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:
769  //  //
770  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
771                       SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
772                       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,
773                       const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
774                       const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact,
775  {                                   const escript::Data& d_dirac,const escript::Data& y_dirac) const
776    {
777       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
778       if (smat==0)
779       {
780        throw FinleyAdapterException("finley only supports its own system matrices.");
781       }
782     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
783     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
784     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 756  void MeshAdapter::addPDEToSystem( Line 790  void MeshAdapter::addPDEToSystem(
790     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
791     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
792     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
793       escriptDataC _d_dirac=d_dirac.getDataC();
794       escriptDataC _y_dirac=y_dirac.getDataC();
795    
796     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
797    
798     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 );
799     checkFinleyError();     checkFinleyError();
800    
801     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
802     checkFinleyError();     checkFinleyError();
803    
804     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
805       checkFinleyError();
806    
807        Finley_Assemble_PDE(mesh->Nodes,mesh->Points, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
808     checkFinleyError();     checkFinleyError();
809  }  }
810    
811  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
812                       escript::Data& mat,                                          escript::Data& mat,
813                       const escript::Data& D,                                          const escript::Data& D,
814                       const escript::Data& d) const                                          const escript::Data& d,
815                                            const escript::Data& d_dirac,
816                                            const bool useHRZ) const
817  {  {
818     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
819     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
820     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
821       escriptDataC _d_dirac=d_dirac.getDataC();
822    
823     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
824    
825     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
826     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     checkFinleyError();
827      
828       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
829       checkFinleyError();
830    
831       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Points,&_mat, &_d_dirac, useHRZ);
832     checkFinleyError();     checkFinleyError();
833    
834  }  }
835    
836    
837  //  //
838  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
839  //  //
840  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
841  {  {
842     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
843    
# Line 799  void MeshAdapter::addPDEToRHS( escript:: Line 846  void MeshAdapter::addPDEToRHS( escript::
846     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
847     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
848     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
849       escriptDataC _y_dirac=y_dirac.getDataC();
850    
851     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 );
852     checkFinleyError();     checkFinleyError();
# Line 808  void MeshAdapter::addPDEToRHS( escript:: Line 856  void MeshAdapter::addPDEToRHS( escript::
856    
857     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 );
858     checkFinleyError();     checkFinleyError();
859    
860       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs , 0, 0, 0, 0, 0, &_y_dirac );
861       checkFinleyError();
862    
863  }  }
864  //  //
865  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
866  //  //
867  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
868                       TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
869                       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,
870                       const escript::Data& d, const escript::Data& y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
871                       const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d, const escript::Data& y,
872  {                                             const escript::Data& d_contact,const escript::Data& y_contact,
873     DataArrayView::ShapeType shape;                                             const escript::Data& d_dirac, const escript::Data& y_dirac) const
874    {
875       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
876       if (tpa==0)
877       {
878        throw FinleyAdapterException("finley only supports its own transport problems.");
879       }
880    
881    
882       DataTypes::ShapeType shape;
883     source.expand();     source.expand();
    escript:: Data tmp(0.0,M.getDataPointShape(),tp.getFunctionSpace(),true);  
884     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
    escriptDataC _tmp=tmp.getDataC();  
885     escriptDataC _M=M.getDataC();     escriptDataC _M=M.getDataC();
886     escriptDataC _A=A.getDataC();     escriptDataC _A=A.getDataC();
887     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 834  void MeshAdapter::addPDEToTransportProbl Line 893  void MeshAdapter::addPDEToTransportProbl
893     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
894     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
895     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
896       escriptDataC _d_dirac=d_dirac.getDataC();
897       escriptDataC _y_dirac=y_dirac.getDataC();
898    
    Finley_Mesh* mesh=m_finleyMesh.get();  
    Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();  
     
899    
900     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_tmp, &_M);     Finley_Mesh* mesh=m_finleyMesh.get();
901     checkFinleyError();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
    /* add mass matix to lumped mass matrix of transport problem */  
    double* tmp_prt=getSampleData(&_tmp,0);  
    int i;  
    int n=Paso_FCTransportProblem_getTotalNumRows(_tp);  
    #pragma omp parallel for private(i) schedule(static)  
    for (i=0;i<n ;++i) _tp->lumped_mass_matrix[i]+=tmp_prt[i];  
902    
903     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, 0, 0, &_D, &_X, &_Y );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
904     checkFinleyError();     checkFinleyError();
905    
906     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->flux_matrix, &_source, 0, &_B, &_C, 0, 0, 0 );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
907     checkFinleyError();     checkFinleyError();
908    
909     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
# Line 859  void MeshAdapter::addPDEToTransportProbl Line 911  void MeshAdapter::addPDEToTransportProbl
911    
912     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
913     checkFinleyError();     checkFinleyError();
914    
915       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
916       checkFinleyError();
917    
918  }  }
919    
920  //  //
# Line 866  void MeshAdapter::addPDEToTransportProbl Line 922  void MeshAdapter::addPDEToTransportProbl
922  //  //
923  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
924  {  {
925    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
926    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
927    if (inDomain!=*this)       if (inDomain!=*this)  
928      throw FinleyAdapterException("Error - Illegal domain of interpolant.");        throw FinleyAdapterException("Error - Illegal domain of interpolant.");
929    if (targetDomain!=*this)     if (targetDomain!=*this)
930      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");        throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
931    
932    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
933    escriptDataC _target=target.getDataC();     escriptDataC _target=target.getDataC();
934    escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
935    switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
936       case(Nodes):     case(Nodes):
937          switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
938             case(Nodes):        case(Nodes):
939             case(ReducedNodes):        case(ReducedNodes):
940             case(DegreesOfFreedom):        case(DegreesOfFreedom):
941             case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
942                 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
943                 break;        break;
944             case(Elements):        case(Elements):
945             case(ReducedElements):        case(ReducedElements):
946                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
947                 break;        break;
948             case(FaceElements):        case(FaceElements):
949             case(ReducedFaceElements):        case(ReducedFaceElements):
950                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
951                 break;        break;
952             case(Points):        case(Points):
953                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
954                 break;        break;
955             case(ContactElementsZero):        case(ContactElementsZero):
956             case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
957                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
958                 break;        break;
959             case(ContactElementsOne):        case(ContactElementsOne):
960             case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
961                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
962                 break;        break;
963             default:        default:
964                 stringstream temp;           stringstream temp;
965                 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();
966                 throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
967                 break;           break;
968          }        }
969          break;        break;
970       case(ReducedNodes):     case(ReducedNodes):
971          switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
972             case(Nodes):        case(Nodes):
973             case(ReducedNodes):        case(ReducedNodes):
974             case(DegreesOfFreedom):        case(DegreesOfFreedom):
975             case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
976                 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
977                 break;        break;
978             case(Elements):        case(Elements):
979             case(ReducedElements):        case(ReducedElements):
980                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
981                 break;        break;
982             case(FaceElements):        case(FaceElements):
983             case(ReducedFaceElements):        case(ReducedFaceElements):
984                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
985                 break;        break;
986             case(Points):        case(Points):
987                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
988                 break;        break;
989             case(ContactElementsZero):        case(ContactElementsZero):
990             case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
991                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
992                 break;        break;
993             case(ContactElementsOne):        case(ContactElementsOne):
994             case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
995                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
996                 break;        break;
997             default:        default:
998                 stringstream temp;           stringstream temp;
999                 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();
1000                 throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
1001                 break;           break;
1002          }        }
1003          break;        break;
1004       case(Elements):     case(Elements):
1005          if (target.getFunctionSpace().getTypeCode()==Elements) {        if (target.getFunctionSpace().getTypeCode()==Elements) {
1006             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1007          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1008             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
1009          } else {        } else {
1010             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
1011          }        }
1012          break;        break;
1013       case(ReducedElements):     case(ReducedElements):
1014          if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1015             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1016          } else {        } else {
1017             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.");
1018          }        }
1019          break;        break;
1020       case(FaceElements):     case(FaceElements):
1021          if (target.getFunctionSpace().getTypeCode()==FaceElements) {        if (target.getFunctionSpace().getTypeCode()==FaceElements) {
1022             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1023          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1024             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1025          } else {        } else {
1026             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1027          }        }
1028          break;        break;
1029       case(ReducedFaceElements):     case(ReducedFaceElements):
1030          if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1031             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1032          } else {        } else {
1033             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.");
1034         }        }
1035         break;        break;
1036       case(Points):     case(Points):
1037          if (target.getFunctionSpace().getTypeCode()==Points) {        if (target.getFunctionSpace().getTypeCode()==Points) {
1038             Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1039          } else {        } else {
1040             throw FinleyAdapterException("Error - No interpolation with data on points possible.");           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1041          }        }
1042          break;        break;
1043       case(ContactElementsZero):     case(ContactElementsZero):
1044       case(ContactElementsOne):     case(ContactElementsOne):
1045          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {        if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1046             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1047          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1048             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1049          } else {        } else {
1050             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1051          }        }
1052          break;        break;
1053       case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1054       case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1055          if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1056             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1057          } else {        } else {
1058             throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1059          }        }
1060          break;        break;
1061       case(DegreesOfFreedom):           case(DegreesOfFreedom):      
1062          switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
1063             case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
1064             case(DegreesOfFreedom):        case(DegreesOfFreedom):
1065                Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1066                break;        break;
1067      
1068             case(Nodes):        case(Nodes):
1069             case(ReducedNodes):        case(ReducedNodes):
1070                if (getMPISize()>1) {        if (getMPISize()>1) {
1071                    escript::Data temp=escript::Data(in);           escript::Data temp=escript::Data(in);
1072                    temp.expand();           temp.expand();
1073                    escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1074                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1075                } else {        } else {
1076                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1077                }        }
1078                break;        break;
1079             case(Elements):        case(Elements):
1080             case(ReducedElements):        case(ReducedElements):
1081                if (getMPISize()>1) {        if (getMPISize()>1) {
1082                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1083                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1084                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1085                } else {        } else {
1086                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1087                }        }
1088                break;        break;
1089             case(FaceElements):        case(FaceElements):
1090             case(ReducedFaceElements):        case(ReducedFaceElements):
1091                if (getMPISize()>1) {        if (getMPISize()>1) {
1092                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1093                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1094                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1095      
1096                } else {        } else {
1097                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1098                }        }
1099                break;        break;
1100             case(Points):        case(Points):
1101                if (getMPISize()>1) {        if (getMPISize()>1) {
1102                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1103                   escriptDataC _in2 = temp.getDataC();           //escriptDataC _in2 = temp.getDataC();
1104                } else {        } else {
1105                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1106                }        }
1107                break;        break;
1108             case(ContactElementsZero):        case(ContactElementsZero):
1109             case(ContactElementsOne):        case(ContactElementsOne):
1110             case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1111             case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1112                if (getMPISize()>1) {        if (getMPISize()>1) {
1113                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1114                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1115                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1116                } else {        } else {
1117                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1118                }        }
1119               break;        break;
1120             default:        default:
1121               stringstream temp;           stringstream temp;
1122               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();
1123               throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
1124               break;           break;
1125          }        }
1126          break;        break;
1127       case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1128         switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
1129            case(Nodes):        case(Nodes):
1130               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1131               break;        break;
1132            case(ReducedNodes):        case(ReducedNodes):
1133                if (getMPISize()>1) {        if (getMPISize()>1) {
1134                    escript::Data temp=escript::Data(in);           escript::Data temp=escript::Data(in);
1135                    temp.expand();           temp.expand();
1136                    escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1137                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1138                } else {        } else {
1139                    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1140                }        }
1141                break;        break;
1142            case(DegreesOfFreedom):        case(DegreesOfFreedom):
1143               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1144               break;        break;
1145            case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
1146               Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1147               break;        break;
1148            case(Elements):        case(Elements):
1149            case(ReducedElements):        case(ReducedElements):
1150                if (getMPISize()>1) {        if (getMPISize()>1) {
1151                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1152                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1153                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1154                } else {        } else {
1155                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1156               }        }
1157               break;        break;
1158            case(FaceElements):        case(FaceElements):
1159            case(ReducedFaceElements):        case(ReducedFaceElements):
1160                if (getMPISize()>1) {        if (getMPISize()>1) {
1161                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1162                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1163                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1164                } else {        } else {
1165                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1166                }        }
1167               break;        break;
1168            case(Points):        case(Points):
1169                if (getMPISize()>1) {        if (getMPISize()>1) {
1170                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1171                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1172                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1173                } else {        } else {
1174                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1175               }        }
1176               break;        break;
1177            case(ContactElementsZero):        case(ContactElementsZero):
1178            case(ContactElementsOne):        case(ContactElementsOne):
1179            case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1180            case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1181                if (getMPISize()>1) {        if (getMPISize()>1) {
1182                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1183                   escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1184                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1185                } else {        } else {
1186                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1187               }        }
1188               break;        break;
1189            default:        default:
1190               stringstream temp;           stringstream temp;
1191               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();
1192               throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
1193               break;           break;
1194         }        }
1195         break;        break;
1196       default:     default:
1197          stringstream temp;        stringstream temp;
1198          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1199          throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1200          break;        break;
1201    }     }
1202    checkFinleyError();     checkFinleyError();
1203  }  }
1204    
1205  //  //
# Line 1151  void MeshAdapter::interpolateOnDomain(es Line 1207  void MeshAdapter::interpolateOnDomain(es
1207  //  //
1208  void MeshAdapter::setToX(escript::Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
1209  {  {
1210    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1211    if (argDomain!=*this)     if (argDomain!=*this)
1212       throw FinleyAdapterException("Error - Illegal domain of data point locations");        throw FinleyAdapterException("Error - Illegal domain of data point locations");
1213    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1214    // in case of values node coordinates we can do the job directly:     // in case of values node coordinates we can do the job directly:
1215    if (arg.getFunctionSpace().getTypeCode()==Nodes) {     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1216       escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1217       Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1218    } else {     } else {
1219       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1220       escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1221       Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1222       // this is then interpolated onto arg:        // this is then interpolated onto arg:
1223       interpolateOnDomain(arg,tmp_data);        interpolateOnDomain(arg,tmp_data);
1224    }     }
1225    checkFinleyError();     checkFinleyError();
1226  }  }
1227    
1228  //  //
# Line 1174  void MeshAdapter::setToX(escript::Data& Line 1230  void MeshAdapter::setToX(escript::Data&
1230  //  //
1231  void MeshAdapter::setToNormal(escript::Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
1232  {  {
1233    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1234    if (normalDomain!=*this)     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1235       throw FinleyAdapterException("Error - Illegal domain of normal locations");     if (normalDomain!=*this)
1236    Finley_Mesh* mesh=m_finleyMesh.get();        throw FinleyAdapterException("Error - Illegal domain of normal locations");
1237    escriptDataC _normal=normal.getDataC();     Finley_Mesh* mesh=m_finleyMesh.get();
1238    switch(normal.getFunctionSpace().getTypeCode()) {     escriptDataC _normal=normal.getDataC();
1239      case(Nodes):     switch(normal.getFunctionSpace().getTypeCode()) {
1240        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");     case(Nodes):
1241        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1242      case(ReducedNodes):     break;
1243        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");     case(ReducedNodes):
1244        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1245      case(Elements):     break;
1246        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");     case(Elements):
1247        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1248      case(ReducedElements):     break;
1249        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");     case(ReducedElements):
1250        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1251      case (FaceElements):     break;
1252        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     case (FaceElements):
1253        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1254      case (ReducedFaceElements):     break;
1255        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     case (ReducedFaceElements):
1256        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1257      case(Points):     break;
1258        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");     case(Points):
1259        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1260      case (ContactElementsOne):     break;
1261      case (ContactElementsZero):     case (ContactElementsOne):
1262        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);     case (ContactElementsZero):
1263        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1264      case (ReducedContactElementsOne):     break;
1265      case (ReducedContactElementsZero):     case (ReducedContactElementsOne):
1266        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);     case (ReducedContactElementsZero):
1267        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1268      case(DegreesOfFreedom):     break;
1269        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");     case(DegreesOfFreedom):
1270        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1271      case(ReducedDegreesOfFreedom):     break;
1272        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");     case(ReducedDegreesOfFreedom):
1273        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1274      default:     break;
1275       default:
1276        stringstream temp;        stringstream temp;
1277        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();
1278        throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1279        break;        break;
1280    }     }
1281    checkFinleyError();     checkFinleyError();
1282  }  }
1283    
1284  //  //
# Line 1229  void MeshAdapter::setToNormal(escript::D Line 1286  void MeshAdapter::setToNormal(escript::D
1286  //  //
1287  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1288  {  {
1289    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1290    if (targetDomain!=*this)     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1291       throw FinleyAdapterException("Error - Illegal domain of interpolation target");     if (targetDomain!=this)
1292          throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1293    
1294    throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
1295  }  }
1296    
1297  //  //
1298  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
1299  //  //
1300  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1301  {  {
1302    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1303    if (argDomain!=*this)     if (argDomain!=*this)
1304       throw FinleyAdapterException("Error - Illegal domain of integration kernel");        throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1305    
1306    double blocktimer_start = blocktimer_time();     double blocktimer_start = blocktimer_time();
1307    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1308    escriptDataC _temp;     escriptDataC _temp;
1309    escript::Data temp;     escript::Data temp;
1310    escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1311    switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1312       case(Nodes):     case(Nodes):
1313          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1314          _temp=temp.getDataC();     _temp=temp.getDataC();
1315          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1316          break;     break;
1317       case(ReducedNodes):     case(ReducedNodes):
1318          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1319          _temp=temp.getDataC();     _temp=temp.getDataC();
1320          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1321          break;     break;
1322       case(Elements):     case(Elements):
1323          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1324          break;     break;
1325       case(ReducedElements):     case(ReducedElements):
1326          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1327          break;     break;
1328       case(FaceElements):     case(FaceElements):
1329          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1330          break;     break;
1331       case(ReducedFaceElements):     case(ReducedFaceElements):
1332          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1333          break;     break;
1334       case(Points):     case(Points):
1335          throw FinleyAdapterException("Error - Integral of data on points is not supported.");     throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1336          break;     break;
1337       case(ContactElementsZero):     case(ContactElementsZero):
1338          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1339          break;     break;
1340       case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1341          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1342          break;     break;
1343       case(ContactElementsOne):     case(ContactElementsOne):
1344          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1345          break;     break;
1346       case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1347          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1348          break;     break;
1349       case(DegreesOfFreedom):     case(DegreesOfFreedom):
1350          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1351          _temp=temp.getDataC();     _temp=temp.getDataC();
1352          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1353          break;     break;
1354       case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1355          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1356          _temp=temp.getDataC();     _temp=temp.getDataC();
1357          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1358          break;     break;
1359       default:     default:
1360          stringstream temp;        stringstream temp;
1361          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();
1362          throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1363          break;        break;
1364    }     }
1365    checkFinleyError();     checkFinleyError();
1366    blocktimer_increment("integrate()", blocktimer_start);     blocktimer_increment("integrate()", blocktimer_start);
1367  }  }
1368    
1369  //  //
# Line 1313  void MeshAdapter::setToIntegrals(std::ve Line 1371  void MeshAdapter::setToIntegrals(std::ve
1371  //  //
1372  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1373  {  {
1374    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1375    if (argDomain!=*this)     if (argDomain!=*this)
1376      throw FinleyAdapterException("Error - Illegal domain of gradient argument");        throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1377    const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1378    if (gradDomain!=*this)     if (gradDomain!=*this)
1379       throw FinleyAdapterException("Error - Illegal domain of gradient");        throw FinleyAdapterException("Error - Illegal domain of gradient");
1380    
1381    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1382    escriptDataC _grad=grad.getDataC();     escriptDataC _grad=grad.getDataC();
1383    escriptDataC nodeDataC;     escriptDataC nodeDataC;
1384    escript::Data temp;     escript::Data temp;
1385    if (getMPISize()>1) {     if (getMPISize()>1) {
1386        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1387          temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1388          nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1389        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1390          temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1391          nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1392        } else {        } else {
1393          nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
1394        }        }
1395    } else {     } else {
1396       nodeDataC = arg.getDataC();        nodeDataC = arg.getDataC();
1397    }     }
1398    switch(grad.getFunctionSpace().getTypeCode()) {     switch(grad.getFunctionSpace().getTypeCode()) {
1399         case(Nodes):     case(Nodes):
1400            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");     throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1401            break;     break;
1402         case(ReducedNodes):     case(ReducedNodes):
1403            throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1404            break;     break;
1405         case(Elements):     case(Elements):
1406            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1407            break;     break;
1408         case(ReducedElements):     case(ReducedElements):
1409            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1410            break;     break;
1411         case(FaceElements):     case(FaceElements):
1412            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1413            break;     break;
1414         case(ReducedFaceElements):     case(ReducedFaceElements):
1415            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1416            break;     break;
1417         case(Points):     case(Points):
1418            throw FinleyAdapterException("Error - Gradient at points is not supported.");     throw FinleyAdapterException("Error - Gradient at points is not supported.");
1419            break;     break;
1420         case(ContactElementsZero):     case(ContactElementsZero):
1421            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1422            break;     break;
1423         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1424            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1425            break;     break;
1426         case(ContactElementsOne):     case(ContactElementsOne):
1427            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1428            break;     break;
1429         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1430            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1431            break;     break;
1432         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1433            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1434            break;     break;
1435         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1436            throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1437            break;     break;
1438         default:     default:
1439            stringstream temp;        stringstream temp;
1440            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();
1441            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1442            break;        break;
1443    }     }
1444    checkFinleyError();     checkFinleyError();
1445  }  }
1446    
1447  //  //
# Line 1391  void MeshAdapter::setToGradient(escript: Line 1449  void MeshAdapter::setToGradient(escript:
1449  //  //
1450  void MeshAdapter::setToSize(escript::Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
1451  {  {
1452    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1453    escriptDataC tmp=size.getDataC();     escriptDataC tmp=size.getDataC();
1454    switch(size.getFunctionSpace().getTypeCode()) {     switch(size.getFunctionSpace().getTypeCode()) {
1455         case(Nodes):     case(Nodes):
1456            throw FinleyAdapterException("Error - Size of nodes is not supported.");     throw FinleyAdapterException("Error - Size of nodes is not supported.");
1457            break;     break;
1458         case(ReducedNodes):     case(ReducedNodes):
1459            throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1460            break;     break;
1461         case(Elements):     case(Elements):
1462            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1463            break;     break;
1464         case(ReducedElements):     case(ReducedElements):
1465            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1466            break;     break;
1467         case(FaceElements):     case(FaceElements):
1468            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1469            break;     break;
1470         case(ReducedFaceElements):     case(ReducedFaceElements):
1471            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1472            break;     break;
1473         case(Points):     case(Points):
1474            throw FinleyAdapterException("Error - Size of point elements is not supported.");     throw FinleyAdapterException("Error - Size of point elements is not supported.");
1475            break;     break;
1476         case(ContactElementsZero):     case(ContactElementsZero):
1477         case(ContactElementsOne):     case(ContactElementsOne):
1478            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1479            break;     break;
1480         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1481         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1482            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1483            break;     break;
1484         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1485            throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1486            break;     break;
1487         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1488            throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1489            break;     break;
1490         default:     default:
1491            stringstream temp;        stringstream temp;
1492            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();
1493            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1494            break;        break;
1495    }     }
1496    checkFinleyError();     checkFinleyError();
1497  }  }
1498    
1499  // sets the location of nodes:  //
1500    // sets the location of nodes
1501    //
1502  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1503  {  {
1504    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1505    escriptDataC tmp;     escriptDataC tmp;
1506    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1507    if (newDomain!=*this)     if (newDomain!=*this)
1508       throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1509    tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1510    Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1511    checkFinleyError();         Finley_Mesh_setCoordinates(mesh,&tmp);
1512  }     } else {
1513           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1514  // saves a data array in openDX format:         tmp = new_x_inter.getDataC();
1515  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const         Finley_Mesh_setCoordinates(mesh,&tmp);
1516  {     }
1517      int MAX_namelength=256;     checkFinleyError();
1518      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);  
1519    
1520      return;  //
1521    // Helper for the save* methods. Extracts optional data variable names and the
1522    // corresponding pointers from python dictionary. Caller must free arrays.
1523    //
1524    void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1525    {
1526       numData = boost::python::extract<int>(arg.attr("__len__")());
1527       /* win32 refactor */
1528       names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1529       data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1530       dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1531    
1532       boost::python::list keys=arg.keys();
1533       for (int i=0; i<numData; ++i) {
1534          string n=boost::python::extract<string>(keys[i]);
1535          escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1536          if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1537             throw FinleyAdapterException("Error: Data must be defined on same Domain");
1538          data[i] = d.getDataC();
1539          dataPtr[i] = &(data[i]);
1540          names[i] = TMPMEMALLOC(n.length()+1, char);
1541          strcpy(names[i], n.c_str());
1542       }
1543  }  }
1544    
1545  // saves a data array in openVTK format:  //
1546  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  // saves mesh and optionally data arrays in openDX format
1547    //
1548    void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1549  {  {
1550      int MAX_namelength=256;     int num_data;
1551      const int num_data=boost::python::extract<int>(arg.attr("__len__")());     char **names;
1552    /* win32 refactor */     escriptDataC *data;
1553    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;  
   }  
1554    
1555    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1556    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);
1557    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;     checkFinleyError();
1558    
1559      boost::python::list keys=arg.keys();     /* win32 refactor */
1560      for (int i=0;i<num_data;++i) {     TMPMEMFREE(data);
1561           std::string n=boost::python::extract<std::string>(keys[i]);     TMPMEMFREE(ptr_data);
1562           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);     for(int i=0; i<num_data; i++)
1563           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)     {
1564               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");        TMPMEMFREE(names[i]);
1565           data[i]=d.getDataC();     }
1566           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);  
1567    
1568    checkFinleyError();     return;
1569    /* win32 refactor */  }
   TMPMEMFREE(c_names);  
   TMPMEMFREE(data);  
   TMPMEMFREE(ptr_data);  
   for(int i=0;i<num_data;i++)  
   {  
     TMPMEMFREE(names[i]);  
   }  
   TMPMEMFREE(names);  
1570    
1571      return;  //
1572    // saves mesh and optionally data arrays in VTK format
1573    //
1574    void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1575    {
1576        boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1577        pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1578                  metadata, metadata_schema, arg);
1579  }  }
1580                                                                                                                                                                            
1581                                                                                                                                                                            bool MeshAdapter::ownSample(int fs_code, index_t id) const
1582  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  {
1583  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  #ifdef ESYS_MPI
1584                        const int row_blocksize,      index_t myFirstNode=0, myLastNode=0, k=0;
1585                        const escript::FunctionSpace& row_functionspace,      index_t* globalNodeIndex=0;
1586                        const int column_blocksize,      Finley_Mesh* mesh_p=m_finleyMesh.get();
1587                        const escript::FunctionSpace& column_functionspace,      if (fs_code == FINLEY_REDUCED_NODES)
1588                        const int type) const      {
1589  {      myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1590      int reduceRowOrder=0;      myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1591      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.");  
1592      }      }
1593      if (column_functionspace.getTypeCode()==DegreesOfFreedom) {      else
1594          reduceColOrder=0;      {
1595      } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1596          reduceColOrder=1;      myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1597      } else {      globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
         throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");  
1598      }      }
1599      // generate matrix:      k=globalNodeIndex[id];
1600            return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1601      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);  #endif
1602      checkFinleyError();      return true;
1603      Paso_SystemMatrix* fsystemMatrix;  }
1604      int trilinos = 0;  
1605      if (trilinos) {  
1606    
1607    //
1608    // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1609    //
1610    ASM_ptr MeshAdapter::newSystemMatrix(
1611                                                     const int row_blocksize,
1612                                                     const escript::FunctionSpace& row_functionspace,
1613                                                     const int column_blocksize,
1614                                                     const escript::FunctionSpace& column_functionspace,
1615                                                     const int type) const
1616    {
1617       int reduceRowOrder=0;
1618       int reduceColOrder=0;
1619       // is the domain right?
1620       const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1621       if (row_domain!=*this)
1622          throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1623       const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1624       if (col_domain!=*this)
1625          throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1626       // is the function space type right
1627       if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1628          reduceRowOrder=0;
1629       } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1630          reduceRowOrder=1;
1631       } else {
1632          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1633       }
1634       if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1635          reduceColOrder=0;
1636       } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1637          reduceColOrder=1;
1638       } else {
1639          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1640       }
1641       // generate matrix:
1642    
1643       Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1644       checkFinleyError();
1645       Paso_SystemMatrix* fsystemMatrix;
1646       int trilinos = 0;
1647       if (trilinos) {
1648  #ifdef TRILINOS  #ifdef TRILINOS
1649        /* Allocation Epetra_VrbMatrix here */        /* Allocation Epetra_VrbMatrix here */
1650  #endif  #endif
1651      }     }
1652      else {     else {
1653        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1654      }     }
1655      checkPasoError();     checkPasoError();
1656      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1657      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1658       return ASM_ptr(sma);
1659    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1660  }  }
1661    
1662    //
1663  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1664  TransportProblemAdapter MeshAdapter::newTransportProblem(  //
1665                        const double theta,  ATP_ptr MeshAdapter::newTransportProblem(
1666                        const double dt_max,                                                           const bool useBackwardEuler,
1667                        const int blocksize,                                                           const int blocksize,
1668                        const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1669                        const int type) const                                                           const int type) const
1670  {  {
1671      int reduceOrder=0;     int reduceOrder=0;
1672      // is the domain right?     // is the domain right?
1673      const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(functionspace.getDomain());     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1674      if (domain!=*this)     if (domain!=*this)
1675            throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1676      // is the function space type right     // is the function space type right
1677      if (functionspace.getTypeCode()==DegreesOfFreedom) {     if (functionspace.getTypeCode()==DegreesOfFreedom) {
1678          reduceOrder=0;        reduceOrder=0;
1679      } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1680          reduceOrder=1;        reduceOrder=1;
1681      } else {     } else {
1682          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1683      }     }
1684      // generate matrix:     // generate matrix:
1685        
1686      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1687      checkFinleyError();     checkFinleyError();
1688      Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1689      transportProblem=Paso_FCTransportProblem_alloc(theta,dt_max,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1690      checkPasoError();     checkPasoError();
1691      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1692      return TransportProblemAdapter(transportProblem,theta,dt_max,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1693       return ATP_ptr(tpa);
1694    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1695  }  }
1696    
1697  //  //
# Line 1639  TransportProblemAdapter MeshAdapter::new Line 1703  TransportProblemAdapter MeshAdapter::new
1703  // 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:
1704  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1705  {  {
1706    switch(functionSpaceCode) {     switch(functionSpaceCode) {
1707         case(Nodes):     case(Nodes):
1708         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1709         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1710            return false;     return false;
1711            break;     break;
1712         case(Elements):     case(Elements):
1713         case(FaceElements):     case(FaceElements):
1714         case(Points):     case(Points):
1715         case(ContactElementsZero):     case(ContactElementsZero):
1716         case(ContactElementsOne):     case(ContactElementsOne):
1717         case(ReducedElements):     case(ReducedElements):
1718         case(ReducedFaceElements):     case(ReducedFaceElements):
1719         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1720         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1721            return true;     return true;
1722            break;     break;
1723         default:     default:
1724            stringstream temp;        stringstream temp;
1725            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;
1726            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1727            break;        break;
1728    }     }
1729    checkFinleyError();     checkFinleyError();
1730    return false;     return false;
1731    }
1732    
1733    bool
1734    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1735    {
1736       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1737        class 1: DOF <-> Nodes
1738        class 2: ReducedDOF <-> ReducedNodes
1739        class 3: Points
1740        class 4: Elements
1741        class 5: ReducedElements
1742        class 6: FaceElements
1743        class 7: ReducedFaceElements
1744        class 8: ContactElementZero <-> ContactElementOne
1745        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1746    
1747       There is also a set of lines. Interpolation is possible down a line but not between lines.
1748       class 1 and 2 belong to all lines so aren't considered.
1749        line 0: class 3
1750        line 1: class 4,5
1751        line 2: class 6,7
1752        line 3: class 8,9
1753    
1754       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1755       eg hasnodes is true if we have at least one instance of Nodes.
1756       */
1757        if (fs.empty())
1758        {
1759            return false;
1760        }
1761        vector<int> hasclass(10);
1762        vector<int> hasline(4);
1763        bool hasnodes=false;
1764        bool hasrednodes=false;
1765        bool hascez=false;
1766        bool hasrcez=false;
1767        for (int i=0;i<fs.size();++i)
1768        {
1769        switch(fs[i])
1770        {
1771        case(Nodes):   hasnodes=true;   // no break is deliberate
1772        case(DegreesOfFreedom):
1773            hasclass[1]=1;
1774            break;
1775        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1776        case(ReducedDegreesOfFreedom):
1777            hasclass[2]=1;
1778            break;
1779        case(Points):
1780            hasline[0]=1;
1781            hasclass[3]=1;
1782            break;
1783        case(Elements):
1784            hasclass[4]=1;
1785            hasline[1]=1;
1786            break;
1787        case(ReducedElements):
1788            hasclass[5]=1;
1789            hasline[1]=1;
1790            break;
1791        case(FaceElements):
1792            hasclass[6]=1;
1793            hasline[2]=1;
1794            break;
1795        case(ReducedFaceElements):
1796            hasclass[7]=1;
1797            hasline[2]=1;
1798            break;
1799        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1800        case(ContactElementsOne):
1801            hasclass[8]=1;
1802            hasline[3]=1;
1803            break;
1804        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1805        case(ReducedContactElementsOne):
1806            hasclass[9]=1;
1807            hasline[3]=1;
1808            break;
1809        default:
1810            return false;
1811        }
1812        }
1813        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1814        // fail if we have more than one leaf group
1815    
1816        if (totlines>1)
1817        {
1818        return false;   // there are at least two branches we can't interpolate between
1819        }
1820        else if (totlines==1)
1821        {
1822        if (hasline[0]==1)      // we have points
1823        {
1824            resultcode=Points;
1825        }
1826        else if (hasline[1]==1)
1827        {
1828            if (hasclass[5]==1)
1829            {
1830            resultcode=ReducedElements;
1831            }
1832            else
1833            {
1834            resultcode=Elements;
1835            }
1836        }
1837        else if (hasline[2]==1)
1838        {
1839            if (hasclass[7]==1)
1840            {
1841            resultcode=ReducedFaceElements;
1842            }
1843            else
1844            {
1845            resultcode=FaceElements;
1846            }
1847        }
1848        else    // so we must be in line3
1849        {
1850            if (hasclass[9]==1)
1851            {
1852            // need something from class 9
1853            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1854            }
1855            else
1856            {
1857            // something from class 8
1858            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1859            }
1860        }
1861        }
1862        else    // totlines==0
1863        {
1864        if (hasclass[2]==1)
1865        {
1866            // something from class 2
1867            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1868        }
1869        else
1870        {   // something from class 1
1871            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1872        }
1873        }
1874        return true;
1875  }  }
1876    
1877  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1878  {  {
1879    switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1880       case(Nodes):     case(Nodes):
1881          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1882             case(Nodes):      case(Nodes):
1883             case(ReducedNodes):      case(ReducedNodes):
1884             case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1885             case(DegreesOfFreedom):      case(DegreesOfFreedom):
1886             case(Elements):      case(Elements):
1887             case(ReducedElements):      case(ReducedElements):
1888             case(FaceElements):      case(FaceElements):
1889             case(ReducedFaceElements):      case(ReducedFaceElements):
1890             case(Points):      case(Points):
1891             case(ContactElementsZero):      case(ContactElementsZero):
1892             case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1893             case(ContactElementsOne):      case(ContactElementsOne):
1894             case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1895                 return true;      return true;
1896             default:      default:
1897                 stringstream temp;            stringstream temp;
1898                 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;
1899                 throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1900          }     }
1901          break;     break;
1902       case(ReducedNodes):     case(ReducedNodes):
1903          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1904             case(ReducedNodes):      case(ReducedNodes):
1905             case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1906             case(Elements):      case(Elements):
1907             case(ReducedElements):      case(ReducedElements):
1908             case(FaceElements):      case(FaceElements):
1909             case(ReducedFaceElements):      case(ReducedFaceElements):
1910             case(Points):      case(Points):
1911             case(ContactElementsZero):      case(ContactElementsZero):
1912             case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1913             case(ContactElementsOne):      case(ContactElementsOne):
1914             case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1915                 return true;      return true;
1916            case(Nodes):      case(Nodes):
1917            case(DegreesOfFreedom):      case(DegreesOfFreedom):
1918               return false;      return false;
1919             default:      default:
1920                 stringstream temp;          stringstream temp;
1921                 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;
1922                 throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1923          }     }
1924          break;     break;
1925       case(Elements):     case(Elements):
1926          if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1927             return true;        return true;
1928          } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1929             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;  
1930          } else {          } else {
1931             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());  
1932          }          }
1933          break;     case(ReducedElements):
1934       case(ReducedDegreesOfFreedom):      if (functionSpaceType_target==ReducedElements) {
1935         switch(functionSpaceType_target) {        return true;
1936            case(ReducedDegreesOfFreedom):      } else {
1937            case(ReducedNodes):            return false;
1938            case(Elements):      }
1939            case(ReducedElements):     case(FaceElements):
1940            case(FaceElements):      if (functionSpaceType_target==FaceElements) {
1941            case(ReducedFaceElements):              return true;
1942            case(Points):      } else if (functionSpaceType_target==ReducedFaceElements) {
1943            case(ContactElementsZero):              return true;
1944            case(ReducedContactElementsZero):      } else {
1945            case(ContactElementsOne):              return false;
1946            case(ReducedContactElementsOne):      }
1947                return true;     case(ReducedFaceElements):
1948            case(Nodes):      if (functionSpaceType_target==ReducedFaceElements) {
1949            case(DegreesOfFreedom):              return true;
1950               return false;      } else {
1951            default:          return false;
1952               stringstream temp;      }
1953               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;     case(Points):
1954               throw FinleyAdapterException(temp.str());      if (functionSpaceType_target==Points) {
1955         }              return true;
1956         break;      } else {
1957       default:              return false;
1958          stringstream temp;      }
1959          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;     case(ContactElementsZero):
1960          throw FinleyAdapterException(temp.str());     case(ContactElementsOne):
1961          break;      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1962    }              return true;
1963    checkFinleyError();      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1964    return false;              return true;
1965        } else {
1966                return false;
1967        }
1968       case(ReducedContactElementsZero):
1969       case(ReducedContactElementsOne):
1970        if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1971                return true;
1972        } else {
1973                return false;
1974        }
1975       case(DegreesOfFreedom):
1976        switch(functionSpaceType_target) {
1977        case(ReducedDegreesOfFreedom):
1978        case(DegreesOfFreedom):
1979        case(Nodes):
1980        case(ReducedNodes):
1981        case(Elements):
1982        case(ReducedElements):
1983        case(Points):
1984        case(FaceElements):
1985        case(ReducedFaceElements):
1986        case(ContactElementsZero):
1987        case(ReducedContactElementsZero):
1988        case(ContactElementsOne):
1989        case(ReducedContactElementsOne):
1990        return true;
1991        default:
1992            stringstream temp;
1993            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1994            throw FinleyAdapterException(temp.str());
1995        }
1996        break;
1997       case(ReducedDegreesOfFreedom):
1998       switch(functionSpaceType_target) {
1999        case(ReducedDegreesOfFreedom):
2000        case(ReducedNodes):
2001        case(Elements):
2002        case(ReducedElements):
2003        case(FaceElements):
2004        case(ReducedFaceElements):
2005        case(Points):
2006        case(ContactElementsZero):
2007        case(ReducedContactElementsZero):
2008        case(ContactElementsOne):
2009        case(ReducedContactElementsOne):
2010        return true;
2011        case(Nodes):
2012        case(DegreesOfFreedom):
2013        return false;
2014        default:
2015            stringstream temp;
2016            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
2017            throw FinleyAdapterException(temp.str());
2018        }
2019        break;
2020       default:
2021          stringstream temp;
2022          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
2023          throw FinleyAdapterException(temp.str());
2024          break;
2025       }
2026       checkFinleyError();
2027       return false;
2028  }  }
2029    
2030  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
2031  {  {
2032      return false;     return false;
2033  }  }
2034    
2035  bool MeshAdapter::operator==(const AbstractDomain& other) const  bool MeshAdapter::operator==(const AbstractDomain& other) const
2036  {  {
2037    const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);     const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2038    if (temp!=0) {     if (temp!=0) {
2039      return (m_finleyMesh==temp->m_finleyMesh);        return (m_finleyMesh==temp->m_finleyMesh);
2040    } else {     } else {
2041      return false;        return false;
2042    }     }
2043  }  }
2044    
2045  bool MeshAdapter::operator!=(const AbstractDomain& other) const  bool MeshAdapter::operator!=(const AbstractDomain& other) const
2046  {  {
2047    return !(operator==(other));     return !(operator==(other));
2048  }  }
2049    
2050  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
2051  {  {
2052     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2053       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2054       checkPasoError();
2055       return out;
2056    }
2057    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2058    {
2059       Finley_Mesh* mesh=m_finleyMesh.get();
2060       int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2061     checkPasoError();     checkPasoError();
2062     return out;     return out;
2063  }  }
2064    
2065  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2066  {  {
2067    return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2068  }  }
2069    
2070  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2071  {  {
2072    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2073  }  }
2074    
2075  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2076  {  {
2077    return function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2078  }  }
2079    
2080  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2081  {  {
2082    int *out=0,i;     int *out = NULL;
2083    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2084    switch (functionSpaceType) {     switch (functionSpaceType) {
2085      case(Nodes):     case(Nodes):
2086        out=mesh->Nodes->Id;     out=mesh->Nodes->Id;
2087        break;     break;
2088      case(ReducedNodes):     case(ReducedNodes):
2089        out=mesh->Nodes->reducedNodesId;     out=mesh->Nodes->reducedNodesId;
2090        break;     break;
2091      case(Elements):     case(Elements):
2092        out=mesh->Elements->Id;     out=mesh->Elements->Id;
2093        break;     break;
2094      case(ReducedElements):     case(ReducedElements):
2095        out=mesh->Elements->Id;     out=mesh->Elements->Id;
2096        break;     break;
2097      case(FaceElements):     case(FaceElements):
2098        out=mesh->FaceElements->Id;     out=mesh->FaceElements->Id;
2099        break;     break;
2100      case(ReducedFaceElements):     case(ReducedFaceElements):
2101        out=mesh->FaceElements->Id;     out=mesh->FaceElements->Id;
2102        break;     break;
2103      case(Points):     case(Points):
2104        out=mesh->Points->Id;     out=mesh->Points->Id;
2105        break;     break;
2106      case(ContactElementsZero):     case(ContactElementsZero):
2107        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2108        break;     break;
2109      case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
2110        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2111        break;     break;
2112      case(ContactElementsOne):     case(ContactElementsOne):
2113        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2114        break;     break;
2115      case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
2116        out=mesh->ContactElements->Id;     out=mesh->ContactElements->Id;
2117        break;     break;
2118      case(DegreesOfFreedom):     case(DegreesOfFreedom):
2119        out=mesh->Nodes->degreesOfFreedomId;     out=mesh->Nodes->degreesOfFreedomId;
2120        break;     break;
2121      case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2122        out=mesh->Nodes->reducedDegreesOfFreedomId;     out=mesh->Nodes->reducedDegreesOfFreedomId;
2123        break;     break;
2124      default:     default:
2125        stringstream temp;        stringstream temp;
2126        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2127        throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2128        break;        break;
2129    }     }
2130    return out;     return out;
2131  }  }
2132  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2133  {  {
2134    int out=0;     int out=0;
2135    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2136    switch (functionSpaceType) {     switch (functionSpaceType) {
2137    case(Nodes):     case(Nodes):
2138      out=mesh->Nodes->Tag[sampleNo];     out=mesh->Nodes->Tag[sampleNo];
2139      break;     break;
2140    case(ReducedNodes):     case(ReducedNodes):
2141      throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2142      break;     break;
2143    case(Elements):     case(Elements):
2144      out=mesh->Elements->Tag[sampleNo];     out=mesh->Elements->Tag[sampleNo];
2145      break;     break;
2146    case(ReducedElements):     case(ReducedElements):
2147      out=mesh->Elements->Tag[sampleNo];     out=mesh->Elements->Tag[sampleNo];
2148      break;     break;
2149    case(FaceElements):     case(FaceElements):
2150      out=mesh->FaceElements->Tag[sampleNo];     out=mesh->FaceElements->Tag[sampleNo];
2151      break;     break;
2152    case(ReducedFaceElements):     case(ReducedFaceElements):
2153      out=mesh->FaceElements->Tag[sampleNo];     out=mesh->FaceElements->Tag[sampleNo];
2154      break;     break;
2155    case(Points):     case(Points):
2156      out=mesh->Points->Tag[sampleNo];     out=mesh->Points->Tag[sampleNo];
2157      break;     break;
2158    case(ContactElementsZero):     case(ContactElementsZero):
2159      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2160      break;     break;
2161    case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
2162      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2163      break;     break;
2164    case(ContactElementsOne):     case(ContactElementsOne):
2165      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2166      break;     break;
2167    case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
2168      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2169      break;     break;
2170    case(DegreesOfFreedom):     case(DegreesOfFreedom):
2171      throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2172      break;     break;
2173    case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2174      throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");     throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2175      break;     break;
2176    default:     default:
2177      stringstream temp;        stringstream temp;
2178      temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2179      throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2180      break;        break;
2181    }     }
2182    return out;     return out;
2183  }  }
2184    
2185    
2186  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
2187  {  {
2188       Finley_Mesh* mesh=m_finleyMesh.get();
2189       escriptDataC tmp=mask.getDataC();
2190       switch(functionSpaceType) {
2191       case(Nodes):
2192       Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2193       break;
2194       case(ReducedNodes):
2195       throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2196       break;
2197       case(DegreesOfFreedom):
2198       throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2199       break;
2200       case(ReducedDegreesOfFreedom):
2201       throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2202       break;
2203       case(Elements):
2204       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2205       break;
2206       case(ReducedElements):
2207       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2208       break;
2209       case(FaceElements):
2210       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2211       break;
2212       case(ReducedFaceElements):
2213       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2214       break;
2215       case(Points):
2216       Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2217       break;
2218       case(ContactElementsZero):
2219       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2220       break;
2221       case(ReducedContactElementsZero):
2222       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2223       break;
2224       case(ContactElementsOne):
2225       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2226       break;
2227       case(ReducedContactElementsOne):
2228       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2229       break;
2230       default:
2231          stringstream temp;
2232          temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2233          throw FinleyAdapterException(temp.str());
2234       }
2235       checkFinleyError();
2236       return;
2237    }
2238    
2239    void MeshAdapter::setTagMap(const string& name,  int tag)
2240    {
2241       Finley_Mesh* mesh=m_finleyMesh.get();
2242       Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2243       checkPasoError();
2244       // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2245    }
2246    
2247    int MeshAdapter::getTag(const string& name) const
2248    {
2249       Finley_Mesh* mesh=m_finleyMesh.get();
2250       int tag=0;
2251       tag=Finley_Mesh_getTag(mesh, name.c_str());
2252       checkPasoError();
2253       // throwStandardException("MeshAdapter::getTag is not implemented.");
2254       return tag;
2255    }
2256    
2257    bool MeshAdapter::isValidTagName(const string& name) const
2258    {
2259       Finley_Mesh* mesh=m_finleyMesh.get();
2260       return Finley_Mesh_isValidTagName(mesh,name.c_str());
2261    }
2262    
2263    string MeshAdapter::showTagNames() const
2264    {
2265       stringstream temp;
2266       Finley_Mesh* mesh=m_finleyMesh.get();
2267       Finley_TagMap* tag_map=mesh->TagMap;
2268       while (tag_map) {
2269          temp << tag_map->name;
2270          tag_map=tag_map->next;
2271          if (tag_map) temp << ", ";
2272       }
2273       return temp.str();
2274    }
2275    
2276    int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2277    {
2278    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2279    escriptDataC tmp=mask.getDataC();    dim_t numTags=0;
2280    switch(functionSpaceType) {    switch(functionSpaceCode) {
2281         case(Nodes):     case(Nodes):
2282            Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);            numTags=mesh->Nodes->numTagsInUse;
2283            break;            break;
2284         case(ReducedNodes):     case(ReducedNodes):
2285            throw FinleyAdapterException("Error - ReducedNodes does not support tags");            throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2286            break;            break;
2287         case(DegreesOfFreedom):     case(DegreesOfFreedom):
2288            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2289            break;            break;
2290         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2291            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2292            break;            break;
2293         case(Elements):     case(Elements):
2294            Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     case(ReducedElements):
2295            break;            numTags=mesh->Elements->numTagsInUse;
2296         case(ReducedElements):            break;
2297            Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     case(FaceElements):
2298       case(ReducedFaceElements):
2299              numTags=mesh->FaceElements->numTagsInUse;
2300              break;
2301       case(Points):
2302              numTags=mesh->Points->numTagsInUse;
2303              break;
2304       case(ContactElementsZero):
2305       case(ReducedContactElementsZero):
2306       case(ContactElementsOne):
2307       case(ReducedContactElementsOne):
2308              numTags=mesh->ContactElements->numTagsInUse;
2309            break;            break;
2310         case(FaceElements):     default:
2311            Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);        stringstream temp;
2312            break;        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2313         case(ReducedFaceElements):        throw FinleyAdapterException(temp.str());
2314            Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);    }
2315            break;    return numTags;
2316         case(Points):  }
2317            Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);  
2318    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2319    {
2320      Finley_Mesh* mesh=m_finleyMesh.get();
2321      index_t* tags=NULL;
2322      switch(functionSpaceCode) {
2323       case(Nodes):
2324              tags=mesh->Nodes->tagsInUse;
2325            break;            break;
2326         case(ContactElementsZero):     case(ReducedNodes):
2327            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2328            break;            break;
2329         case(ReducedContactElementsZero):     case(DegreesOfFreedom):
2330            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2331            break;            break;
2332         case(ContactElementsOne):     case(ReducedDegreesOfFreedom):
2333            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2334            break;            break;
2335         case(ReducedContactElementsOne):     case(Elements):
2336            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);     case(ReducedElements):
2337              tags=mesh->Elements->tagsInUse;
2338              break;
2339       case(FaceElements):
2340       case(ReducedFaceElements):
2341              tags=mesh->FaceElements->tagsInUse;
2342              break;
2343       case(Points):
2344              tags=mesh->Points->tagsInUse;
2345              break;
2346       case(ContactElementsZero):
2347       case(ReducedContactElementsZero):
2348       case(ContactElementsOne):
2349       case(ReducedContactElementsOne):
2350              tags=mesh->ContactElements->tagsInUse;
2351            break;            break;
2352         default:     default:
2353            stringstream temp;        stringstream temp;
2354            temp << "Error - Finley does not know anything about function space type " << functionSpaceType;        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2355            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2356    }    }
2357    checkFinleyError();    return tags;
   return;  
2358  }  }
2359    
2360  void MeshAdapter::setTagMap(const std::string& name,  int tag)  
2361    bool MeshAdapter::canTag(int functionSpaceCode) const
2362  {  {
2363    Finley_Mesh* mesh=m_finleyMesh.get();    switch(functionSpaceCode) {
2364    Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     case(Nodes):
2365    checkPasoError();     case(Elements):
2366    // throwStandardException("MeshAdapter::set TagMap is not implemented.");     case(ReducedElements):
2367       case(FaceElements):
2368       case(ReducedFaceElements):
2369       case(Points):
2370       case(ContactElementsZero):
2371       case(ReducedContactElementsZero):
2372       case(ContactElementsOne):
2373       case(ReducedContactElementsOne):
2374              return true;
2375       case(ReducedNodes):
2376       case(DegreesOfFreedom):
2377       case(ReducedDegreesOfFreedom):
2378          return false;
2379       default:
2380        return false;
2381      }
2382  }  }
2383    
2384  int MeshAdapter::getTag(const std::string& name) const  AbstractDomain::StatusType MeshAdapter::getStatus() const
2385  {  {
2386    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2387    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;  
2388  }  }
2389    
2390  bool MeshAdapter::isValidTagName(const std::string& name) const  int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2391  {  {
2392      
2393    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2394    return Finley_Mesh_isValidTagName(mesh,name.c_str());    int order =-1;
2395      switch(functionSpaceCode) {
2396       case(Nodes):
2397       case(DegreesOfFreedom):
2398              order=mesh->approximationOrder;
2399              break;
2400       case(ReducedNodes):
2401       case(ReducedDegreesOfFreedom):
2402              order=mesh->reducedApproximationOrder;
2403              break;
2404       case(Elements):
2405       case(FaceElements):
2406       case(Points):
2407       case(ContactElementsZero):
2408       case(ContactElementsOne):
2409              order=mesh->integrationOrder;
2410              break;
2411       case(ReducedElements):
2412       case(ReducedFaceElements):
2413       case(ReducedContactElementsZero):
2414       case(ReducedContactElementsOne):
2415              order=mesh->reducedIntegrationOrder;
2416              break;
2417       default:
2418          stringstream temp;
2419          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2420          throw FinleyAdapterException(temp.str());
2421      }
2422      return order;
2423  }  }
2424    
2425  std::string MeshAdapter::showTagNames() const  bool MeshAdapter::supportsContactElements() const
2426  {  {
2427    stringstream temp;    return true;
2428    Finley_Mesh* mesh=m_finleyMesh.get();  }
2429    Finley_TagMap* tag_map=mesh->TagMap;  
2430    while (tag_map) {  ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2431       temp << tag_map->name;  {
2432       tag_map=tag_map->next;    m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2433       if (tag_map) temp << ", ";  }
2434    }  
2435    return temp.str();  ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2436    {
2437      Finley_ReferenceElementSet_dealloc(m_refSet);
2438  }  }
2439    
2440    void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2441    {
2442          const int dim = getDim();
2443          int numPoints=boost::python::extract<int>(points.attr("__len__")());
2444          int numTags=boost::python::extract<int>(tags.attr("__len__")());
2445          Finley_Mesh* mesh=m_finleyMesh.get();
2446          
2447          if  ( (numTags > 0) and ( numPoints !=  numTags ) )
2448         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2449          
2450          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2451          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2452          
2453          for (int i=0;i<numPoints;++i) {
2454           int tag_id=-1;
2455           int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2456           if  ( numComps !=   dim ) {
2457                   stringstream temp;          
2458                   temp << "Error - illegal number of components " << numComps << " for point " << i;              
2459                   throw FinleyAdapterException(temp.str());
2460           }
2461           points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2462           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2463           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2464          
2465           if ( numTags > 0) {
2466                  boost::python::extract<string> ex_str(tags[i]);
2467              if  ( ex_str.check() ) {
2468                  tag_id=getTag( ex_str());
2469              } else {
2470                   boost::python::extract<int> ex_int(tags[i]);
2471                   if ( ex_int.check() ) {
2472                       tag_id=ex_int();
2473                   } else {
2474                    stringstream temp;          
2475                        temp << "Error - unable to extract tag for point " << i;
2476                    throw FinleyAdapterException(temp.str());
2477                  }
2478              }
2479           }      
2480               tags_ptr[i]=tag_id;
2481          }
2482          
2483          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2484          checkPasoError();
2485          
2486          TMPMEMFREE(points_ptr);
2487          TMPMEMFREE(tags_ptr);
2488    }
2489    
2490    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2491    {  
2492        boost::python::list points =  boost::python::list();
2493        boost::python::list tags =  boost::python::list();
2494        points.append(point);
2495        tags.append(tag);
2496        addDiracPoints(points, tags);
2497    }
2498    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2499    {
2500            boost::python::list points =   boost::python::list();
2501            boost::python::list tags =   boost::python::list();
2502            points.append(point);
2503            tags.append(tag);
2504            addDiracPoints(points, tags);
2505    }  
2506  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1375  
changed lines
  Added in v.3522

  ViewVC Help
Powered by ViewVC 1.1.26