/[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 3755 by caltinay, Thu Jan 5 06:51:31 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  
 {  
     return m_mpiInfo->size;  
 }  
   
 int RipleyDomain::getMPIRank() const  
29  {  {
30      return m_mpiInfo->rank;      throw RipleyException("loadMesh() not implemented");
31  }  }
32    
33  void RipleyDomain::MPIBarrier() const  escript::Domain_ptr RipleyDomain::readMesh(const string& filename)
34  {  {
35  #ifdef ESYS_MPI      throw RipleyException("readMesh() not implemented");
     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);  
57      }      }
58  #endif      return false;
   
     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++;  
         }  
     }  
   
     // 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);  
             break;  
         case FaceElements:  
             Ripley_Assemble_getSize(m_nodes, mesh->FaceElements, &tmp);  
119              break;              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);  
120      }      }
 }  
121    
122  bool RipleyDomain::ownSample(int fs_code, index_t id) const      stringstream msg;
123  {      msg << "getDataShape(): Unsupported function space type " << fsType
124  #ifdef ESYS_MPI          << " for " << getDescription();
125      index_t myFirstNode, myLastNode, k;      throw RipleyException(msg.str());
     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);  
     }  
     return static_cast<bool>((myFirstNode <= k) && (k < myLastNode));  
 #endif  
     return true;  
126  }  }
127    
128    string RipleyDomain::showTagNames() const
 //  
 // 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");  
   
     // generate matrix  
     Paso_SystemMatrixPattern* fsystemMatrixPattern = getPattern(reduceOrder, reduceOrder);  
     Paso_TransportProblem* transportProblem = Paso_TransportProblem_alloc(  
             useBackwardEuler, fsystemMatrixPattern, blocksize);  
     checkPasoError();  
     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);  
     escript::AbstractTransportProblem* atp=new TransportProblemAdapter(  
             transportProblem, useBackwardEuler, blocksize, functionspace);  
     return escript::ATP_ptr(atp);  
 }  
   
 // returns true if data on the function space is considered as being cell  
 // 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  void RipleyDomain::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const      switch(fsType) {
 {  
     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        // generate matrix
639        Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
640        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    void RipleyDomain::addPDEToSystem(
650            escript::AbstractSystemMatrix& mat, escript::Data& rhs,
651            const escript::Data& A, const escript::Data& B, const escript::Data& C,
652            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
653            const escript::Data& d, const escript::Data& y,
654            const escript::Data& d_contact, const escript::Data& y_contact,
655            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        if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
665            throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
666    
667        //TODO: more input param checks (shape, function space etc)
668    
669        Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
670    
671        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::setNewX(const escript::Data& arg)
688    {
689        throw RipleyException("setNewX(): Operation not supported");
690    }
691    
692    //protected
693    void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
694    {
695        const dim_t numComp = in.getDataPointSize();
696        const dim_t dpp = in.getNumDataPointsPerSample();
697        out.requireWrite();
698    #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      // Rearrange elements in order to try and bring them closer to memory  //protected
712      // locations of the nodes (distributed shared memory).  void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
713      optimizeElementOrdering();  {
714        const dim_t numComp = in.getDataPointSize();
715        out.requireWrite();
716    #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  /* useful DEBUG:  //protected
724      {  void RipleyDomain::updateTagsInUse(int fsType) const
725          printf("prepare: global DOF : %d\n",newGlobalNumDOFs);  {
726          IndexPair range = m_nodes->getGlobalIdRange();      IndexVector* tagsInUse=NULL;
727          printf("prepare: global node id range = %d :%d\n",range.first,range.second);      const IndexVector* tags=NULL;
728          range = m_nodes->getIdRange();      switch(fsType) {
729          printf("prepare: local node id range = %d :%d\n",range.first,range.second);          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      // create the global indices      // gather global unique tag values from tags into tagsInUse
748      IndexVector maskReducedNodes(m_nodes->getNumNodes(), -1);      tagsInUse->clear();
749      markNodes(maskReducedNodes, 0);      index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
     IndexVector indexReducedNodes = packMask(maskReducedNodes);  
750    
751      IndexVector nodeDistribution(m_mpiInfo->size+1);      while (true) {
752      m_nodes->createDenseNodeLabeling(nodeDistribution, distribution);          // find smallest value bigger than lastFoundValue
753      m_nodes->createDenseReducedDOFLabeling(maskReducedNodes);          minFoundValue = INDEX_T_MAX;
754      m_nodes->createDenseReducedNodeLabeling(maskReducedNodes);  #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 the missing mappings          // if we found a new value add it to the tagsInUse vector
775      m_nodes->createNodeFileMappings(indexReducedNodes, distribution,          if (minFoundValue < INDEX_T_MAX) {
776                                      nodeDistribution);              tagsInUse->push_back(minFoundValue);
777                lastFoundValue = minFoundValue;
778            } else
779                break;
780        }
781    }
782    
783      updateTagsInUse();  //
784    // the methods that follow have to be implemented by the subclasses
785    //
786    
787    string RipleyDomain::getDescription() const
788    {
789        throw RipleyException("getDescription() not implemented");
790  }  }
791    
792  void RipleyDomain::optimizeElementOrdering()  void RipleyDomain::write(const string& filename) const
793  {  {
794      m_elements->optimizeOrdering();      throw RipleyException("write() not implemented");
     m_faceElements->optimizeOrdering();  
     m_points->optimizeOrdering();  
795  }  }
796    
797  void RipleyDomain::updateTagsInUse()  void RipleyDomain::dump(const string& filename) const
798  {  {
799      m_nodes->updateTagsInUse();      throw RipleyException("dump() not implemented");
     m_elements->updateTagsInUse();  
     m_faceElements->updateTagsInUse();  
     m_points->updateTagsInUse();  
800  }  }
801    
802  void RipleyDomain::createColoring(const IndexVector &node_localDOF_map)  const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
803  {  {
804      m_elements->createColoring(node_localDOF_map);      throw RipleyException("borrowSampleReferenceIDs() not implemented");
     m_faceElements->createColoring(node_localDOF_map);  
     m_points->createColoring(node_localDOF_map);  
805  }  }
806    
807  void RipleyDomain::distributeByRankOfDOF(const IndexVector &dofDistribution)  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
808  {  {
809      RankVector mpiRankOfDOF = m_nodes->getOwnerOfDOFs(dofDistribution);      throw RipleyException("interpolateACross() not implemented");
810    }
811    
812      // First the elements are redistributed according to mpiRankOfDOF.  bool RipleyDomain::probeInterpolationACross(int fsType_source,
813      // At the input the Node tables refer to the local labeling of          const escript::AbstractDomain&, int fsType_target) const
814      // the nodes while at the output they refer to the global labeling  {
815      // which is rectified in the next step      throw RipleyException("probeInterpolationACross() not implemented");
816      m_elements->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());  }
     m_faceElements->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());  
     m_points->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());  
817    
818      // resolve the node IDs  void RipleyDomain::setToNormal(escript::Data& normal) const
819      resolveNodeIds();  {
820        throw RipleyException("setToNormal() not implemented");
821    }
822    
823      // create a local labeling of the DOFs  void RipleyDomain::setToSize(escript::Data& size) const
824      IndexPair dofRange = m_nodes->getDOFRange();  {
825      dim_t len = dofRange.second - dofRange.first + 1;      throw RipleyException("setToSize() not implemented");
826    }
827    
828      // local mask for used nodes  void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
829      IndexVector tmp_node_localDOF_mask(len, -1);  {
830      IndexVector tmp_node_localDOF_map(m_nodes->getNumNodes(), -1);      throw RipleyException("setToGradient() not implemented");
831    }
832    
833  #pragma omp parallel for schedule(static)  void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
834      for (dim_t n = 0; n < m_nodes->getNumNodes(); n++) {  {
835  #ifdef BOUNDS_CHECK      throw RipleyException("setToIntegrals() not implemented");
836          if ((m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first) >= len  }
                 || (m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first) < 0) {  
             printf("BOUNDS_CHECK %s %d\n", __FILE__, __LINE__);  
             exit(1);  
         }  
 #endif  
         tmp_node_localDOF_mask[m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first] = n;  
     }  
837    
838      dim_t numDOFs = 0;  bool RipleyDomain::ownSample(int fsType, index_t id) const
839      for (dim_t n = 0; n < len; n++) {  {
840          if (tmp_node_localDOF_mask[n] >= 0) {      throw RipleyException("ownSample() not implemented");
             tmp_node_localDOF_mask[n] = numDOFs;  
             numDOFs++;  
         }  
     }  
 #pragma omp parallel for  
     for (dim_t n = 0; n < m_nodes->getNumNodes(); n++) {  
         tmp_node_localDOF_map[n] = tmp_node_localDOF_mask[m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first];  
     }  
     // create element coloring  
     createColoring(tmp_node_localDOF_map);  
841  }  }
842    
843  /**************************************************************  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
844     Check whether there is any node which has no vertex. In case          const escript::Data& D, const escript::Data& d,
845     such a node exists, we don't use ParMetis since it requires          const escript::Data& d_dirac, const bool useHRZ) const
846     that every node has at least 1 vertex (at line 129 of file  {
847     "xyzpart.c" in parmetis 3.1.1, variable "nvtxs" would be 0 if      throw RipleyException("addPDEToLumpedSystem() not implemented");
    any node has no vertex).  
  **************************************************************/  
 #ifdef USE_PARMETIS  
 static int Check_Inputs_For_Parmetis(dim_t mpiSize, dim_t rank,  
                         const vector<dim_t> &distribution, MPI_Comm *comm)  
 {  
     dim_t i, len;  
     int ret_val = 1;  
   
     if (rank == 0) {  
         for (i = 0; i < mpiSize; i++) {  
             len = distribution[i + 1] - distribution[i];  
             if (len == 0) {  
                 ret_val = 0;  
                 break;  
             }  
         }  
     }  
     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;  
848  }  }
 #endif  
849    
850  void RipleyDomain::optimizeDOFDistribution(RankVector &distribution)  void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
851            const escript::Data& Y, const escript::Data& y,
852            const escript::Data& y_contact, const escript::Data& y_dirac) const
853  {  {
854      const Esys_MPI_rank myRank = m_mpiInfo->rank;      throw RipleyException("addPDEToRHS() not implemented");
855      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)]);  
             }  
         }  
856    
857          IndexMatrix indexList(myNumVertices);  void RipleyDomain::addPDEToTransportProblem(
858          // ksteube CSR of DOF IDs          escript::AbstractTransportProblem& tp,
859          // create the adjacency structure xadj and adjncy          escript::Data& source, const escript::Data& M,
860          // insert contributions from element matrices into columns of indexList          const escript::Data& A, const escript::Data& B, const escript::Data& C,
861          m_elements->insertIntoIndexMatrixNoMainDiagonal(indexList,          const escript::Data& D, const escript::Data& X, const escript::Data& Y,
862                  m_nodes->getGlobalDegreesOfFreedom(),          const escript::Data& d, const escript::Data& y,
863                  m_nodes->getGlobalDegreesOfFreedom(),          const escript::Data& d_contact, const escript::Data& y_contact,
864                  myFirstVertex, myLastVertex);          const escript::Data& d_dirac, const escript::Data& y_dirac) const
865          m_faceElements->insertIntoIndexMatrixNoMainDiagonal(indexList,  {
866                  m_nodes->getGlobalDegreesOfFreedom(),      throw RipleyException("addPDEToTransportProblem() not implemented");
867                  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  
868    
869          Paso_Pattern_free(pattern);  escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
870            const int blocksize, const escript::FunctionSpace& functionspace,
871            const int type) const
872    {
873        throw RipleyException("newTransportProblem() not implemented");
874    }
875    
876          // create a new distribution and labeling of the DOF  Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
877          RankVector newDistribution(mpiSize+1, 0);                                                     bool reducedColOrder) const
878  #pragma omp parallel  {
879          {      throw RipleyException("getPattern() not implemented");
880              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];  
             }  
         }  
881    
882          // recvbuf will be the concatenation of each CPU's contribution  dim_t RipleyDomain::getNumDataPointsGlobal() const
883          // to newDistribution  {
884          dim_t *recvbuf = TMPMEMALLOC(mpiSize * mpiSize, dim_t);      throw RipleyException("getNumDataPointsGlobal() not implemented");
885    }
886    
887  #ifdef ESYS_MPI  IndexVector RipleyDomain::getNumNodesPerDim() const
888          MPI_Allgather(&newDistribution[0], mpiSize, MPI_INT, recvbuf, mpiSize, MPI_INT, m_mpiInfo->comm);  {
889  #else      throw RipleyException("getNumNodesPerDim() not implemented");
890          for (dim_t i = 0; i < mpiSize; ++i)  }
             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);  
891    
892          // now the overlap needs to be created by sending the partition  IndexVector RipleyDomain::getNumElementsPerDim() const
893          // around  {
894          m_nodes->resetGlobalDegreesOfFreedom(newGlobalDOFid, distribution);      throw RipleyException("getNumElementsPerDim() not implemented");
895          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);  
896    
897      Esys_MPIInfo_noError(m_mpiInfo);  IndexVector RipleyDomain::getNumFacesPerBoundary() const
898      checkPasoError();  {
899        throw RipleyException("getNumFacesPerBoundary() not implemented");
900    }
901    
902      /* shift new labeling to create a global id */  IndexVector RipleyDomain::getNodeDistribution() const
903  #pragma omp parallel for  {
904      for (dim_t i = 0; i < myNumVertices; ++i)      throw RipleyException("getNodeDistribution() not implemented");
905          newGlobalDOFid[i] += myFirstVertex;  }
906    
907      /* distribute new labeling to other processors */  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
908      m_nodes->resetGlobalDegreesOfFreedom(newGlobalDOFid, distribution);  {
909  #if 0      throw RipleyException("getFirstCoordAndSpacing() not implemented");
     for (i = 0; i < m_nodes->getNumNodes(); ++i)  
         printf("%d ", m_nodes->getGlobalDegreesOfFreedom()[i]);  
     printf("\n");  
 #endif  
910  }  }
911    
912  void RipleyDomain::resolveNodeIds()  dim_t RipleyDomain::getNumFaceElements() const
913  {  {
914      // find the minimum and maximum node ID used by elements      throw RipleyException("getNumFaceElements() not implemented");
915      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);  
916    
917  #ifdef Ripley_TRACE  dim_t RipleyDomain::getNumElements() const
918      index_t globalMinId, globalMaxId;  {
919  #ifdef ESYS_MPI      throw RipleyException("getNumElements() not implemented");
920      index_t idRange[2], globalIdRange[2];  }
     idRange[0] = -minId;  
     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;  
     }  
921    
922      // create a new node file  dim_t RipleyDomain::getNumNodes() const
923      m_nodes = m_nodes->gatherGlobal(newLocalToGlobalNodeLabels);  {
924        throw RipleyException("getNumNodes() not implemented");
925    }
926    
927      // relabel nodes of the elements  dim_t RipleyDomain::getNumDOF() const
928      relabelElementNodes(globalToNewLocalNodeLabels, minId);  {
929        throw RipleyException("getNumDOF() not implemented");
930  }  }
931    
932  void RipleyDomain::relabelElementNodes(const IndexVector &newNode, index_t offset)  void RipleyDomain::assembleCoordinates(escript::Data& arg) const
933  {  {
934      m_elements->relabelNodes(newNode, offset);      throw RipleyException("assembleCoordinates() not implemented");
     m_faceElements->relabelNodes(newNode, offset);  
     m_points->relabelNodes(newNode, offset);  
935  }  }
936    
937  void RipleyDomain::markNodes(IndexVector &mask, index_t offset)  void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
938            const escript::Data& A, const escript::Data& B, const escript::Data& C,
939            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
940            const escript::Data& d, const escript::Data& y) const
941  {  {
942      m_elements->markNodes(mask, offset);      throw RipleyException("assemblePDESingle() not implemented");
     m_faceElements->markNodes(mask, offset);  
     m_points->markNodes(mask, offset);  
943  }  }
944    
945  void RipleyDomain::createMappings(const IndexVector &dofDistribution, const IndexVector &nodeDistribution)  void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
946            const escript::Data& A, const escript::Data& B, const escript::Data& C,
947            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
948            const escript::Data& d, const escript::Data& y) const
949  {  {
950      IndexVector maskReducedNodes(m_nodes->getNumNodes(), -1);      throw RipleyException("assemblePDESystem() not implemented");
951      markNodes(maskReducedNodes, 0);  }
952    
953      IndexVector indexReducedNodes = packMask(maskReducedNodes);  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
954      m_nodes->createNodeFileMappings(indexReducedNodes, dofDistribution,  {
955                                      nodeDistribution);      throw RipleyException("interpolateNodesOnElements() not implemented");
956  }  }
957    
958  Paso_SystemMatrixPattern *RipleyDomain::makePattern(bool reduce_row_order, bool reduce_col_order) const  void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
959  {  {
960      Paso_Connector *colConnector, *rowConnector;      throw RipleyException("interpolateNodesOnFaces() not implemented");
961      NodeMapping_ptr colMap, rowMap;  }
     Paso_Distribution *colDistribution = NULL, *rowDistribution = NULL;  
962    
963      if (reduce_col_order) {  void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
964          colMap = m_nodes->getReducedDegreesOfFreedomMapping();  {
965          colDistribution = m_nodes->getReducedDegreesOfFreedomDistribution();      throw RipleyException("nodesToDOF() not implemented");
966          colConnector = m_nodes->getReducedDegreesOfFreedomConnector();  }
     } else {  
         colMap = m_nodes->getDegreesOfFreedomMapping();  
         colDistribution = m_nodes->getDegreesOfFreedomDistribution();  
         colConnector = m_nodes->getDegreesOfFreedomConnector();  
     }  
967    
968      if (reduce_row_order) {  void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
969          rowMap = m_nodes->getReducedDegreesOfFreedomMapping();  {
970          rowDistribution = m_nodes->getReducedDegreesOfFreedomDistribution();      throw RipleyException("dofToNodes() not implemented");
         rowConnector = m_nodes->getReducedDegreesOfFreedomConnector();  
     } else {  
         rowMap = m_nodes->getDegreesOfFreedomMapping();  
         rowDistribution = m_nodes->getDegreesOfFreedomDistribution();  
         rowConnector = m_nodes->getDegreesOfFreedomConnector();  
     }  
   
     //  insert contributions from element matrices into columns in indexList  
     IndexMatrix indexList(rowMap->numTargets);  
     m_elements->insertIntoIndexMatrix(indexList, rowMap->target, colMap->target);  
     m_faceElements->insertIntoIndexMatrix(indexList, rowMap->target, colMap->target);  
     m_points->insertIntoIndexMatrix(indexList, rowMap->target, colMap->target);  
   
     // create patterns  
     Paso_Pattern *mainPattern = createPatternFromIndexMatrix(indexList, 0,  
             Paso_Distribution_getMyNumComponents(rowDistribution),  
             0, Paso_Distribution_getMyNumComponents(colDistribution), 0);  
     Paso_Pattern *colCouplePattern = createPatternFromIndexMatrix(indexList, 0,  
             Paso_Distribution_getMyNumComponents(rowDistribution),  
             Paso_Distribution_getMyNumComponents(colDistribution),  
             colMap->numTargets,  
             -Paso_Distribution_getMyNumComponents(colDistribution));  
     Paso_Pattern *rowCouplePattern = createPatternFromIndexMatrix(indexList,  
             Paso_Distribution_getMyNumComponents(rowDistribution),  
             rowMap->numTargets, 0,  
             Paso_Distribution_getMyNumComponents(colDistribution), 0);  
   
     // if everything is in order we can create the return value  
     Paso_SystemMatrixPattern *out = Paso_SystemMatrixPattern_alloc(  
             MATRIX_FORMAT_DEFAULT, rowDistribution, colDistribution,  
             mainPattern, colCouplePattern, rowCouplePattern,  
             colConnector, rowConnector);  
   
     // clean up  
     Paso_Pattern_free(mainPattern);  
     Paso_Pattern_free(colCouplePattern);  
     Paso_Pattern_free(rowCouplePattern);  
     Esys_MPIInfo_noError(m_mpiInfo);  
     return out;  
 }  
   
 Paso_SystemMatrixPattern *RipleyDomain::getPattern(bool reduce_row_order, bool reduce_col_order) const  
 {  
     return Paso_SystemMatrixPattern_getReference(makePattern(reduce_row_order, reduce_col_order));  
     /*  
     Paso_SystemMatrixPattern *out = NULL;  
   
     // make sure that the requested pattern is available  
     if (reduce_row_order) {  
         if (reduce_col_order) {  
             if (m_reducedReducedPattern == NULL)  
                 m_reducedReducedPattern = makePattern(reduce_row_order, reduce_col_order);  
         } else {  
             if (m_reducedFullPattern == NULL)  
                 m_reducedFullPattern = makePattern(reduce_row_order, reduce_col_order);  
         }  
     } else {  
         if (reduce_col_order) {  
             if (m_fullReducedPattern == NULL)  
                 m_fullReducedPattern = makePattern(reduce_row_order, reduce_col_order);  
         } else {  
             if (m_fullFullPattern == NULL)  
                 m_fullFullPattern = makePattern(reduce_row_order, reduce_col_order);  
         }  
     }  
     if (Ripley_noError()) {  
         if (reduce_row_order) {  
             if (reduce_col_order) {  
                 out = Paso_SystemMatrixPattern_getReference(m_reducedReducedPattern);  
             } else {  
                 out = Paso_SystemMatrixPattern_getReference(m_reducedFullPattern);  
             }  
         } else {  
             if (reduce_col_order) {  
                 out = Paso_SystemMatrixPattern_getReference(m_fullReducedPattern);  
             } else {  
                 out = Paso_SystemMatrixPattern_getReference(m_fullFullPattern);  
             }  
         }  
     }  
     return out;  
     */  
971  }  }
972    
973  } // end of namespace ripley  } // end of namespace ripley

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

  ViewVC Help
Powered by ViewVC 1.1.26