/[escript]/trunk/ripley/src/RipleyDomain.cpp
ViewVC logotype

Diff of /trunk/ripley/src/RipleyDomain.cpp

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

revision 3670 by caltinay, Wed Nov 16 02:21:18 2011 UTC revision 3756 by caltinay, Fri Jan 6 02:35:19 2012 UTC
# Line 12  Line 12 
12  *******************************************************/  *******************************************************/
13    
14  #include <ripley/RipleyDomain.h>  #include <ripley/RipleyDomain.h>
 #include <ripley/IndexList.h>  
 #include <ripley/Util.h>  
 extern "C" {  
 #include "esysUtils/blocktimer.h"  
 }  
 #include <escript/Data.h>  
15  #include <escript/DataFactory.h>  #include <escript/DataFactory.h>
16    #include <escript/FunctionSpaceFactory.h>
17    #include <pasowrap/SystemMatrixAdapter.h>
18    #include <pasowrap/TransportProblemAdapter.h>
19    
 #ifdef USE_NETCDF  
 #include <netcdfcpp.h>  
 #endif  
   
 #ifdef USE_PARMETIS  
 #include <parmetis.h>  
 #endif  
   
 #include <boost/python/import.hpp>  
 #include <boost/python/tuple.hpp>  
20  #include <iomanip>  #include <iomanip>
21    
22  using namespace std;  using namespace std;
23    using paso::SystemMatrixAdapter;
24    using paso::TransportProblemAdapter;
25    
26  namespace ripley {  namespace ripley {
27    
28  //  escript::Domain_ptr RipleyDomain::loadMesh(const string& filename)
 // define the static constants  
 RipleyDomain::FunctionSpaceNamesMapType RipleyDomain::m_functionSpaceTypeNames;  
   
 RipleyDomain::RipleyDomain(const string& name, dim_t numDim, Esys_MPIInfo* mpiInfo) :  
     m_name(name)  
 {  
     setFunctionSpaceTypeNames();  
     m_mpiInfo = Esys_MPIInfo_getReference(mpiInfo);  
     m_nodes.reset(new NodeFile(numDim, m_mpiInfo));  
     if (numDim==3) {  
         m_elements.reset(new ElementFile(Hex8, m_mpiInfo));  
         m_faceElements.reset(new ElementFile(Rec4, m_mpiInfo));  
     } else if (numDim==2) {  
         m_elements.reset(new ElementFile(Rec4, m_mpiInfo));  
         m_faceElements.reset(new ElementFile(Line2, m_mpiInfo));  
     } else  
         throw RipleyException("RipleyDomain: Invalid number of dimensions");  
     m_points.reset(new ElementFile(Point1, m_mpiInfo));  
     m_fullFullPattern = NULL;  
     m_fullReducedPattern = NULL;  
     m_reducedFullPattern = NULL;  
     m_reducedReducedPattern = NULL;  
 }  
   
 RipleyDomain::~RipleyDomain()  
 {  
     Esys_MPIInfo_free(m_mpiInfo);  
 }  
   
 int RipleyDomain::getMPISize() const  
29  {  {
30      return m_mpiInfo->size;      throw RipleyException("loadMesh() not implemented");
31  }  }
32    
33  int RipleyDomain::getMPIRank() const  escript::Domain_ptr RipleyDomain::readMesh(const string& filename)
34  {  {
35      return m_mpiInfo->rank;      throw RipleyException("readMesh() not implemented");
 }  
   
 void RipleyDomain::MPIBarrier() const  
 {  
 #ifdef ESYS_MPI  
     MPI_Barrier(m_mpiInfo->comm);  
 #endif  
36  }  }
37    
38  bool RipleyDomain::onMasterProcessor() const  RipleyDomain::RipleyDomain(dim_t dim) :
39        m_numDim(dim),
40        m_status(0)
41  {  {
42      return m_mpiInfo->rank == 0;      m_mpiInfo = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
43  }  }
44    
45  #ifdef ESYS_MPI  RipleyDomain::~RipleyDomain()
 MPI_Comm  
 #else  
 unsigned int  
 #endif  
 RipleyDomain::getMPIComm() const  
 {  
 #ifdef ESYS_MPI  
     return m_mpiInfo->comm;  
 #else  
     return 0;  
 #endif  
 }  
   
 void RipleyDomain::write(const string& filename) const  
 {  
     if (m_mpiInfo->size > 1)  
         throw RipleyException("write: only single processor runs are supported.");  
   
     // open file  
     ofstream f(filename.c_str());  
     if (!f.is_open()) {  
         stringstream msg;  
         msg << "write: Opening file " << filename << " for writing failed.";  
         throw RipleyException(msg.str());  
     }  
   
     // write header  
     f << m_name << endl;  
   
     f.precision(15);  
     f.setf(ios::scientific, ios::floatfield);  
   
     // write nodes  
     dim_t numDim = getDim();  
     f << numDim << "D-Nodes " << m_nodes->getNumNodes() << endl;  
     for (dim_t i = 0; i < m_nodes->getNumNodes(); i++) {  
         f << m_nodes->getIdVector()[i] << " "  
             << m_nodes->getGlobalDegreesOfFreedom()[i] << " "  
             << m_nodes->getTagVector()[i];  
         for (dim_t j = 0; j < numDim; j++)  
             f << " " << setw(20) << m_nodes->getCoordinates()[INDEX2(j, i, numDim)];  
         f << endl;  
     }  
   
     // write all element types  
     const ElementFile_ptr eFiles[3] = { m_elements, m_faceElements, m_points };  
   
     for (size_t k=0; k<3; k++) {  
         f << eFiles[k]->getName() << " " << eFiles[k]->getNumElements() << endl;  
         dim_t NN = eFiles[k]->getNumNodes();  
         for (dim_t i = 0; i < eFiles[k]->getNumElements(); i++) {  
             f << eFiles[k]->getIdVector()[i] << " "  
                 << eFiles[k]->getTagVector()[i];  
             for (dim_t j = 0; j < NN; j++)  
                 f << " " << m_nodes->getIdVector()[eFiles[k]->getNodes()[INDEX2(j, i, NN)]];  
             f << endl;  
         }  
     }  
   
     // write tags  
     if (m_tagMap.size() > 0) {  
         f << "Tags" << endl;  
         TagMap::const_iterator it;  
         for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {  
             f << it->first << " " << it->second << endl;  
         }  
     }  
     f.close();  
 #ifdef Ripley_TRACE  
     cout << "mesh " << m_name << " has been written to file " << filename << endl;  
 #endif  
 }  
   
 void RipleyDomain::Print_Mesh_Info(const bool full) const  
46  {  {
47      cout << "Ripley_PrintMesh_Info running on CPU " <<      Esys_MPIInfo_free(m_mpiInfo);
         m_mpiInfo->rank << " of " << m_mpiInfo->size << endl;  
     cout << "\tMesh name '" << m_name << "'" << endl;  
   
     // write nodes  
     int numDim = getDim();  
     cout << "\tNodes: " << numDim << "D-Nodes " << m_nodes->getNumNodes() << endl;  
     if (full) {  
         cout << "\t     Id   Tag  gDOF   gNI grDfI  grNI:  Coordinates" << endl;  
         cout.precision(15);  
         cout.setf(ios::scientific, ios::floatfield);  
         for (int i = 0; i < m_nodes->getNumNodes(); i++) {  
             cout << "\t  " << setw(5) << m_nodes->getIdVector()[i] << " "  
                << setw(5) << m_nodes->getTagVector()[i] << " "  
                << setw(5) << m_nodes->getGlobalDegreesOfFreedom()[i] << " "  
                << setw(5) << m_nodes->getGlobalNodesIndex()[i] << " "  
                << setw(5) << m_nodes->getGlobalReducedDOFIndex()[i] << " "  
                << setw(5) << m_nodes->getGlobalReducedNodesIndex()[i] << ": ";  
             for (int j = 0; j < numDim; j++)  
                 cout << " " << m_nodes->getCoordinates()[INDEX2(j, i, numDim)];  
             cout << endl;  
         }  
     }  
   
     // write all element types  
     const char *eNames[3] = { "Elements", "Face elements", "Points" };  
     const ElementFile_ptr eFiles[3] = { m_elements, m_faceElements, m_points };  
   
     for (size_t k=0; k<3; k++) {  
         int mine = 0, overlap = 0;  
         for (dim_t i = 0; i < eFiles[k]->getNumElements(); i++) {  
             if (eFiles[k]->getOwnerVector()[i] == m_mpiInfo->rank)  
                 mine++;  
             else  
                 overlap++;  
         }  
         cout << "\t" << eNames[k] << ": " << eFiles[k]->getName()  
             << " " << eFiles[k]->getNumElements()  
             << " (TypeId=" << eFiles[k]->getTypeId() << ") owner="  
             << mine << " overlap=" << overlap << endl;  
         int NN = eFiles[k]->getNumNodes();  
         if (full) {  
             cout << "\t     Id   Tag Owner Color:  Nodes" << endl;  
             for (int i = 0; i < eFiles[k]->getNumElements(); i++) {  
                 cout << "\t  " << setw(5) << eFiles[k]->getIdVector()[i] << " "  
                    << setw(5) << eFiles[k]->getTagVector()[i] << " "  
                    << setw(5) << eFiles[k]->getOwnerVector()[i] << " "  
                    << setw(5) << eFiles[k]->getColorVector()[i] << ": ";  
                 for (int j = 0; j < NN; j++)  
                     cout << " " << setw(5) << m_nodes->getIdVector()[eFiles[k]->getNodes()[INDEX2(j, i, NN)]];  
                 cout << endl;  
             }  
         }  
     }  
   
   
     // write tags  
     if (m_tagMap.size() > 0) {  
         cout << "\tTags:" << endl;  
         TagMap::const_iterator it;  
         for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {  
             cout << "\t  " << setw(5) << it->second << " "  
                 << it->first << endl;  
         }  
     }  
48  }  }
49    
50  void RipleyDomain::dump(const string& fileName) const  bool RipleyDomain::operator==(const AbstractDomain& other) const
51  {  {
52  #ifdef USE_NETCDF      const RipleyDomain* o=dynamic_cast<const RipleyDomain*>(&other);
53      NcDim* ncdim = NULL;      if (o) {
54      int mpi_size = m_mpiInfo->size;          return (m_tagMap==o->m_tagMap && m_nodeTags==o->m_nodeTags
55      int mpi_rank = m_mpiInfo->rank;                  && m_elementTags==o->m_elementTags
56      int numDim = getDim();                  && m_faceTags==o->m_faceTags);
   
 /* Incoming token indicates it's my turn to write */  
 #ifdef ESYS_MPI  
     MPI_Status status;  
     if (mpi_rank>0) {  
         int dummy;  
         MPI_Recv(&dummy, 0, MPI_INT, mpi_rank-1, 81800, m_mpiInfo->comm, &status);  
     }  
 #endif  
   
     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),  
                                                       mpi_size, mpi_rank);  
   
     // NetCDF error handler  
     NcError err(NcError::verbose_nonfatal);  
     // Create the file.  
     NcFile dataFile(newFileName, NcFile::Replace);  
     const string msgPrefix("dump: netCDF operation failed - ");  
   
     // check if writing was successful  
     if (!dataFile.is_valid())  
         throw RipleyException(msgPrefix+"Open file for output");  
   
     const size_t numTags = m_tagMap.size();  
   
     if (numTags > 0)  
         if (! (ncdim = dataFile.add_dim("dim_Tags", numTags)) )  
             throw RipleyException(msgPrefix+"add_dim(dim_Tags)");  
   
     // Attributes: MPI size, MPI rank, Name, order, reduced_order  
     if (!dataFile.add_att("mpi_size", mpi_size))  
         throw RipleyException(msgPrefix+"add_att(mpi_size)");  
     if (!dataFile.add_att("mpi_rank", mpi_rank))  
         throw RipleyException(msgPrefix+"add_att(mpi_rank)");  
     if (!dataFile.add_att("Name", m_name.c_str()))  
         throw RipleyException(msgPrefix+"add_att(Name)");  
     if (!dataFile.add_att("numDim", numDim))  
         throw RipleyException(msgPrefix+"add_att(numDim)");  
     if (!dataFile.add_att("order", 1))  
         throw RipleyException(msgPrefix+"add_att(order)");  
     if (!dataFile.add_att("reduced_order", 1))  
         throw RipleyException(msgPrefix+"add_att(reduced_order)");  
     if (!dataFile.add_att("num_Tags", static_cast<int>(numTags)))  
         throw RipleyException(msgPrefix+"add_att(num_Tags)");  
   
     m_nodes->dumpToNetCDF(dataFile);  
     m_elements->dumpToNetCDF(dataFile, "Elements");  
     m_faceElements->dumpToNetCDF(dataFile, "FaceElements");  
     m_points->dumpToNetCDF(dataFile, "Points");  
   
     // // // // // TagMap // // // // //  
     if (numTags > 0) {  
   
         // Copy tag keys into temp array  
         vector<int> Tags_keys;  
         TagMap::const_iterator it;  
         for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {  
             Tags_keys.push_back(it->second);  
         }  
   
         // Tags_keys  
         NcVar *ncVar;  
         if (! (ncVar = dataFile.add_var("Tags_keys", ncInt, ncdim)) )  
             throw RipleyException(msgPrefix+"add_var(Tags_keys)");  
         if (!ncVar->put(&Tags_keys[0], numTags))  
             throw RipleyException(msgPrefix+"put(Tags_keys)");  
   
         // 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  
         int i = 0;  
         for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {  
             stringstream name;  
             name << "Tags_name_" << i;  
             if (!dataFile.add_att(name.str().c_str(), it->first.c_str()) )  
                 throw RipleyException(msgPrefix+"add_att(Tags_names_XX)");  
             i++;  
         }  
57      }      }
58        return false;
     // Send token to next MPI process so he can take his turn  
 #ifdef ESYS_MPI  
     if (mpi_rank<mpi_size-1) {  
         int dummy = 0;  
         MPI_Send(&dummy, 0, MPI_INT, mpi_rank+1, 81800, m_mpiInfo->comm);  
     }  
 #endif  
   
     // netCDF file is closed by destructor of NcFile object  
 #else  
     throw RipleyException("dump: not configured with netCDF. Please contact your installation manager.");  
 #endif  /* USE_NETCDF */  
 }  
   
 string RipleyDomain::getDescription() const  
 {  
     return "RipleyMesh";  
 }  
   
 string RipleyDomain::functionSpaceTypeAsString(int functionSpaceType) const  
 {  
     FunctionSpaceNamesMapType::iterator loc;  
     loc=m_functionSpaceTypeNames.find(functionSpaceType);  
     if (loc!=m_functionSpaceTypeNames.end())  
         return loc->second;  
   
     return "Invalid function space type code";  
 }  
   
 bool RipleyDomain::isValidFunctionSpaceType(int functionSpaceType) const  
 {  
     FunctionSpaceNamesMapType::iterator loc;  
     loc=m_functionSpaceTypeNames.find(functionSpaceType);  
     return (loc!=m_functionSpaceTypeNames.end());  
 }  
   
 void RipleyDomain::setFunctionSpaceTypeNames()  
 {  
     if (m_functionSpaceTypeNames.size() > 0)  
         return;  
   
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Ripley_DegreesOfFreedom"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Ripley_ReducedDegreesOfFreedom"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(Nodes,"Ripley_Nodes"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(ReducedNodes,"Ripley_Reduced_Nodes"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(Elements,"Ripley_Elements"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(ReducedElements,"Ripley_Reduced_Elements"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(FaceElements,"Ripley_Face_Elements"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Ripley_Reduced_Face_Elements"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(Points,"Ripley_Points"));  
 }  
   
 int RipleyDomain::getContinuousFunctionCode() const  
 {  
     return Nodes;  
 }  
   
 int RipleyDomain::getReducedContinuousFunctionCode() const  
 {  
     return ReducedNodes;  
 }  
   
 int RipleyDomain::getFunctionCode() const  
 {  
     return Elements;  
 }  
   
 int RipleyDomain::getReducedFunctionCode() const  
 {  
     return ReducedElements;  
 }  
   
 int RipleyDomain::getFunctionOnBoundaryCode() const  
 {  
     return FaceElements;  
 }  
   
 int RipleyDomain::getReducedFunctionOnBoundaryCode() const  
 {  
     return ReducedFaceElements;  
 }  
   
 int RipleyDomain::getFunctionOnContactZeroCode() const  
 {  
     throw RipleyException("Ripley does not support contact elements");  
 }  
   
 int RipleyDomain::getReducedFunctionOnContactZeroCode() const  
 {  
     throw RipleyException("Ripley does not support contact elements");  
 }  
   
 int RipleyDomain::getFunctionOnContactOneCode() const  
 {  
     throw RipleyException("Ripley does not support contact elements");  
 }  
   
 int RipleyDomain::getReducedFunctionOnContactOneCode() const  
 {  
     throw RipleyException("Ripley does not support contact elements");  
 }  
   
 int RipleyDomain::getSolutionCode() const  
 {  
     return DegreesOfFreedom;  
 }  
   
 int RipleyDomain::getReducedSolutionCode() const  
 {  
     return ReducedDegreesOfFreedom;  
 }  
   
 int RipleyDomain::getDiracDeltaFunctionsCode() const  
 {  
     return Points;  
 }  
   
 //  
 // returns the spatial dimension of the mesh  
 //  
 int RipleyDomain::getDim() const  
 {  
     return m_nodes->getNumDim();  
 }  
   
 //  
 // returns the number of data points summed across all MPI processes  
 //  
 int RipleyDomain::getNumDataPointsGlobal() const  
 {  
     return m_nodes->getGlobalNumNodes();  
59  }  }
60    
61  //  bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
 // returns the number of data points per sample and the number of samples  
 // needed to represent data on a parts of the mesh.  
 //  
 pair<int,int> RipleyDomain::getDataShape(int functionSpaceCode) const  
62  {  {
63      int numDataPointsPerSample=0;      switch (fsType) {
     int numSamples=0;  
     switch (functionSpaceCode) {  
         case Nodes:  
             numDataPointsPerSample=1;  
             numSamples=m_nodes->getNumNodes();  
             break;  
         case ReducedNodes:  
             numDataPointsPerSample=1;  
             numSamples=m_nodes->getNumReducedNodes();  
             break;  
         case Elements:  
             numSamples=m_elements->getNumElements();  
             numDataPointsPerSample=m_elements->getNumLocalDim()+1;  
             break;  
         case ReducedElements:  
             numSamples=m_elements->getNumElements();  
             numDataPointsPerSample=(m_elements->getNumLocalDim()==0)?0:1;  
             break;  
         case FaceElements:  
             numDataPointsPerSample=m_faceElements->getNumLocalDim()+1;  
             numSamples=m_faceElements->getNumElements();  
             break;  
         case ReducedFaceElements:  
             numDataPointsPerSample=(m_faceElements->getNumLocalDim()==0)?0:1;  
             numSamples=m_faceElements->getNumElements();  
             break;  
         case Points:  
             numDataPointsPerSample=1;  
             numSamples=m_points->getNumElements();  
             break;  
64          case DegreesOfFreedom:          case DegreesOfFreedom:
             numDataPointsPerSample=1;  
             numSamples=m_nodes->getNumDegreesOfFreedom();  
             break;  
65          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
             numDataPointsPerSample=1;  
             numSamples=m_nodes->getNumReducedDegreesOfFreedom();  
             break;  
         default:  
             stringstream temp;  
             temp << "Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();  
             throw RipleyException(temp.str());  
     }  
     return pair<int,int>(numDataPointsPerSample,numSamples);  
 }  
   
 //  
 // adds linear PDE of second order into a given stiffness matrix and right hand side  
 //  
 void RipleyDomain::addPDEToSystem(  
         escript::AbstractSystemMatrix& mat, escript::Data& rhs,  
         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& d, const escript::Data& y,  
         const escript::Data& d_contact, const escript::Data& y_contact,  
         const escript::Data& d_dirac,const escript::Data& y_dirac) const  
 {  
     if (!d_contact.isEmpty() || !y_contact.isEmpty())  
         throw RipleyException("Ripley does not support contact elements");  
   
     SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);  
     if (smat==NULL)  
         throw RipleyException("Ripley only accepts its own system matrices");  
 /*  
     escriptDataC _rhs=rhs.getDataC();  
     escriptDataC _A =A.getDataC();  
     escriptDataC _B=B.getDataC();  
     escriptDataC _C=C.getDataC();  
     escriptDataC _D=D.getDataC();  
     escriptDataC _X=X.getDataC();  
     escriptDataC _Y=Y.getDataC();  
     escriptDataC _d=d.getDataC();  
     escriptDataC _y=y.getDataC();  
     escriptDataC _d_dirac=d_dirac.getDataC();  
     escriptDataC _y_dirac=y_dirac.getDataC();  
   
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     Ripley_Assemble_PDE(m_nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y);  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y);  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes,mesh->Points, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d_dirac, 0, &_y_dirac);  
     checkPasoError();  
 */  
 }  
   
 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,  
         const escript::Data& D, const escript::Data& d,  
         const escript::Data& d_dirac, const bool useHRZ) const  
 {  
 /*  
     escriptDataC _mat=mat.getDataC();  
     escriptDataC _D=D.getDataC();  
     escriptDataC _d=d.getDataC();  
     escriptDataC _d_dirac=d_dirac.getDataC();  
   
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     Ripley_Assemble_LumpedSystem(m_nodes, mesh->Elements, &_mat, &_D, useHRZ);  
     checkPasoError();  
   
     Ripley_Assemble_LumpedSystem(m_nodes, mesh->FaceElements, &_mat, &_d, useHRZ);  
     checkPasoError();  
   
     Ripley_Assemble_LumpedSystem(m_nodes, mesh->FaceElements, &_mat, &_d_dirac, useHRZ);  
     checkPasoError();  
 */  
 }  
   
   
 //  
 // adds linear PDE of second order into the right hand side only  
 //  
 void RipleyDomain::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  
 {  
     if (!y_contact.isEmpty())  
         throw RipleyException("Ripley does not support contact elements");  
 /*  
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
   
     escriptDataC _rhs=rhs.getDataC();  
     escriptDataC _X=X.getDataC();  
     escriptDataC _Y=Y.getDataC();  
     escriptDataC _y=y.getDataC();  
     escriptDataC _y_dirac=y_dirac.getDataC();  
     Ripley_Assemble_PDE(m_nodes, mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes, mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes, mesh->Points, 0, &_rhs, 0, 0, 0, 0, 0, &_y_dirac );  
     checkPasoError();  
 */  
 }  
   
 //  
 // adds PDE of second order into a transport problem  
 //  
 void RipleyDomain::addPDEToTransportProblem(  
         escript::AbstractTransportProblem& tp,  
         escript::Data& source, const escript::Data& M,  
         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& d, const escript::Data& y,  
         const escript::Data& d_contact, const escript::Data& y_contact,  
         const escript::Data& d_dirac, const escript::Data& y_dirac) const  
 {  
     if (!d_contact.isEmpty() || !y_contact.isEmpty())  
         throw RipleyException("Ripley does not support contact elements");  
     TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);  
     if (tpa==NULL)  
         throw RipleyException("Ripley only accepts its own transport problems");  
 /*  
     DataTypes::ShapeType shape;  
     source.expand();  
     escriptDataC _source=source.getDataC();  
     escriptDataC _M=M.getDataC();  
     escriptDataC _A=A.getDataC();  
     escriptDataC _B=B.getDataC();  
     escriptDataC _C=C.getDataC();  
     escriptDataC _D=D.getDataC();  
     escriptDataC _X=X.getDataC();  
     escriptDataC _Y=Y.getDataC();  
     escriptDataC _d=d.getDataC();  
     escriptDataC _y=y.getDataC();  
     escriptDataC _d_dirac=d_dirac.getDataC();  
     escriptDataC _y_dirac=y_dirac.getDataC();  
   
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();  
     Ripley_Assemble_PDE(m_nodes, mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0);  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes, mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y);  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes, mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y);  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes, mesh->Points, _tp->transport_matrix, &_source, 0, 0, 0, &_d_dirac, 0, &_y_dirac);  
     checkPasoError();  
 */  
 }  
   
 //  
 // interpolates data between different function spaces  
 //  
 void RipleyDomain::interpolateOnDomain(escript::Data& target,  
                                       const escript::Data& in) const  
 {  
     const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));  
     const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));  
     if (inDomain != *this)  
         throw RipleyException("Illegal domain of interpolant");  
     if (targetDomain != *this)  
         throw RipleyException("Illegal domain of interpolation target");  
 /*  
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     escriptDataC _target=target.getDataC();  
     escriptDataC _in=in.getDataC();  
   
     switch (in.getFunctionSpace().getTypeCode()) {  
66          case Nodes:          case Nodes:
             switch (target.getFunctionSpace().getTypeCode()) {  
                 case Nodes:  
                 case ReducedNodes:  
                 case DegreesOfFreedom:  
                 case ReducedDegreesOfFreedom:  
                     Ripley_Assemble_CopyNodalData(m_nodes, &_target, &_in);  
                     break;  
                 case Elements:  
                 case ReducedElements:  
                     Ripley_Assemble_interpolate(m_nodes, mesh->Elements, &_in, &_target);  
                     break;  
                 case FaceElements:  
                 case ReducedFaceElements:  
                     Ripley_Assemble_interpolate(m_nodes, mesh->FaceElements, &_in, &_target);  
                     break;  
                 case Points:  
                     Ripley_Assemble_interpolate(m_nodes, mesh->Points, &_in, &_target);  
                     break;  
                 default:  
                     stringstream temp;  
                     temp << "Interpolation on Domain: Ripley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
                     throw RipleyException(temp.str());  
             }  
             break;  
67          case ReducedNodes:          case ReducedNodes:
             switch (target.getFunctionSpace().getTypeCode()) {  
                 case Nodes:  
                 case ReducedNodes:  
                 case DegreesOfFreedom:  
                 case ReducedDegreesOfFreedom:  
                     Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);  
                     break;  
                 case Elements:  
                 case ReducedElements:  
                     Ripley_Assemble_interpolate(m_nodes,mesh->Elements,&_in,&_target);  
                     break;  
                 case FaceElements:  
                 case ReducedFaceElements:  
                     Ripley_Assemble_interpolate(m_nodes,mesh->FaceElements,&_in,&_target);  
                     break;  
                 case Points:  
                     Ripley_Assemble_interpolate(m_nodes,mesh->Points,&_in,&_target);  
                     break;  
                 default:  
                     stringstream temp;  
                     temp << "Interpolation on Domain: Ripley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
                     throw RipleyException(temp.str());  
             }  
             break;  
68          case Elements:          case Elements:
             if (target.getFunctionSpace().getTypeCode()==Elements) {  
                 Ripley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
             } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
                 Ripley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);  
             } else {  
                 throw RipleyException("No interpolation with data on elements possible.");  
             }  
             break;  
69          case ReducedElements:          case ReducedElements:
             if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
                 Ripley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
             } else {  
                 throw RipleyException("No interpolation with data on elements with reduced integration order possible.");  
             }  
             break;  
70          case FaceElements:          case FaceElements:
             if (target.getFunctionSpace().getTypeCode()==FaceElements) {  
                 Ripley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
             } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {  
                 Ripley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);  
             } else {  
                 throw RipleyException("No interpolation with data on face elements possible.");  
             }  
             break;  
71          case ReducedFaceElements:          case ReducedFaceElements:
             if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {  
                 Ripley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
             } else {  
                 throw RipleyException("No interpolation with data on face elements with reduced integration order possible.");  
             }  
             break;  
72          case Points:          case Points:
73              if (target.getFunctionSpace().getTypeCode()==Points) {              return true;
                 Ripley_Assemble_CopyElementData(mesh->Points,&_target,&_in);  
             } else {  
                 throw RipleyException("No interpolation with data on points possible.");  
             }  
             break;  
         case DegreesOfFreedom:  
             switch (target.getFunctionSpace().getTypeCode()) {  
                 case ReducedDegreesOfFreedom:  
                 case DegreesOfFreedom:  
                     Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);  
                     break;  
   
                 case Nodes:  
                 case ReducedNodes:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in);  
                         temp.expand();  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in2);  
                     } else {  
                         Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);  
                     }  
                     break;  
                 case Elements:  
                 case ReducedElements:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in, continuousFunction(*this) );  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_interpolate(m_nodes, mesh->Elements,&_in2,&_target);  
                     } else {  
                         Ripley_Assemble_interpolate(m_nodes, mesh->Elements,&_in,&_target);  
                     }  
                     break;  
                 case FaceElements:  
                 case ReducedFaceElements:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in, continuousFunction(*this) );  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_interpolate(m_nodes, mesh->FaceElements,&_in2,&_target);  
                     } else {  
                         Ripley_Assemble_interpolate(m_nodes, mesh->FaceElements,&_in,&_target);  
                     }  
                     break;  
                 case Points:  
                     if (getMPISize()>1) {  
                         //escript::Data temp=escript::Data(in, continuousFunction(*this) );  
                         //escriptDataC _in2 = temp.getDataC();  
                     } else {  
                         Ripley_Assemble_interpolate(m_nodes,mesh->Points,&_in,&_target);  
                     }  
                     break;  
                 default:  
                     stringstream temp;  
                     temp << "Interpolation On Domain: Ripley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
                     throw RipleyException(temp.str());  
             }  
             break;  
         case ReducedDegreesOfFreedom:  
             switch (target.getFunctionSpace().getTypeCode()) {  
                 case Nodes:  
                     throw RipleyException("Ripley does not support interpolation from reduced degrees of freedom to mesh nodes.");  
                     break;  
                 case ReducedNodes:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in);  
                         temp.expand();  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in2);  
                     } else {  
                         Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);  
                     }  
                     break;  
                 case DegreesOfFreedom:  
                     throw RipleyException("Ripley does not support interpolation from reduced degrees of freedom to degrees of freedom");  
                 case ReducedDegreesOfFreedom:  
                     Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);  
                     break;  
                 case Elements:  
                 case ReducedElements:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in, reducedContinuousFunction(*this) );  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_interpolate(m_nodes, mesh->Elements,&_in2,&_target);  
                     } else {  
                         Ripley_Assemble_interpolate(m_nodes, mesh->Elements,&_in,&_target);  
                     }  
                     break;  
                 case FaceElements:  
                 case ReducedFaceElements:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in, reducedContinuousFunction(*this) );  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_interpolate(m_nodes,mesh->FaceElements,&_in2,&_target);  
                     } else {  
                         Ripley_Assemble_interpolate(m_nodes,mesh->FaceElements,&_in,&_target);  
                     }  
                     break;  
                 case Points:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in, reducedContinuousFunction(*this) );  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_interpolate(m_nodes, mesh->Points,&_in2,&_target);  
                     } else {  
                         Ripley_Assemble_interpolate(m_nodes, mesh->Points,&_in,&_target);  
                     }  
                     break;  
                 default:  
                     stringstream temp;  
                     temp << "Interpolation On Domain: Ripley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
                     throw RipleyException(temp.str());  
             }  
             break;  
74          default:          default:
75              stringstream temp;              break;
             temp << "Interpolation On Domain: Ripley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();  
             throw RipleyException(temp.str());  
     }  
 */  
     checkPasoError();  
 }  
   
 //  
 // copies the locations of sample points into arg:  
 //  
 void RipleyDomain::setToX(escript::Data& arg) const  
 {  
     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(*(arg.getFunctionSpace().getDomain()));  
     if (argDomain != *this)  
         throw RipleyException("setToX: Illegal domain of data point locations");  
   
     // do we need to interpolate?  
     if (arg.getFunctionSpace().getTypeCode()==Nodes) {  
         m_nodes->assembleCoordinates(arg);  
     } else {  
         escript::Data tmp_data=Vector(0.0, continuousFunction(*this), true);  
         m_nodes->assembleCoordinates(tmp_data);  
         // this is then interpolated onto arg:  
         interpolateOnDomain(arg, tmp_data);  
76      }      }
77        return false;
78  }  }
79    
80  //  string RipleyDomain::functionSpaceTypeAsString(int fsType) const
 // returns the normal vectors at the location of data points as a Data object  
 //  
 void RipleyDomain::setToNormal(escript::Data& normal) const  
81  {  {
82  /*      switch (fsType) {
83      const RipleyDomain& normalDomain=dynamic_cast<const RipleyDomain&>(*(normal.getFunctionSpace().getDomain()));          case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
84      if (normalDomain!=*this)          case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
85          throw RipleyException("Illegal domain of normal locations");          case Nodes: return "Ripley_Nodes";
86      Ripley_Mesh* mesh=m_ripleyMesh.get();          case ReducedNodes: return "Ripley_Reduced_Nodes";
87      escriptDataC _normal=normal.getDataC();          case Elements: return "Ripley_Elements";
88      switch(normal.getFunctionSpace().getTypeCode()) {          case ReducedElements: return "Ripley_Reduced_Elements";
89          case Nodes:          case FaceElements: return "Ripley_Face_Elements";
90              throw RipleyException("Ripley does not support surface normal vectors for nodes");          case ReducedFaceElements: return "Ripley_Reduced_Face_Elements";
91          case ReducedNodes:          case Points: return "Ripley_Points";
             throw RipleyException("Ripley does not support surface normal vectors for reduced nodes");  
         case Elements:  
             throw RipleyException("Ripley does not support surface normal vectors for elements");  
         case ReducedElements:  
             throw RipleyException("Ripley does not support surface normal vectors for elements with reduced integration order");  
         case FaceElements:  
             Ripley_Assemble_setNormal(m_nodes,mesh->FaceElements,&_normal);  
             break;  
         case ReducedFaceElements:  
             Ripley_Assemble_setNormal(m_nodes,mesh->FaceElements,&_normal);  
             break;  
         case Points:  
             throw RipleyException("Ripley does not support surface normal vectors for point elements");  
         case DegreesOfFreedom:  
             throw RipleyException("Ripley does not support surface normal vectors for degrees of freedom.");  
         case ReducedDegreesOfFreedom:  
             throw RipleyException("Ripley does not support surface normal vectors for reduced degrees of freedom.");  
92          default:          default:
93              stringstream temp;              break;
             temp << "Normal Vectors: Ripley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();  
             throw RipleyException(temp.str());  
94      }      }
95  */      return "Invalid function space type code";
 }  
   
 //  
 // interpolates data to other domain  
 //  
 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const  
 {  
     escript::const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();  
     const RipleyDomain* targetDomain=dynamic_cast<const RipleyDomain*>(targetDomain_p.get());  
     if (targetDomain!=this)  
         throw RipleyException("Illegal domain of interpolation target");  
   
     throw RipleyException("Ripley does not allow interpolation across domains yet");  
96  }  }
97    
98  //  pair<int,int> RipleyDomain::getDataShape(int fsType) const
 // calculates the integral of a function arg  
 //  
 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const  
99  {  {
100  /*      const int ptsPerSample = (m_numDim==2 ? 4 : 8);
101      const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(*(arg.getFunctionSpace().getDomain()));      switch (fsType) {
     if (argDomain!=*this)  
         throw RipleyException("Illegal domain of integration kernel");  
   
     double blocktimer_start = blocktimer_time();  
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     escriptDataC _temp;  
     escript::Data temp;  
     escriptDataC _arg=arg.getDataC();  
     switch(arg.getFunctionSpace().getTypeCode()) {  
102          case Nodes:          case Nodes:
103          case ReducedNodes:          case ReducedNodes: //FIXME: reduced
104                return pair<int,int>(1, getNumNodes());
105          case DegreesOfFreedom:          case DegreesOfFreedom:
106          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom: //FIXME: reduced
107              temp=escript::Data( arg, escript::function(*this) );              return pair<int,int>(1, getNumDOF());
             _temp=temp.getDataC();  
             Ripley_Assemble_integrate(m_nodes,mesh->Elements,&_temp,&integrals[0]);  
             break;  
108          case Elements:          case Elements:
109          case ReducedElements:              return pair<int,int>(ptsPerSample, getNumElements());
             Ripley_Assemble_integrate(m_nodes,mesh->Elements,&_arg,&integrals[0]);  
             break;  
110          case FaceElements:          case FaceElements:
111          case ReducedFaceElements:              return pair<int,int>(ptsPerSample/2, getNumFaceElements());
             Ripley_Assemble_integrate(m_nodes,mesh->FaceElements,&_arg,&integrals[0]);  
             break;  
         case Points:  
             throw RipleyException("Integral of data on points is not supported.");  
         default:  
             stringstream temp;  
             temp << "Integrals: Ripley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();  
             throw RipleyException(temp.str());  
     }  
     blocktimer_increment("integrate()", blocktimer_start);  
 */  
 }  
   
 //  
 // calculates the gradient of arg  
 //  
 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const  
 {  
 /*  
     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(*(arg.getFunctionSpace().getDomain()));  
     if (argDomain!=*this)  
         throw RipleyException("Illegal domain of gradient argument");  
     const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(*(grad.getFunctionSpace().getDomain()));  
     if (gradDomain!=*this)  
         throw RipleyException("Illegal domain of gradient");  
   
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     escriptDataC _grad=grad.getDataC();  
     escriptDataC nodeDataC;  
     escript::Data temp;  
     if (getMPISize()>1) {  
         if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {  
             temp=escript::Data(arg, continuousFunction(*this) );  
             nodeDataC = temp.getDataC();  
         } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {  
             temp=escript::Data(arg, reducedContinuousFunction(*this) );  
             nodeDataC = temp.getDataC();  
         } else {  
             nodeDataC = arg.getDataC();  
         }  
     } else {  
         nodeDataC = arg.getDataC();  
     }  
     switch(grad.getFunctionSpace().getTypeCode()) {  
         case Nodes:  
             throw RipleyException("Gradient at nodes is not supported.");  
         case ReducedNodes:  
             throw RipleyException("Gradient at reduced nodes is not supported.");  
         case Elements:  
             Ripley_Assemble_gradient(m_nodes,mesh->Elements,&_grad,&nodeDataC);  
             break;  
112          case ReducedElements:          case ReducedElements:
113              Ripley_Assemble_gradient(m_nodes,mesh->Elements,&_grad,&nodeDataC);              return pair<int,int>(1, getNumElements());
             break;  
         case FaceElements:  
             Ripley_Assemble_gradient(m_nodes,mesh->FaceElements,&_grad,&nodeDataC);  
             break;  
114          case ReducedFaceElements:          case ReducedFaceElements:
115              Ripley_Assemble_gradient(m_nodes,mesh->FaceElements,&_grad,&nodeDataC);              return pair<int,int>(1, getNumFaceElements());
             break;  
116          case Points:          case Points:
117              throw RipleyException("Gradient at points is not supported");              return pair<int,int>(1, 0); //FIXME: dirac
         case DegreesOfFreedom:  
             throw RipleyException("Gradient at degrees of freedom is not supported");  
         case ReducedDegreesOfFreedom:  
             throw RipleyException("Gradient at reduced degrees of freedom is not supported");  
118          default:          default:
             stringstream temp;  
             temp << "Gradient: Ripley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();  
             throw RipleyException(temp.str());  
     }  
 */  
 }  
   
 //  
 // returns the size of elements  
 //  
 void RipleyDomain::setToSize(escript::Data& size) const  
 {  
 /*  
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     escriptDataC tmp=size.getDataC();  
     switch(size.getFunctionSpace().getTypeCode()) {  
         case Nodes:  
             throw RipleyException("Size of nodes is not supported");  
         case ReducedNodes:  
             throw RipleyException("Size of reduced nodes is not supported");  
         case Elements:  
             Ripley_Assemble_getSize(m_nodes, mesh->Elements, &tmp);  
             break;  
         case ReducedElements:  
             Ripley_Assemble_getSize(m_nodes, mesh->Elements, &tmp);  
119              break;              break;
         case FaceElements:  
             Ripley_Assemble_getSize(m_nodes, mesh->FaceElements, &tmp);  
             break;  
         case ReducedFaceElements:  
             Ripley_Assemble_getSize(m_nodes, mesh->FaceElements, &tmp);  
             break;  
         case Points:  
             throw RipleyException("Size of point elements is not supported");  
         case DegreesOfFreedom:  
             throw RipleyException("Size of degrees of freedom is not supported");  
         case ReducedDegreesOfFreedom:  
             throw RipleyException("Size of reduced degrees of freedom is not supported");  
         default:  
             stringstream temp;  
             temp << "Element size: Ripley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();  
             throw RipleyException(temp.str());  
     }  
 */  
 }  
   
 //  
 // sets the location of nodes  
 //  
 void RipleyDomain::setNewX(const escript::Data& new_x)  
 {  
     const RipleyDomain& newDomain=dynamic_cast<const RipleyDomain&>(*(new_x.getFunctionSpace().getDomain()));  
     if (newDomain!=*this)  
         throw RipleyException("Illegal domain of new node locations");  
   
     if (new_x.getFunctionSpace() == continuousFunction(*this)) {  
         m_nodes->setCoordinates(new_x);  
     } else {  
         escript::Data new_x_inter=escript::Data(new_x, continuousFunction(*this) );  
         m_nodes->setCoordinates(new_x_inter);  
     }  
 }  
   
 bool RipleyDomain::ownSample(int fs_code, index_t id) const  
 {  
 #ifdef ESYS_MPI  
     index_t myFirstNode, myLastNode, k;  
     if (fs_code == ReducedNodes) {  
         myFirstNode = m_nodes->getFirstReducedNode();  
         myLastNode = m_nodes->getLastReducedNode();  
         k = m_nodes->getIndexForGlobalReducedNode(id);  
     } else {  
         myFirstNode = m_nodes->getFirstNode();  
         myLastNode = m_nodes->getLastNode();  
         k = m_nodes->getIndexForGlobalNode(id);  
120      }      }
     return static_cast<bool>((myFirstNode <= k) && (k < myLastNode));  
 #endif  
     return true;  
 }  
   
   
 //  
 // creates a SystemMatrixAdapter stiffness matrix and initializes it with zeros  
 //  
 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,  
         const escript::FunctionSpace& row_functionspace,  
         const int column_blocksize,  
         const escript::FunctionSpace& column_functionspace,  
         const int type) const  
 {  
     bool reduceRowOrder=false;  
     bool reduceColOrder=false;  
     // is the domain right?  
     const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));  
     if (row_domain!=*this)  
         throw RipleyException("Domain of row function space does not match the domain of matrix generator");  
     const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));  
     if (col_domain!=*this)  
         throw RipleyException("Domain of columnn function space does not match the domain of matrix generator");  
     // is the function space type right  
     if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {  
         reduceRowOrder=true;  
     } else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)  
         throw RipleyException("Illegal function space type for system matrix rows");  
     if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {  
         reduceColOrder=true;  
     } else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)  
         throw RipleyException("Illegal function space type for system matrix columns");  
   
     // generate matrix:  
     Paso_SystemMatrixPattern* fsystemMatrixPattern=getPattern(reduceRowOrder, reduceColOrder);  
     checkPasoError();  
     Paso_SystemMatrix* fsystemMatrix = Paso_SystemMatrix_alloc(type,  
             fsystemMatrixPattern, row_blocksize, column_blocksize, FALSE);  
     checkPasoError();  
     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);  
     SystemMatrixAdapter* sma = new SystemMatrixAdapter(fsystemMatrix,  
             row_blocksize, row_functionspace, column_blocksize, column_functionspace);  
     return escript::ASM_ptr(sma);  
 }  
   
 //  
 // creates a TransportProblemAdapter  
 //  
 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,  
         const int blocksize, const escript::FunctionSpace& functionspace,  
         const int type) const  
 {  
     int reduceOrder=0;  
     // is the domain right?  
     const RipleyDomain& domain=dynamic_cast<const RipleyDomain&>(*(functionspace.getDomain()));  
     if (domain!=*this)  
         throw RipleyException("Domain of function space does not match the domain of transport problem generator");  
     // is the function space type right  
     if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {  
         reduceOrder=1;  
     } else if (functionspace.getTypeCode()!=DegreesOfFreedom)  
         throw RipleyException("Illegal function space type for system matrix rows");  
121    
122      // generate matrix      stringstream msg;
123      Paso_SystemMatrixPattern* fsystemMatrixPattern = getPattern(reduceOrder, reduceOrder);      msg << "getDataShape(): Unsupported function space type " << fsType
124      Paso_TransportProblem* transportProblem = Paso_TransportProblem_alloc(          << " for " << getDescription();
125              useBackwardEuler, fsystemMatrixPattern, blocksize);      throw RipleyException(msg.str());
     checkPasoError();  
     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);  
     escript::AbstractTransportProblem* atp=new TransportProblemAdapter(  
             transportProblem, useBackwardEuler, blocksize, functionspace);  
     return escript::ATP_ptr(atp);  
126  }  }
127    
128  // returns true if data on the function space is considered as being cell  string RipleyDomain::showTagNames() const
 // centered  
 bool RipleyDomain::isCellOriented(int functionSpaceCode) const  
129  {  {
130      switch(functionSpaceCode) {      stringstream ret;
131          case Nodes:      TagMap::const_iterator it;
132          case DegreesOfFreedom:      for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
133          case ReducedDegreesOfFreedom:          if (it!=m_tagMap.begin()) ret << ", ";
134              return false;          ret << it->first;
         case Elements:  
         case FaceElements:  
         case Points:  
         case ReducedElements:  
         case ReducedFaceElements:  
             return true;  
         default:  
             stringstream temp;  
             temp << "Ripley does not know anything about function space type " << functionSpaceCode;  
             throw RipleyException(temp.str());  
135      }      }
136      return false;      return ret.str();
137  }  }
138    
139  bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const  bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
140  {  {
141     /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]     /*
142      class 1: DOF <-> Nodes      The idea is to use equivalence classes (i.e. types which can be
143      class 2: ReducedDOF <-> ReducedNodes      interpolated back and forth):
144      class 3: Points      class 0: DOF <-> Nodes
145      class 4: Elements      class 1: ReducedDOF <-> ReducedNodes
146      class 5: ReducedElements      class 2: Points
147      class 6: FaceElements      class 3: Elements
148      class 7: ReducedFaceElements      class 4: ReducedElements
149      class 8: ContactElementZero <-> ContactElementOne      class 5: FaceElements
150      class 9: ReducedContactElementZero <-> ReducedContactElementOne      class 6: ReducedFaceElements
151    
152     There is also a set of lines. Interpolation is possible down a line but not between lines.      There is also a set of lines. Interpolation is possible down a line but not
153     class 1 and 2 belong to all lines so aren't considered.      between lines.
154      line 0: class 3      class 0 and 1 belong to all lines so aren't considered.
155      line 1: class 4,5      line 0: class 2
156      line 2: class 6,7      line 1: class 3,4
157      line 3: class 8,9      line 2: class 5,6
158    
159      For classes with multiple members (eg class 2) we have vars to record if      For classes with multiple members (eg class 1) we have vars to record if
160      there is at least one instance. e.g. hasnodes is true if we have at least      there is at least one instance. e.g. hasnodes is true if we have at least
161      one instance of Nodes.      one instance of Nodes.
162      */      */
163      if (fs.empty())      if (fs.empty())
164          return false;          return false;
165      vector<int> hasclass(10);      vector<bool> hasclass(7, false);
166      vector<int> hasline(4);      vector<int> hasline(3, 0);
167      bool hasnodes=false;      bool hasnodes=false;
168      bool hasrednodes=false;      bool hasrednodes=false;
169      for (int i=0;i<fs.size();++i) {      for (size_t i=0; i<fs.size(); ++i) {
170          switch (fs[i]) {          switch (fs[i]) {
171              case Nodes: hasnodes=true; // no break is deliberate              case Nodes: hasnodes=true; // fall through
172              case DegreesOfFreedom:              case DegreesOfFreedom:
173                  hasclass[1]=1;                  hasclass[0]=true;
174                  break;                  break;
175              case ReducedNodes: hasrednodes=true; // no break is deliberate              case ReducedNodes: hasrednodes=true; // fall through
176              case ReducedDegreesOfFreedom:              case ReducedDegreesOfFreedom:
177                  hasclass[2]=1;                  hasclass[1]=true;
178                  break;                  break;
179              case Points:              case Points:
180                  hasline[0]=1;                  hasline[0]=1;
181                  hasclass[3]=1;                  hasclass[2]=true;
182                  break;                  break;
183              case Elements:              case Elements:
184                  hasclass[4]=1;                  hasclass[3]=true;
185                  hasline[1]=1;                  hasline[1]=1;
186                  break;                  break;
187              case ReducedElements:              case ReducedElements:
188                  hasclass[5]=1;                  hasclass[4]=true;
189                  hasline[1]=1;                  hasline[1]=1;
190                  break;                  break;
191              case FaceElements:              case FaceElements:
192                  hasclass[6]=1;                  hasclass[5]=true;
193                  hasline[2]=1;                  hasline[2]=1;
194                  break;                  break;
195              case ReducedFaceElements:              case ReducedFaceElements:
196                  hasclass[7]=1;                  hasclass[6]=true;
197                  hasline[2]=1;                  hasline[2]=1;
198                  break;                  break;
199              default:              default:
200                  return false;                  return false;
201          }          }
202      }      }
203      int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];      int numLines=hasline[0]+hasline[1]+hasline[2];
204    
205      // fail if we have more than one leaf group      // fail if we have more than one leaf group
206      // = there are at least two branches we can't interpolate between      // = there are at least two branches we can't interpolate between
207      if (totlines>1)      if (numLines > 1)
208          return false;          return false;
209      else if (totlines==1) {      else if (numLines==1) {
210          // we have points          // we have points
211          if (hasline[0]==1) {          if (hasline[0]==1)
212              resultcode=Points;              resultcode=Points;
213          } else if (hasline[1]==1) {          else if (hasline[1]==1) {
214              if (hasclass[5]==1) {              if (hasclass[4])
215                  resultcode=ReducedElements;                  resultcode=ReducedElements;
216              } else {              else
217                  resultcode=Elements;                  resultcode=Elements;
218              }          } else { // hasline[2]==1
219          } else if (hasline[2]==1) {              if (hasclass[6])
             if (hasclass[7]==1) {  
220                  resultcode=ReducedFaceElements;                  resultcode=ReducedFaceElements;
221              } else {              else
222                  resultcode=FaceElements;                  resultcode=FaceElements;
             }  
         } else { // so we must be in line3  
             throw RipleyException("Programmer Error - choosing between contact elements - we should never get here");  
223          }          }
224      } else { // totlines==0      } else { // numLines==0
225          if (hasclass[2]==1) {          if (hasclass[1])
226              // something from class 2              // something from class 1
227              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
228          } else { // something from class 1          else // something from class 0
229              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
         }  
230      }      }
231      return true;      return true;
232  }  }
233    
234  bool RipleyDomain::probeInterpolationOnDomain(int functionSpaceType_source,  bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
235                                               int functionSpaceType_target) const                                                int fsType_target) const
236  {  {
237      if (functionSpaceType_target != Nodes &&      if (fsType_target != Nodes &&
238              functionSpaceType_target != ReducedNodes &&              fsType_target != ReducedNodes &&
239              functionSpaceType_target != ReducedDegreesOfFreedom &&              fsType_target != ReducedDegreesOfFreedom &&
240              functionSpaceType_target != DegreesOfFreedom &&              fsType_target != DegreesOfFreedom &&
241              functionSpaceType_target != Elements &&              fsType_target != Elements &&
242              functionSpaceType_target != ReducedElements &&              fsType_target != ReducedElements &&
243              functionSpaceType_target != FaceElements &&              fsType_target != FaceElements &&
244              functionSpaceType_target != ReducedFaceElements &&              fsType_target != ReducedFaceElements &&
245              functionSpaceType_target != Points) {              fsType_target != Points) {
246          stringstream temp;          stringstream msg;
247          temp << "Interpolation On Domain: Ripley does not know anything about function space type " << functionSpaceType_target;          msg << "probeInterpolationOnDomain(): Invalid functionspace type "
248          throw RipleyException(temp.str());              << fsType_target << " for " << getDescription();
249            throw RipleyException(msg.str());
250      }      }
251    
252      switch (functionSpaceType_source) {      switch (fsType_source) {
253          case Nodes:          case Nodes:
254          case DegreesOfFreedom:          case DegreesOfFreedom:
255              return true;              return true;
256          case ReducedNodes:          case ReducedNodes:
257          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
258              return (functionSpaceType_target != Nodes &&              return (fsType_target != Nodes &&
259                      functionSpaceType_target != DegreesOfFreedom);                      fsType_target != DegreesOfFreedom);
260          case Elements:          case Elements:
261              return (functionSpaceType_target==Elements ||              return (fsType_target==Elements ||
262                      functionSpaceType_target==ReducedElements);                      fsType_target==ReducedElements);
263          case ReducedElements:          case ReducedElements:
264              return (functionSpaceType_target==ReducedElements);              return (fsType_target==ReducedElements);
265          case FaceElements:          case FaceElements:
266              return (functionSpaceType_target==FaceElements ||              return (fsType_target==FaceElements ||
267                      functionSpaceType_target==ReducedFaceElements);                      fsType_target==ReducedFaceElements);
268          case ReducedFaceElements:          case ReducedFaceElements:
269              return (functionSpaceType_target==ReducedFaceElements);              return (fsType_target==ReducedFaceElements);
270          case Points:          case Points:
271              return (functionSpaceType_target==Points);              return (fsType_target==Points);
272    
273          default:          default: {
274              stringstream temp;              stringstream msg;
275              temp << "Interpolation On Domain: Ripley does not know anything about function space type " << functionSpaceType_source;              msg << "probeInterpolationOnDomain(): Invalid functionspace type "
276              throw RipleyException(temp.str());                  << fsType_source << " for " << getDescription();
277                throw RipleyException(msg.str());
278            }
279      }      }
     return false;  
280  }  }
281    
282  bool RipleyDomain::probeInterpolationACross(int functionSpaceType_source,  void RipleyDomain::interpolateOnDomain(escript::Data& target,
283          const AbstractDomain& targetDomain, int functionSpaceType_target) const                                         const escript::Data& in) const
284  {  {
285      return false;      const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
286  }      const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
287        if (inDomain != *this)
288            throw RipleyException("Illegal domain of interpolant");
289        if (targetDomain != *this)
290            throw RipleyException("Illegal domain of interpolation target");
291    
292  bool RipleyDomain::operator==(const AbstractDomain& other) const      stringstream msg;
293  {      msg << "interpolateOnDomain() not implemented for function space "
294      const RipleyDomain* temp=dynamic_cast<const RipleyDomain*>(&other);          << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
295      if (temp != NULL) {          << " -> "
296          return (m_name==temp->m_name && m_nodes==temp->m_nodes &&          << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
297                  m_elements==temp->m_elements &&  
298                  m_faceElements==temp->m_faceElements &&      const int inFS = in.getFunctionSpace().getTypeCode();
299                  m_points==temp->m_points);      const int outFS = target.getFunctionSpace().getTypeCode();
300    
301        // simplest case: 1:1 copy
302        if (inFS==outFS) {
303            copyData(target, *const_cast<escript::Data*>(&in));
304        // not allowed: reduced->non-reduced
305        } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
306                && (outFS==Nodes || outFS==DegreesOfFreedom)) {
307            throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");
308        } else if ((inFS==Elements && outFS==ReducedElements)
309                || (inFS==FaceElements && outFS==ReducedFaceElements)) {
310            averageData(target, *const_cast<escript::Data*>(&in));
311        } else {
312            switch (inFS) {
313                case Nodes:
314                case ReducedNodes: //FIXME: reduced
315                    switch (outFS) {
316                        case DegreesOfFreedom:
317                        case ReducedDegreesOfFreedom: //FIXME: reduced
318                            if (getMPISize()==1)
319                                copyData(target, *const_cast<escript::Data*>(&in));
320                            else
321                                nodesToDOF(target,*const_cast<escript::Data*>(&in));
322                            break;
323    
324                        case Nodes:
325                        case ReducedNodes: //FIXME: reduced
326                            copyData(target, *const_cast<escript::Data*>(&in));
327                            break;
328    
329                        case Elements:
330                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
331                            break;
332    
333                        case ReducedElements:
334                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
335                            break;
336    
337                        case FaceElements:
338                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
339                            break;
340    
341                        case ReducedFaceElements:
342                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
343                            break;
344    
345                        default:
346                            throw RipleyException(msg.str());
347                    }
348                    break;
349    
350                case DegreesOfFreedom:
351                case ReducedDegreesOfFreedom: //FIXME: reduced
352                    switch (outFS) {
353                        case Nodes:
354                        case ReducedNodes: //FIXME: reduced
355                            if (getMPISize()==1)
356                                copyData(target, *const_cast<escript::Data*>(&in));
357                            else
358                                dofToNodes(target, *const_cast<escript::Data*>(&in));
359                            break;
360    
361                        case DegreesOfFreedom:
362                        case ReducedDegreesOfFreedom: //FIXME: reduced
363                            copyData(target, *const_cast<escript::Data*>(&in));
364                            break;
365    
366                        case Elements:
367                        case ReducedElements:
368                            if (getMPISize()==1) {
369                                interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
370                            } else {
371                                escript::Data contIn=escript::Data(in, continuousFunction(*this));
372                                interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
373                            }
374                            break;
375    
376                        case FaceElements:
377                        case ReducedFaceElements:
378                            if (getMPISize()==1) {
379                                interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
380                            } else {
381                                escript::Data contIn=escript::Data(in, continuousFunction(*this));
382                                interpolateNodesOnElements(target, contIn, outFS==ReducedFaceElements);
383                            }
384                            break;
385    
386                        default:
387                            throw RipleyException(msg.str());
388                    }
389                    break;
390                default:
391                    throw RipleyException(msg.str());
392            }
393      }      }
     return false;  
394  }  }
395    
396  bool RipleyDomain::operator!=(const AbstractDomain& other) const  escript::Data RipleyDomain::getX() const
397  {  {
398      return !(operator==(other));      return escript::continuousFunction(*this).getX();
399  }  }
400    
401  int RipleyDomain::getSystemMatrixTypeId(const int solver,  escript::Data RipleyDomain::getNormal() const
         const int preconditioner, const int package, const bool symmetry) const  
402  {  {
403      int out=Paso_SystemMatrix_getSystemMatrixTypeId(      return escript::functionOnBoundary(*this).getNormal();
             SystemMatrixAdapter::mapOptionToPaso(solver),  
             SystemMatrixAdapter::mapOptionToPaso(preconditioner),  
             SystemMatrixAdapter::mapOptionToPaso(package),  
             symmetry?1:0, m_mpiInfo);  
     checkPasoError();  
     return out;  
404  }  }
405    
406  int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,  escript::Data RipleyDomain::getSize() const
         const int package, const bool symmetry) const  
407  {  {
408      int out=Paso_TransportProblem_getTypeId(      return escript::function(*this).getSize();
             SystemMatrixAdapter::mapOptionToPaso(solver),  
             SystemMatrixAdapter::mapOptionToPaso(preconditioner),  
             SystemMatrixAdapter::mapOptionToPaso(package),  
             symmetry?1:0, m_mpiInfo);  
     checkPasoError();  
     return out;  
409  }  }
410    
411  escript::Data RipleyDomain::getX() const  void RipleyDomain::setToX(escript::Data& arg) const
412  {  {
413      return continuousFunction(*this).getX();      const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
414  }              *(arg.getFunctionSpace().getDomain()));
415        if (argDomain != *this)
416            throw RipleyException("setToX: Illegal domain of data point locations");
417        if (!arg.isExpanded())
418            throw RipleyException("setToX: Expanded Data object expected");
419    
420  escript::Data RipleyDomain::getNormal() const      if (arg.getFunctionSpace().getTypeCode()==Nodes) {
421  {          assembleCoordinates(arg);
422      return functionOnBoundary(*this).getNormal();      } else {
423            // interpolate the result
424            escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
425            assembleCoordinates(contData);
426            interpolateOnDomain(arg, contData);
427        }
428  }  }
429    
430  escript::Data RipleyDomain::getSize() const  bool RipleyDomain::isCellOriented(int fsType) const
431  {  {
432      return escript::function(*this).getSize();      switch(fsType) {
433            case Nodes:
434            case DegreesOfFreedom:
435            case ReducedDegreesOfFreedom:
436                return false;
437            case Elements:
438            case FaceElements:
439            case Points:
440            case ReducedElements:
441            case ReducedFaceElements:
442                return true;
443            default:
444                break;
445        }
446        stringstream msg;
447        msg << "isCellOriented(): Illegal function space type " << fsType
448            << " on " << getDescription();
449        throw RipleyException(msg.str());
450  }  }
451    
452  const int* RipleyDomain::borrowSampleReferenceIDs(int functionSpaceType) const  bool RipleyDomain::canTag(int fsType) const
453  {  {
454      switch (functionSpaceType) {      switch(fsType) {
455          case Nodes:          case Nodes:
             return &m_nodes->getIdVector()[0];  
         case ReducedNodes:  
             return &m_nodes->getReducedNodesId()[0];  
456          case Elements:          case Elements:
457          case ReducedElements:          case ReducedElements:
             return &m_elements->getIdVector()[0];  
458          case FaceElements:          case FaceElements:
459          case ReducedFaceElements:          case ReducedFaceElements:
460              return &m_faceElements->getIdVector()[0];              return true;
         case Points:  
             return &m_points->getIdVector()[0];  
461          case DegreesOfFreedom:          case DegreesOfFreedom:
             return &m_nodes->getDegreesOfFreedomId()[0];  
462          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
463              return &m_nodes->getReducedDegreesOfFreedomId()[0];          case Points:
464                return false;
465          default:          default:
466              stringstream msg;              break;
             msg << "borrowSampleReferenceIDs: Invalid function space type " << functionSpaceType << " for domain: " << getDescription();  
             throw RipleyException(msg.str());  
467      }      }
468        stringstream msg;
469        msg << "canTag(): Illegal function space type " << fsType << " on "
470            << getDescription();
471        throw RipleyException(msg.str());
472  }  }
473    
474  int RipleyDomain::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& cMask) const
475  {  {
476      throw RipleyException("not implemented");      IndexVector* target=NULL;
477  }      dim_t num=0;
478    
479        switch(fsType) {
 void RipleyDomain::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const  
 {  
     switch (functionSpaceType) {  
480          case Nodes:          case Nodes:
481              m_nodes->setTags(newTag, mask);              num=getNumNodes();
482                target=&m_nodeTags;
483              break;              break;
         case ReducedNodes:  
             throw RipleyException("ReducedNodes does not support tags");  
         case DegreesOfFreedom:  
             throw RipleyException("DegreesOfFreedom does not support tags");  
         case ReducedDegreesOfFreedom:  
             throw RipleyException("ReducedDegreesOfFreedom does not support tags");  
484          case Elements:          case Elements:
             m_elements->setTags(newTag, mask);  
             break;  
485          case ReducedElements:          case ReducedElements:
486              m_elements->setTags(newTag, mask);              num=getNumElements();
487                target=&m_elementTags;
488              break;              break;
489          case FaceElements:          case FaceElements:
             m_faceElements->setTags(newTag, mask);  
             break;  
490          case ReducedFaceElements:          case ReducedFaceElements:
491              m_faceElements->setTags(newTag, mask);              num=getNumFaceElements();
492              break;              target=&m_faceTags;
         case Points:  
             m_points->setTags(newTag, mask);  
493              break;              break;
494          default:          default: {
495              stringstream msg;              stringstream msg;
496              msg << "Ripley does not know anything about function space type " << functionSpaceType;              msg << "setTags(): not implemented for "
497                    << functionSpaceTypeAsString(fsType);
498              throw RipleyException(msg.str());              throw RipleyException(msg.str());
499            }
500        }
501        if (target->size() != num) {
502            target->assign(num, -1);
503      }      }
 }  
   
 void RipleyDomain::setTagMap(const string& name, int tag)  
 {  
     m_tagMap[name] = tag;  
 }  
   
 int RipleyDomain::getTag(const string& name) const  
 {  
     return m_tagMap.find(name)->second;  
 }  
   
 bool RipleyDomain::isValidTagName(const string& name) const  
 {  
     return (m_tagMap.find(name)!=m_tagMap.end());  
 }  
504    
505  string RipleyDomain::showTagNames() const      escript::Data& mask=*const_cast<escript::Data*>(&cMask);
506  {  #pragma omp parallel for
507      stringstream ret;      for (index_t i=0; i<num; i++) {
508      TagMap::const_iterator it;          if (mask.getSampleDataRO(i)[0] > 0) {
509      for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {              (*target)[i]=newTag;
510          if (it!=m_tagMap.begin()) ret << ", ";          }
         ret << it->first;  
511      }      }
512      return ret.str();      updateTagsInUse(fsType);
513  }  }
514    
515  int RipleyDomain::getNumberOfTagsInUse(int functionSpaceCode) const  int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
516  {  {
517      dim_t numTags=0;      switch(fsType) {
     switch (functionSpaceCode) {  
518          case Nodes:          case Nodes:
519              numTags=m_nodes->getNumberOfTagsInUse();              if (m_nodeTags.size() > sampleNo)
520                    return m_nodeTags[sampleNo];
521              break;              break;
         case ReducedNodes:  
             throw RipleyException("ReducedNodes does not support tags");  
         case DegreesOfFreedom:  
             throw RipleyException("DegreesOfFreedom does not support tags");  
         case ReducedDegreesOfFreedom:  
             throw RipleyException("ReducedDegreesOfFreedom does not support tags");  
522          case Elements:          case Elements:
523          case ReducedElements:          case ReducedElements:
524              numTags=m_elements->getNumberOfTagsInUse();              if (m_elementTags.size() > sampleNo)
525                    return m_elementTags[sampleNo];
526              break;              break;
527          case FaceElements:          case FaceElements:
528          case ReducedFaceElements:          case ReducedFaceElements:
529              numTags=m_faceElements->getNumberOfTagsInUse();              if (m_faceTags.size() > sampleNo)
530              break;                  return m_faceTags[sampleNo];
         case Points:  
             numTags=m_points->getNumberOfTagsInUse();  
531              break;              break;
532          default:          default: {
533              stringstream msg;              stringstream msg;
534              msg << "Ripley does not know anything about function space type " << functionSpaceCode;              msg << "getTagFromSampleNo(): not implemented for "
535                    << functionSpaceTypeAsString(fsType);
536              throw RipleyException(msg.str());              throw RipleyException(msg.str());
537            }
538      }      }
539      return numTags;      return -1;
540  }  }
541    
542  const int* RipleyDomain::borrowListOfTagsInUse(int functionSpaceCode) const  int RipleyDomain::getNumberOfTagsInUse(int fsType) const
543  {  {
544      switch (functionSpaceCode) {      switch(fsType) {
545          case Nodes:          case Nodes:
546              return &m_nodes->getTagsInUse()[0];              return m_nodeTagsInUse.size();
         case ReducedNodes:  
             throw RipleyException("ReducedNodes does not support tags");  
         case DegreesOfFreedom:  
             throw RipleyException("DegreesOfFreedom does not support tags");  
         case ReducedDegreesOfFreedom:  
             throw RipleyException("ReducedDegreesOfFreedom does not support tags");  
547          case Elements:          case Elements:
548          case ReducedElements:          case ReducedElements:
549              return &m_elements->getTagsInUse()[0];              return m_elementTagsInUse.size();
550          case FaceElements:          case FaceElements:
551          case ReducedFaceElements:          case ReducedFaceElements:
552              return &m_faceElements->getTagsInUse()[0];              return m_faceTagsInUse.size();
553          case Points:          default: {
             return &m_points->getTagsInUse()[0];  
         default:  
554              stringstream msg;              stringstream msg;
555              msg << "Ripley does not know anything about function space type " << functionSpaceCode;              msg << "getNumberOfTagsInUse(): not implemented for "
556                    << functionSpaceTypeAsString(fsType);
557              throw RipleyException(msg.str());              throw RipleyException(msg.str());
558            }
559      }      }
560  }  }
561    
562  bool RipleyDomain::canTag(int functionSpaceCode) const  const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
563  {  {
564      switch(functionSpaceCode) {      switch(fsType) {
565          case Nodes:          case Nodes:
566                return &m_nodeTagsInUse[0];
567          case Elements:          case Elements:
568          case ReducedElements:          case ReducedElements:
569                return &m_elementTagsInUse[0];
570          case FaceElements:          case FaceElements:
571          case ReducedFaceElements:          case ReducedFaceElements:
572          case Points:              return &m_faceTagsInUse[0];
573              return true;          default: {
574                stringstream msg;
575                msg << "borrowListOfTagsInUse(): not implemented for "
576                    << functionSpaceTypeAsString(fsType);
577                throw RipleyException(msg.str());
578            }
579      }      }
     return false;  
580  }  }
581    
582  escript::AbstractDomain::StatusType RipleyDomain::getStatus() const  void RipleyDomain::Print_Mesh_Info(const bool full) const
583  {  {
584      return m_nodes->getStatus();      cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
585            << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
586        cout << "Number of dimensions: " << m_numDim << endl;
587    
588        // write tags
589        if (m_tagMap.size() > 0) {
590            cout << "Tags:" << endl;
591            TagMap::const_iterator it;
592            for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
593                cout << "  " << setw(5) << it->second << " "
594                    << it->first << endl;
595            }
596        }
597  }  }
598    
599  void RipleyDomain::prepare(bool optimize)  int RipleyDomain::getSystemMatrixTypeId(const int solver,
600            const int preconditioner, const int package, const bool symmetry) const
601  {  {
602      resolveNodeIds();      return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
603                package, symmetry, m_mpiInfo);
604    }
605    
606      // first step is to distribute the elements according to a global  int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
607      // distribution of DOF          const int package, const bool symmetry) const
608      // - create dense labeling for the DOFs  {
609      dim_t newGlobalNumDOFs=m_nodes->createDenseDOFLabeling();      return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
610                package, symmetry, m_mpiInfo);
611      // - create a distribution of the global DOFs and determine  }
     //   the MPI_rank controlling the DOFs on this processor  
     IndexVector distribution(m_mpiInfo->size+1);  
     Esys_MPIInfo_setDistribution(m_mpiInfo, 0, newGlobalNumDOFs-1, &distribution[0]);  
     checkPasoError();  
   
     // Now the mesh is redistributed according to the mpiRankOfDOF vector.  
     // This will redistribute the Nodes and Elements including overlap and will  
     // create an element coloring but will not create any mappings (see later  
     // in this function).  
     distributeByRankOfDOF(distribution);  
   
     // At this stage we are able to start an optimization of the DOF  
     // distribution using ParMetis. On return distribution is altered and new  
     // DOF IDs have been assigned.  
     if (optimize) {  
         if (m_mpiInfo->size > 1) {  
             optimizeDOFDistribution(distribution);  
             distributeByRankOfDOF(distribution);  
         }  
612    
613          // optimize the local labeling of the degrees of freedom  escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
614          optimizeDOFLabeling(distribution);          const escript::FunctionSpace& row_functionspace,
615      }          const int column_blocksize,
616            const escript::FunctionSpace& column_functionspace,
617            const int type) const
618    {
619        bool reduceRowOrder=false;
620        bool reduceColOrder=false;
621        // is the domain right?
622        const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
623        if (row_domain!=*this)
624            throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
625        const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
626        if (col_domain!=*this)
627            throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
628        // is the function space type right?
629        if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
630            reduceRowOrder=true;
631        else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
632            throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
633        if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
634            reduceColOrder=true;
635        else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
636            throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
637    
638      // Rearrange elements in order to try and bring them closer to memory      // generate matrix
639      // locations of the nodes (distributed shared memory).      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
640      optimizeElementOrdering();      Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
641                row_blocksize, column_blocksize, FALSE);
642        paso::checkPasoError();
643        Paso_SystemMatrixPattern_free(pattern);
644        escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
645                    row_functionspace, column_blocksize, column_functionspace));
646        return sma;
647    }
648    
649  /* useful DEBUG:  void RipleyDomain::addPDEToSystem(
650      {          escript::AbstractSystemMatrix& mat, escript::Data& rhs,
651          printf("prepare: global DOF : %d\n",newGlobalNumDOFs);          const escript::Data& A, const escript::Data& B, const escript::Data& C,
652          IndexPair range = m_nodes->getGlobalIdRange();          const escript::Data& D, const escript::Data& X, const escript::Data& Y,
653          printf("prepare: global node id range = %d :%d\n",range.first,range.second);          const escript::Data& d, const escript::Data& y,
654          range = m_nodes->getIdRange();          const escript::Data& d_contact, const escript::Data& y_contact,
655          printf("prepare: local node id range = %d :%d\n",range.first,range.second);          const escript::Data& d_dirac,const escript::Data& y_dirac) const
656      }  {
657  */      if (!d_contact.isEmpty() || !y_contact.isEmpty())
658            throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
659    
660        paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
661        if (!sma)
662            throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
663    
664      // create the global indices      if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
665      IndexVector maskReducedNodes(m_nodes->getNumNodes(), -1);          throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
     markNodes(maskReducedNodes, 0);  
     IndexVector indexReducedNodes = packMask(maskReducedNodes);  
666    
667      IndexVector nodeDistribution(m_mpiInfo->size+1);      //TODO: more input param checks (shape, function space etc)
     m_nodes->createDenseNodeLabeling(nodeDistribution, distribution);  
     m_nodes->createDenseReducedDOFLabeling(maskReducedNodes);  
     m_nodes->createDenseReducedNodeLabeling(maskReducedNodes);  
668    
669      // create the missing mappings      Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
     m_nodes->createNodeFileMappings(indexReducedNodes, distribution,  
                                     nodeDistribution);  
670    
671      updateTagsInUse();      if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
672            throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
673    
674        const int numEq=S->logical_row_block_size;
675        const int numComp=S->logical_col_block_size;
676    
677        if (numEq != numComp)
678            throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
679        //TODO: more system matrix checks
680    
681        if (numEq==1)
682            assemblePDESingle(S, rhs, A, B, C, D, X, Y, d, y);
683        else
684            assemblePDESystem(S, rhs, A, B, C, D, X, Y, d, y);
685  }  }
686    
687  void RipleyDomain::optimizeElementOrdering()  void RipleyDomain::setNewX(const escript::Data& arg)
688  {  {
689      m_elements->optimizeOrdering();      throw RipleyException("setNewX(): Operation not supported");
     m_faceElements->optimizeOrdering();  
     m_points->optimizeOrdering();  
690  }  }
691    
692  void RipleyDomain::updateTagsInUse()  //protected
693    void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
694  {  {
695      m_nodes->updateTagsInUse();      const dim_t numComp = in.getDataPointSize();
696      m_elements->updateTagsInUse();      const dim_t dpp = in.getNumDataPointsPerSample();
697      m_faceElements->updateTagsInUse();      out.requireWrite();
698      m_points->updateTagsInUse();  #pragma omp parallel for
699        for (index_t i=0; i<in.getNumSamples(); i++) {
700            const double* src = in.getSampleDataRO(i);
701            double* dest = out.getSampleDataRW(i);
702            for (index_t c=0; c<numComp; c++) {
703                double res=0.;
704                for (index_t q=0; q<dpp; q++)
705                    res+=src[c+q*numComp];
706                *dest++ = res/dpp;
707            }
708        }
709  }  }
710    
711  void RipleyDomain::createColoring(const IndexVector &node_localDOF_map)  //protected
712    void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
713  {  {
714      m_elements->createColoring(node_localDOF_map);      const dim_t numComp = in.getDataPointSize();
715      m_faceElements->createColoring(node_localDOF_map);      out.requireWrite();
716      m_points->createColoring(node_localDOF_map);  #pragma omp parallel for
717        for (index_t i=0; i<in.getNumSamples(); i++) {
718            const double* src = in.getSampleDataRO(i);
719            copy(src, src+numComp, out.getSampleDataRW(i));
720        }
721  }  }
722    
723  void RipleyDomain::distributeByRankOfDOF(const IndexVector &dofDistribution)  //protected
724    void RipleyDomain::updateTagsInUse(int fsType) const
725  {  {
726      RankVector mpiRankOfDOF = m_nodes->getOwnerOfDOFs(dofDistribution);      IndexVector* tagsInUse=NULL;
727        const IndexVector* tags=NULL;
728        switch(fsType) {
729            case Nodes:
730                tags=&m_nodeTags;
731                tagsInUse=&m_nodeTagsInUse;
732                break;
733            case Elements:
734            case ReducedElements:
735                tags=&m_elementTags;
736                tagsInUse=&m_elementTagsInUse;
737                break;
738            case FaceElements:
739            case ReducedFaceElements:
740                tags=&m_faceTags;
741                tagsInUse=&m_faceTagsInUse;
742                break;
743            default:
744                return;
745        }
746    
747      // First the elements are redistributed according to mpiRankOfDOF.      // gather global unique tag values from tags into tagsInUse
748      // At the input the Node tables refer to the local labeling of      tagsInUse->clear();
749      // the nodes while at the output they refer to the global labeling      index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
     // which is rectified in the next step  
     m_elements->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());  
     m_faceElements->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());  
     m_points->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());  
750    
751      // resolve the node IDs      while (true) {
752      resolveNodeIds();          // find smallest value bigger than lastFoundValue
753            minFoundValue = INDEX_T_MAX;
754    #pragma omp parallel private(local_minFoundValue)
755            {
756                local_minFoundValue = minFoundValue;
757    #pragma omp for schedule(static) nowait
758                for (size_t i = 0; i < tags->size(); i++) {
759                    const index_t v = (*tags)[i];
760                    if ((v > lastFoundValue) && (v < local_minFoundValue))
761                        local_minFoundValue = v;
762                }
763    #pragma omp critical
764                {
765                    if (local_minFoundValue < minFoundValue)
766                        minFoundValue = local_minFoundValue;
767                }
768            }
769    #ifdef ESYS_MPI
770            local_minFoundValue = minFoundValue;
771            MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
772    #endif
773    
774      // create a local labeling of the DOFs          // if we found a new value add it to the tagsInUse vector
775      IndexPair dofRange = m_nodes->getDOFRange();          if (minFoundValue < INDEX_T_MAX) {
776      dim_t len = dofRange.second - dofRange.first + 1;              tagsInUse->push_back(minFoundValue);
777                lastFoundValue = minFoundValue;
778            } else
779                break;
780        }
781    }
782    
783      // local mask for used nodes  //protected
784      IndexVector tmp_node_localDOF_mask(len, -1);  Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
785      IndexVector tmp_node_localDOF_map(m_nodes->getNumNodes(), -1);          const IndexVector& index, const dim_t M, const dim_t N) const
786    {
787        // paso will manage the memory
788        index_t* indexC = MEMALLOC(index.size(), index_t);
789        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
790        copy(index.begin(), index.end(), indexC);
791        copy(ptr.begin(), ptr.end(), ptrC);
792        return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
793    }
794    
795  #pragma omp parallel for schedule(static)  //protected
796      for (dim_t n = 0; n < m_nodes->getNumNodes(); n++) {  Paso_Pattern* RipleyDomain::createMainPattern() const
797  #ifdef BOUNDS_CHECK  {
798          if ((m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first) >= len      IndexVector ptr(1,0);
799                  || (m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first) < 0) {      IndexVector index;
800              printf("BOUNDS_CHECK %s %d\n", __FILE__, __LINE__);  
801              exit(1);      for (index_t i=0; i<getNumDOF(); i++) {
802          }          // add the DOF itself
803  #endif          index.push_back(i);
804          tmp_node_localDOF_mask[m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first] = n;          const dim_t num=insertNeighbourNodes(index, i);
805            ptr.push_back(ptr.back()+num+1);
806      }      }
807    
808      dim_t numDOFs = 0;      return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
809      for (dim_t n = 0; n < len; n++) {  }
810          if (tmp_node_localDOF_mask[n] >= 0) {  
811              tmp_node_localDOF_mask[n] = numDOFs;  //protected
812              numDOFs++;  void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
813          }          const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
814    {
815        IndexVector ptr(1,0);
816        IndexVector index;
817        for (index_t i=0; i<getNumDOF(); i++) {
818            index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
819            ptr.push_back(ptr.back()+colIndices[i].size());
820      }      }
821  #pragma omp parallel for  
822      for (dim_t n = 0; n < m_nodes->getNumNodes(); n++) {      const dim_t M=ptr.size()-1;
823          tmp_node_localDOF_map[n] = tmp_node_localDOF_mask[m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first];      *colPattern=createPasoPattern(ptr, index, M, N);
824    
825        IndexVector rowPtr(1,0);
826        IndexVector rowIndex;
827        for (dim_t id=0; id<N; id++) {
828            dim_t n=0;
829            for (dim_t i=0; i<M; i++) {
830                for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
831                    if (index[j]==id) {
832                        rowIndex.push_back(i);
833                        n++;
834                        break;
835                    }
836                }
837            }
838            rowPtr.push_back(rowPtr.back()+n);
839      }      }
840      // create element coloring  
841      createColoring(tmp_node_localDOF_map);      // M and N are now reversed
842        *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
843  }  }
844    
845  /**************************************************************  //protected
846     Check whether there is any node which has no vertex. In case  void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
847     such a node exists, we don't use ParMetis since it requires         const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
848     that every node has at least 1 vertex (at line 129 of file         dim_t num_Sol, const vector<double>& array) const
849     "xyzpart.c" in parmetis 3.1.1, variable "nvtxs" would be 0 if  {
850     any node has no vertex).      const dim_t numMyCols = mat->pattern->mainPattern->numInput;
851   **************************************************************/      const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
852  #ifdef USE_PARMETIS  
853  static int Check_Inputs_For_Parmetis(dim_t mpiSize, dim_t rank,      const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
854                          const vector<dim_t> &distribution, MPI_Comm *comm)      const index_t* mainBlock_index = mat->mainBlock->pattern->index;
855  {      double* mainBlock_val = mat->mainBlock->val;
856      dim_t i, len;      const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
857      int ret_val = 1;      const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
858        double* col_coupleBlock_val = mat->col_coupleBlock->val;
859      if (rank == 0) {      const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
860          for (i = 0; i < mpiSize; i++) {      const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
861              len = distribution[i + 1] - distribution[i];      double* row_coupleBlock_val = mat->row_coupleBlock->val;
862              if (len == 0) {  
863                  ret_val = 0;      for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
864                  break;          // down columns of array
865            const dim_t j_Eq = nodes_Eq[k_Eq];
866            const dim_t i_row = j_Eq;
867    //printf("row:%d\n", i_row);
868            // only look at the matrix rows stored on this processor
869            if (i_row < numMyRows) {
870                for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
871                    const dim_t i_col = nodes_Sol[k_Sol];
872                    if (i_col < numMyCols) {
873                        for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
874                            if (mainBlock_index[k] == i_col) {
875                                mainBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
876                                break;
877                            }
878                        }
879                    } else {
880                        for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
881    //cout << "col:" << i_col-numMyCols << " colIdx:" << col_coupleBlock_index[k] << endl;
882                            if (col_coupleBlock_index[k] == i_col - numMyCols) {
883                                col_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
884                                break;
885                            }
886                        }
887                    }
888                }
889            } else {
890                for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
891                    // across rows of array
892                    const dim_t i_col = nodes_Sol[k_Sol];
893                    if (i_col < numMyCols) {
894                        for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
895                             k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
896                        {
897    //cout << "col:" << i_col << " rowIdx:" << row_coupleBlock_index[k] << endl;
898                            if (row_coupleBlock_index[k] == i_col) {
899                                row_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
900                                break;
901                            }
902                        }
903                    }
904              }              }
905          }          }
906      }      }
     MPI_Bcast(&ret_val, 1, MPI_INTEGER, 0, *comm);  
     if (ret_val == 0)  
         printf("INFO: Not using ParMetis since some nodes have no vertex!\n");  
     return ret_val;  
907  }  }
 #endif  
908    
909  void RipleyDomain::optimizeDOFDistribution(RankVector &distribution)  //
910    // the methods that follow have to be implemented by the subclasses
911    //
912    
913    string RipleyDomain::getDescription() const
914  {  {
915      const Esys_MPI_rank myRank = m_mpiInfo->rank;      throw RipleyException("getDescription() not implemented");
916      dim_t mpiSize = m_mpiInfo->size;  }
     dim_t dim = getDim();  
   
     // first step is to distribute the elements according to a global X of DOF  
     const index_t myFirstVertex = distribution[myRank];  
     const index_t myLastVertex = distribution[myRank + 1];  
     const dim_t myNumVertices = myLastVertex - myFirstVertex;  
     const dim_t globalNumVertices = distribution[mpiSize];  
     dim_t len = 0;  
     for (dim_t p = 0; p < mpiSize; ++p)  
         len = MAX(len, distribution[p+1] - distribution[p]);  
     index_t *partition = TMPMEMALLOC(len, index_t);  
     float *xyz = TMPMEMALLOC(myNumVertices * dim, float);  
     if (!(Esys_checkPtr(partition) || Esys_checkPtr(xyz))) {  
         // set the coordinates  
         // it is assumed that at least one node on this processor provides a  
         // coordinate  
 #pragma omp parallel for  
         for (dim_t i = 0; i < m_nodes->getNumNodes(); ++i) {  
             dim_t k = m_nodes->getGlobalDegreesOfFreedom()[i] - myFirstVertex;  
             if ((k >= 0) && (k < myNumVertices)) {  
                 for (dim_t j = 0; j < dim; ++j)  
                     xyz[k * dim + j] = (float)(m_nodes->getCoordinates()[INDEX2(j, i, dim)]);  
             }  
         }  
917    
918          IndexMatrix indexList(myNumVertices);  void RipleyDomain::write(const string& filename) const
919          // ksteube CSR of DOF IDs  {
920          // create the adjacency structure xadj and adjncy      throw RipleyException("write() not implemented");
921          // insert contributions from element matrices into columns of indexList  }
         m_elements->insertIntoIndexMatrixNoMainDiagonal(indexList,  
                 m_nodes->getGlobalDegreesOfFreedom(),  
                 m_nodes->getGlobalDegreesOfFreedom(),  
                 myFirstVertex, myLastVertex);  
         m_faceElements->insertIntoIndexMatrixNoMainDiagonal(indexList,  
                 m_nodes->getGlobalDegreesOfFreedom(),  
                 m_nodes->getGlobalDegreesOfFreedom(),  
                 myFirstVertex, myLastVertex);  
         m_points->insertIntoIndexMatrixNoMainDiagonal(indexList,  
                 m_nodes->getGlobalDegreesOfFreedom(),  
                 m_nodes->getGlobalDegreesOfFreedom(),  
                 myFirstVertex, myLastVertex);  
   
         // create the local matrix pattern  
         Paso_Pattern *pattern = createPatternFromIndexMatrix(  
                 indexList, 0, myNumVertices, 0, globalNumVertices, 0);  
   
 #ifdef USE_PARMETIS  
         if (mpiSize > 1 && Check_Inputs_For_Parmetis(mpiSize, myRank, distribution, &(m_mpiInfo->comm)) > 0) {  
             int wgtflag = 0;  
             int numflag = 0; /* pattern->ptr is C style: starting from 0 instead of 1 */  
             int ncon = 1;  
             int edgecut;  
             int options[2];  
             float *tpwgts = TMPMEMALLOC(ncon * mpiSize, float);  
             float *ubvec = TMPMEMALLOC(ncon, float);  
             for (dim_t i = 0; i < ncon * mpiSize; i++)  
                 tpwgts[i] = 1.0 / (float)mpiSize;  
             for (dim_t i = 0; i < ncon; i++)  
                 ubvec[i] = 1.05;  
             options[0] = 3;  
             options[1] = 15;  
             ParMETIS_V3_PartGeomKway(&distribution[0], pattern->ptr,  
                     pattern->index, NULL, NULL, &wgtflag, &numflag,  
                     &dim, xyz, &ncon, &mpiSize, tpwgts, ubvec, options,  
                     &edgecut, partition, // new CPU ownership of elements  
                     &m_mpiInfo->comm);  
             TMPMEMFREE(ubvec);  
             TMPMEMFREE(tpwgts);  
         } else {  
             for (dim_t i = 0; i < myNumVertices; ++i)  
                 partition[i] = 0; // CPU 0 owns it  
         }  
 #else  
         for (dim_t i = 0; i < myNumVertices; ++i)  
             partition[i] = myRank; // CPU myRank owns it  
 #endif  
922    
923          Paso_Pattern_free(pattern);  void RipleyDomain::dump(const string& filename) const
924    {
925        throw RipleyException("dump() not implemented");
926    }
927    
928          // create a new distribution and labeling of the DOF  const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
929          RankVector newDistribution(mpiSize+1, 0);  {
930  #pragma omp parallel      throw RipleyException("borrowSampleReferenceIDs() not implemented");
931          {  }
             RankVector loc_partition_count(mpiSize, 0);  
 #pragma omp for  
             for (dim_t i = 0; i < myNumVertices; ++i)  
                 loc_partition_count[partition[i]]++;  
 #pragma omp critical  
             {  
                 for (dim_t i = 0; i < mpiSize; ++i)  
                     newDistribution[i] += loc_partition_count[i];  
             }  
         }  
932    
933          // recvbuf will be the concatenation of each CPU's contribution  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
934          // to newDistribution  {
935          dim_t *recvbuf = TMPMEMALLOC(mpiSize * mpiSize, dim_t);      throw RipleyException("interpolateACross() not implemented");
936    }
937    
938  #ifdef ESYS_MPI  bool RipleyDomain::probeInterpolationACross(int fsType_source,
939          MPI_Allgather(&newDistribution[0], mpiSize, MPI_INT, recvbuf, mpiSize, MPI_INT, m_mpiInfo->comm);          const escript::AbstractDomain&, int fsType_target) const
940  #else  {
941          for (dim_t i = 0; i < mpiSize; ++i)      throw RipleyException("probeInterpolationACross() not implemented");
942              recvbuf[i] = newDistribution[i];  }
 #endif  
         newDistribution[0] = 0;  
         IndexVector newGlobalDOFid(len);  
         for (Esys_MPI_rank rank = 0; rank < mpiSize; rank++) {  
             int c = 0;  
             for (dim_t i = 0; i < myRank; ++i)  
                 c += recvbuf[rank + mpiSize * i];  
             for (dim_t i = 0; i < myNumVertices; ++i) {  
                 if (rank == partition[i]) {  
                     newGlobalDOFid[i] = newDistribution[rank] + c;  
                     c++;  
                 }  
             }  
             for (dim_t i = myRank + 1; i < mpiSize; ++i)  
                 c += recvbuf[rank + mpiSize * i];  
             newDistribution[rank + 1] = newDistribution[rank] + c;  
         }  
         TMPMEMFREE(recvbuf);  
943    
944          // now the overlap needs to be created by sending the partition  void RipleyDomain::setToNormal(escript::Data& normal) const
945          // around  {
946          m_nodes->resetGlobalDegreesOfFreedom(newGlobalDOFid, distribution);      throw RipleyException("setToNormal() not implemented");
947          for (dim_t i = 0; i < mpiSize+1; ++i)  }
             distribution[i] = newDistribution[i];  
     }  
     TMPMEMFREE(partition);  
     TMPMEMFREE(xyz);  
 }  
   
 void RipleyDomain::optimizeDOFLabeling(const IndexVector &distribution)  
 {  
     Esys_MPI_rank myRank = m_mpiInfo->rank;  
     index_t myFirstVertex = distribution[myRank];  
     index_t myLastVertex = distribution[myRank + 1];  
     dim_t myNumVertices = myLastVertex - myFirstVertex;  
     dim_t len = 0;  
     for (dim_t p = 0; p < m_mpiInfo->size; ++p)  
         len = MAX(len, distribution[p + 1] - distribution[p]);  
   
     IndexMatrix indexList(myNumVertices);  
   
     // create the adjacency structure xadj and adjncy  
     // insert contributions from element matrices into columns of indexList  
     m_elements->insertIntoIndexMatrixNoMainDiagonal(indexList,  
             m_nodes->getGlobalDegreesOfFreedom(),  
             m_nodes->getGlobalDegreesOfFreedom(),  
             myFirstVertex, myLastVertex);  
     m_faceElements->insertIntoIndexMatrixNoMainDiagonal(indexList,  
             m_nodes->getGlobalDegreesOfFreedom(),  
             m_nodes->getGlobalDegreesOfFreedom(),  
             myFirstVertex, myLastVertex);  
     m_points->insertIntoIndexMatrixNoMainDiagonal(indexList,  
             m_nodes->getGlobalDegreesOfFreedom(),  
             m_nodes->getGlobalDegreesOfFreedom(),  
             myFirstVertex, myLastVertex);  
   
     // create the local matrix pattern  
     Paso_Pattern *pattern = createPatternFromIndexMatrix(indexList, 0,  
             myNumVertices, myFirstVertex, myLastVertex, -myFirstVertex);  
   
     IndexVector newGlobalDOFid(len);  
     Paso_Pattern_reduceBandwidth(pattern, &newGlobalDOFid[0]);  
     Paso_Pattern_free(pattern);  
948    
949      Esys_MPIInfo_noError(m_mpiInfo);  void RipleyDomain::setToSize(escript::Data& size) const
950      checkPasoError();  {
951        throw RipleyException("setToSize() not implemented");
952    }
953    
954      /* shift new labeling to create a global id */  void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
955  #pragma omp parallel for  {
956      for (dim_t i = 0; i < myNumVertices; ++i)      throw RipleyException("setToGradient() not implemented");
957          newGlobalDOFid[i] += myFirstVertex;  }
958    
959      /* distribute new labeling to other processors */  void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
960      m_nodes->resetGlobalDegreesOfFreedom(newGlobalDOFid, distribution);  {
961  #if 0      throw RipleyException("setToIntegrals() not implemented");
     for (i = 0; i < m_nodes->getNumNodes(); ++i)  
         printf("%d ", m_nodes->getGlobalDegreesOfFreedom()[i]);  
     printf("\n");  
 #endif  
962  }  }
963    
964  void RipleyDomain::resolveNodeIds()  bool RipleyDomain::ownSample(int fsType, index_t id) const
965  {  {
966      // find the minimum and maximum node ID used by elements      throw RipleyException("ownSample() not implemented");
967      IndexPair minMaxId = m_elements->getNodeRange();  }
     index_t minId = minMaxId.first;  
     index_t maxId = minMaxId.second;  
     minMaxId = m_faceElements->getNodeRange();  
     minId = MIN(minId, minMaxId.first);  
     maxId = MAX(maxId, minMaxId.second);  
     minMaxId = m_points->getNodeRange();  
     minId = MIN(minId, minMaxId.first);  
     maxId = MAX(maxId, minMaxId.second);  
968    
969  #ifdef Ripley_TRACE  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
970      index_t globalMinId, globalMaxId;          const escript::Data& D, const escript::Data& d,
971  #ifdef ESYS_MPI          const escript::Data& d_dirac, const bool useHRZ) const
972      index_t idRange[2], globalIdRange[2];  {
973      idRange[0] = -minId;      throw RipleyException("addPDEToLumpedSystem() not implemented");
974      idRange[1] = maxId;  }
     MPI_Allreduce(idRange, globalIdRange, 2, MPI_INT, MPI_MAX, m_mpiInfo->comm);  
     globalMinId = -globalIdRange[0];  
     globalMaxId = globalIdRange[1];  
 #else  
     globalMinId = minId;  
     globalMaxId = maxId;  
 #endif // ESYS_MPI  
     cout << "Node id range used by elements is " << globalMinId << ":"  
         << globalMaxId << endl;  
 #endif // Ripley_TRACE  
   
     // this is only true if we have no elements at all  
     if (minId > maxId) {  
         maxId = -1;  
         minId = 0;  
     }  
   
     // allocate mappings for new local node labeling to global node labeling  
     // (newLocalToGlobalNodeLabels) and global node labeling to the new local  
     // node labeling. globalToNewLocalNodeLabels[i-minId] is the new local id  
     // of global node i  
     dim_t len = (maxId >= minId) ? maxId - minId + 1 : 0;  
     IndexVector globalToNewLocalNodeLabels(len, -1);  
   
     // mark the nodes referred by elements in globalToNewLocalNodeLabels  
     // which is currently used as a mask  
     markNodes(globalToNewLocalNodeLabels, minId);  
   
     // create a local labeling newLocalToGlobalNodeLabels of the local  
     // nodes by packing the mask globalToNewLocalNodeLabels  
     IndexVector newLocalToGlobalNodeLabels = packMask(globalToNewLocalNodeLabels);  
     const dim_t newNumNodes = newLocalToGlobalNodeLabels.size();  
   
     // invert the new labeling and shift the index newLocalToGlobalNodeLabels  
     // to global node IDs  
 #pragma omp parallel for schedule(static)  
     for (dim_t n = 0; n < newNumNodes; n++) {  
 #ifdef BOUNDS_CHECK  
         if (n >= len || n < 0) {  
             printf("BOUNDS_CHECK %s %d n=%d\n", __FILE__, __LINE__, n);  
             exit(1);  
         }  
         if (newLocalToGlobalNodeLabels[n] >= len || newLocalToGlobalNodeLabels[n] < 0) {  
             printf("BOUNDS_CHECK %s %d n=%d\n", __FILE__, __LINE__, n);  
             exit(1);  
         }  
 #endif  
         globalToNewLocalNodeLabels[newLocalToGlobalNodeLabels[n]] = n;  
         newLocalToGlobalNodeLabels[n] += minId;  
     }  
975    
976      // create a new node file  void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
977      m_nodes = m_nodes->gatherGlobal(newLocalToGlobalNodeLabels);          const escript::Data& Y, const escript::Data& y,
978            const escript::Data& y_contact, const escript::Data& y_dirac) const
979    {
980        throw RipleyException("addPDEToRHS() not implemented");
981    }
982    
983      // relabel nodes of the elements  void RipleyDomain::addPDEToTransportProblem(
984      relabelElementNodes(globalToNewLocalNodeLabels, minId);          escript::AbstractTransportProblem& tp,
985            escript::Data& source, const escript::Data& M,
986            const escript::Data& A, const escript::Data& B, const escript::Data& C,
987            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
988            const escript::Data& d, const escript::Data& y,
989            const escript::Data& d_contact, const escript::Data& y_contact,
990            const escript::Data& d_dirac, const escript::Data& y_dirac) const
991    {
992        throw RipleyException("addPDEToTransportProblem() not implemented");
993  }  }
994    
995  void RipleyDomain::relabelElementNodes(const IndexVector &newNode, index_t offset)  escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
996            const int blocksize, const escript::FunctionSpace& functionspace,
997            const int type) const
998  {  {
999      m_elements->relabelNodes(newNode, offset);      throw RipleyException("newTransportProblem() not implemented");
     m_faceElements->relabelNodes(newNode, offset);  
     m_points->relabelNodes(newNode, offset);  
1000  }  }
1001    
1002  void RipleyDomain::markNodes(IndexVector &mask, index_t offset)  Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
1003                                                       bool reducedColOrder) const
1004  {  {
1005      m_elements->markNodes(mask, offset);      throw RipleyException("getPattern() not implemented");
     m_faceElements->markNodes(mask, offset);  
     m_points->markNodes(mask, offset);  
1006  }  }
1007    
1008  void RipleyDomain::createMappings(const IndexVector &dofDistribution, const IndexVector &nodeDistribution)  dim_t RipleyDomain::getNumDataPointsGlobal() const
1009  {  {
1010      IndexVector maskReducedNodes(m_nodes->getNumNodes(), -1);      throw RipleyException("getNumDataPointsGlobal() not implemented");
1011      markNodes(maskReducedNodes, 0);  }
1012    
1013      IndexVector indexReducedNodes = packMask(maskReducedNodes);  IndexVector RipleyDomain::getNumNodesPerDim() const
1014      m_nodes->createNodeFileMappings(indexReducedNodes, dofDistribution,  {
1015                                      nodeDistribution);      throw RipleyException("getNumNodesPerDim() not implemented");
1016  }  }
1017    
1018  Paso_SystemMatrixPattern *RipleyDomain::makePattern(bool reduce_row_order, bool reduce_col_order) const  IndexVector RipleyDomain::getNumElementsPerDim() const
1019  {  {
1020      Paso_Connector *colConnector, *rowConnector;      throw RipleyException("getNumElementsPerDim() not implemented");
1021      NodeMapping_ptr colMap, rowMap;  }
     Paso_Distribution *colDistribution = NULL, *rowDistribution = NULL;  
1022    
1023      if (reduce_col_order) {  IndexVector RipleyDomain::getNumFacesPerBoundary() const
1024          colMap = m_nodes->getReducedDegreesOfFreedomMapping();  {
1025          colDistribution = m_nodes->getReducedDegreesOfFreedomDistribution();      throw RipleyException("getNumFacesPerBoundary() not implemented");
1026          colConnector = m_nodes->getReducedDegreesOfFreedomConnector();  }
     } else {  
         colMap = m_nodes->getDegreesOfFreedomMapping();  
         colDistribution = m_nodes->getDegreesOfFreedomDistribution();  
         colConnector = m_nodes->getDegreesOfFreedomConnector();  
     }  
1027    
1028      if (reduce_row_order) {  IndexVector RipleyDomain::getNodeDistribution() const
1029          rowMap = m_nodes->getReducedDegreesOfFreedomMapping();  {
1030          rowDistribution = m_nodes->getReducedDegreesOfFreedomDistribution();      throw RipleyException("getNodeDistribution() not implemented");
1031          rowConnector = m_nodes->getReducedDegreesOfFreedomConnector();  }
1032      } else {  
1033          rowMap = m_nodes->getDegreesOfFreedomMapping();  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1034          rowDistribution = m_nodes->getDegreesOfFreedomDistribution();  {
1035          rowConnector = m_nodes->getDegreesOfFreedomConnector();      throw RipleyException("getFirstCoordAndSpacing() not implemented");
1036      }  }
1037    
1038      //  insert contributions from element matrices into columns in indexList  dim_t RipleyDomain::getNumFaceElements() const
1039      IndexMatrix indexList(rowMap->numTargets);  {
1040      m_elements->insertIntoIndexMatrix(indexList, rowMap->target, colMap->target);      throw RipleyException("getNumFaceElements() not implemented");
1041      m_faceElements->insertIntoIndexMatrix(indexList, rowMap->target, colMap->target);  }
1042      m_points->insertIntoIndexMatrix(indexList, rowMap->target, colMap->target);  
1043    dim_t RipleyDomain::getNumElements() const
1044      // create patterns  {
1045      Paso_Pattern *mainPattern = createPatternFromIndexMatrix(indexList, 0,      throw RipleyException("getNumElements() not implemented");
1046              Paso_Distribution_getMyNumComponents(rowDistribution),  }
1047              0, Paso_Distribution_getMyNumComponents(colDistribution), 0);  
1048      Paso_Pattern *colCouplePattern = createPatternFromIndexMatrix(indexList, 0,  dim_t RipleyDomain::getNumNodes() const
1049              Paso_Distribution_getMyNumComponents(rowDistribution),  {
1050              Paso_Distribution_getMyNumComponents(colDistribution),      throw RipleyException("getNumNodes() not implemented");
1051              colMap->numTargets,  }
1052              -Paso_Distribution_getMyNumComponents(colDistribution));  
1053      Paso_Pattern *rowCouplePattern = createPatternFromIndexMatrix(indexList,  dim_t RipleyDomain::getNumDOF() const
1054              Paso_Distribution_getMyNumComponents(rowDistribution),  {
1055              rowMap->numTargets, 0,      throw RipleyException("getNumDOF() not implemented");
1056              Paso_Distribution_getMyNumComponents(colDistribution), 0);  }
1057    
1058      // if everything is in order we can create the return value  dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1059      Paso_SystemMatrixPattern *out = Paso_SystemMatrixPattern_alloc(  {
1060              MATRIX_FORMAT_DEFAULT, rowDistribution, colDistribution,      throw RipleyException("insertNeighbourNodes() not implemented");
1061              mainPattern, colCouplePattern, rowCouplePattern,  }
1062              colConnector, rowConnector);  
1063    void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1064      // clean up  {
1065      Paso_Pattern_free(mainPattern);      throw RipleyException("assembleCoordinates() not implemented");
1066      Paso_Pattern_free(colCouplePattern);  }
1067      Paso_Pattern_free(rowCouplePattern);  
1068      Esys_MPIInfo_noError(m_mpiInfo);  void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1069      return out;          const escript::Data& A, const escript::Data& B, const escript::Data& C,
1070  }          const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1071            const escript::Data& d, const escript::Data& y) const
1072  Paso_SystemMatrixPattern *RipleyDomain::getPattern(bool reduce_row_order, bool reduce_col_order) const  {
1073  {      throw RipleyException("assemblePDESingle() not implemented");
1074      return Paso_SystemMatrixPattern_getReference(makePattern(reduce_row_order, reduce_col_order));  }
1075      /*  
1076      Paso_SystemMatrixPattern *out = NULL;  void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1077            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1078      // make sure that the requested pattern is available          const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1079      if (reduce_row_order) {          const escript::Data& d, const escript::Data& y) const
1080          if (reduce_col_order) {  {
1081              if (m_reducedReducedPattern == NULL)      throw RipleyException("assemblePDESystem() not implemented");
1082                  m_reducedReducedPattern = makePattern(reduce_row_order, reduce_col_order);  }
1083          } else {  
1084              if (m_reducedFullPattern == NULL)  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1085                  m_reducedFullPattern = makePattern(reduce_row_order, reduce_col_order);  {
1086          }      throw RipleyException("interpolateNodesOnElements() not implemented");
1087      } else {  }
1088          if (reduce_col_order) {  
1089              if (m_fullReducedPattern == NULL)  void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
1090                  m_fullReducedPattern = makePattern(reduce_row_order, reduce_col_order);  {
1091          } else {      throw RipleyException("interpolateNodesOnFaces() not implemented");
1092              if (m_fullFullPattern == NULL)  }
1093                  m_fullFullPattern = makePattern(reduce_row_order, reduce_col_order);  
1094          }  void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1095      }  {
1096      if (Ripley_noError()) {      throw RipleyException("nodesToDOF() not implemented");
1097          if (reduce_row_order) {  }
1098              if (reduce_col_order) {  
1099                  out = Paso_SystemMatrixPattern_getReference(m_reducedReducedPattern);  void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1100              } else {  {
1101                  out = Paso_SystemMatrixPattern_getReference(m_reducedFullPattern);      throw RipleyException("dofToNodes() not implemented");
             }  
         } else {  
             if (reduce_col_order) {  
                 out = Paso_SystemMatrixPattern_getReference(m_fullReducedPattern);  
             } else {  
                 out = Paso_SystemMatrixPattern_getReference(m_fullFullPattern);  
             }  
         }  
     }  
     return out;  
     */  
1102  }  }
1103    
1104  } // end of namespace ripley  } // end of namespace ripley

Legend:
Removed from v.3670  
changed lines
  Added in v.3756

  ViewVC Help
Powered by ViewVC 1.1.26