/[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 3690 by caltinay, Wed Nov 16 02:21:18 2011 UTC revision 3691 by caltinay, Wed Nov 23 23:07:37 2011 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  
 {  
     return m_mpiInfo->rank;  
 }  
   
 void RipleyDomain::MPIBarrier() const  
 {  
 #ifdef ESYS_MPI  
     MPI_Barrier(m_mpiInfo->comm);  
 #endif  
 }  
   
 bool RipleyDomain::onMasterProcessor() const  
 {  
     return m_mpiInfo->rank == 0;  
 }  
   
 #ifdef ESYS_MPI  
 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  
 {  
     cout << "Ripley_PrintMesh_Info running on CPU " <<  
         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;  
         }  
     }  
 }  
   
 void RipleyDomain::dump(const string& fileName) const  
 {  
 #ifdef USE_NETCDF  
     NcDim* ncdim = NULL;  
     int mpi_size = m_mpiInfo->size;  
     int mpi_rank = m_mpiInfo->rank;  
     int numDim = getDim();  
   
 /* Incoming token indicates it's my turn to write */  
 #ifdef ESYS_MPI  
     MPI_Status status;  
     if (mpi_rank>0) {  
         int dummy;  
         MPI_Recv(&dummy, 0, MPI_INT, mpi_rank-1, 81800, m_mpiInfo->comm, &status);  
     }  
 #endif  
   
     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),  
                                                       mpi_size, mpi_rank);  
   
     // NetCDF error handler  
     NcError err(NcError::verbose_nonfatal);  
     // Create the file.  
     NcFile dataFile(newFileName, NcFile::Replace);  
     const string msgPrefix("dump: netCDF operation failed - ");  
   
     // check if writing was successful  
     if (!dataFile.is_valid())  
         throw RipleyException(msgPrefix+"Open file for output");  
   
     const size_t numTags = m_tagMap.size();  
   
     if (numTags > 0)  
         if (! (ncdim = dataFile.add_dim("dim_Tags", numTags)) )  
             throw RipleyException(msgPrefix+"add_dim(dim_Tags)");  
   
     // Attributes: MPI size, MPI rank, Name, order, reduced_order  
     if (!dataFile.add_att("mpi_size", mpi_size))  
         throw RipleyException(msgPrefix+"add_att(mpi_size)");  
     if (!dataFile.add_att("mpi_rank", mpi_rank))  
         throw RipleyException(msgPrefix+"add_att(mpi_rank)");  
     if (!dataFile.add_att("Name", m_name.c_str()))  
         throw RipleyException(msgPrefix+"add_att(Name)");  
     if (!dataFile.add_att("numDim", numDim))  
         throw RipleyException(msgPrefix+"add_att(numDim)");  
     if (!dataFile.add_att("order", 1))  
         throw RipleyException(msgPrefix+"add_att(order)");  
     if (!dataFile.add_att("reduced_order", 1))  
         throw RipleyException(msgPrefix+"add_att(reduced_order)");  
     if (!dataFile.add_att("num_Tags", static_cast<int>(numTags)))  
         throw RipleyException(msgPrefix+"add_att(num_Tags)");  
   
     m_nodes->dumpToNetCDF(dataFile);  
     m_elements->dumpToNetCDF(dataFile, "Elements");  
     m_faceElements->dumpToNetCDF(dataFile, "FaceElements");  
     m_points->dumpToNetCDF(dataFile, "Points");  
   
     // // // // // TagMap // // // // //  
     if (numTags > 0) {  
   
         // Copy tag keys into temp array  
         vector<int> Tags_keys;  
         TagMap::const_iterator it;  
         for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {  
             Tags_keys.push_back(it->second);  
         }  
   
         // Tags_keys  
         NcVar *ncVar;  
         if (! (ncVar = dataFile.add_var("Tags_keys", ncInt, ncdim)) )  
             throw RipleyException(msgPrefix+"add_var(Tags_keys)");  
         if (!ncVar->put(&Tags_keys[0], numTags))  
             throw RipleyException(msgPrefix+"put(Tags_keys)");  
   
         // Tags_names_*  
         // This is an array of strings, it should be stored as an array but  
         // instead I have hacked in one attribute per string because the NetCDF  
         // manual doesn't tell how to do an array of strings  
         int i = 0;  
         for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {  
             stringstream name;  
             name << "Tags_name_" << i;  
             if (!dataFile.add_att(name.str().c_str(), it->first.c_str()) )  
                 throw RipleyException(msgPrefix+"add_att(Tags_names_XX)");  
             i++;  
         }  
     }  
   
     // 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  
29  {  {
30      FunctionSpaceNamesMapType::iterator loc;      throw RipleyException("loadMesh() not implemented");
     loc=m_functionSpaceTypeNames.find(functionSpaceType);  
     return (loc!=m_functionSpaceTypeNames.end());  
31  }  }
32    
33  void RipleyDomain::setFunctionSpaceTypeNames()  escript::Domain_ptr RipleyDomain::readMesh(const string& filename)
34  {  {
35      if (m_functionSpaceTypeNames.size() > 0)      throw RipleyException("readMesh() not implemented");
         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"));  
36  }  }
37    
38  int RipleyDomain::getContinuousFunctionCode() const  RipleyDomain::RipleyDomain(dim_t dim) :
39        m_numDim(dim),
40        m_status(0)
41  {  {
42      return Nodes;      m_mpiInfo = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
43  }  }
44    
45  int RipleyDomain::getReducedContinuousFunctionCode() const  RipleyDomain::~RipleyDomain()
 {  
     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  
46  {  {
47      return m_nodes->getGlobalNumNodes();      Esys_MPIInfo_free(m_mpiInfo);
48  }  }
49    
50  //  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  
51  {  {
52      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;  
53          case DegreesOfFreedom:          case DegreesOfFreedom:
             numDataPointsPerSample=1;  
             numSamples=m_nodes->getNumDegreesOfFreedom();  
             break;  
54          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()) {  
55          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;  
56          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;  
57          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;  
58          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;  
59          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;  
60          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;  
61          case Points:          case Points:
62              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;  
63          default:          default:
64              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);  
65      }      }
66        return false;
67  }  }
68    
69  //  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  
70  {  {
71  /*      switch (fsType) {
72      const RipleyDomain& normalDomain=dynamic_cast<const RipleyDomain&>(*(normal.getFunctionSpace().getDomain()));          case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
73      if (normalDomain!=*this)          case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
74          throw RipleyException("Illegal domain of normal locations");          case Nodes: return "Ripley_Nodes";
75      Ripley_Mesh* mesh=m_ripleyMesh.get();          case ReducedNodes: return "Ripley_Reduced_Nodes";
76      escriptDataC _normal=normal.getDataC();          case Elements: return "Ripley_Elements";
77      switch(normal.getFunctionSpace().getTypeCode()) {          case ReducedElements: return "Ripley_Reduced_Elements";
78          case Nodes:          case FaceElements: return "Ripley_Face_Elements";
79              throw RipleyException("Ripley does not support surface normal vectors for nodes");          case ReducedFaceElements: return "Ripley_Reduced_Face_Elements";
80          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.");  
81          default:          default:
82              stringstream temp;              break;
             temp << "Normal Vectors: Ripley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();  
             throw RipleyException(temp.str());  
83      }      }
84  */      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");  
85  }  }
86    
87  //  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  
88  {  {
89  /*      const int ptsPerSample = (m_numDim==2 ? 4 : 8);
90      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()) {  
91          case Nodes:          case Nodes:
92          case ReducedNodes:              return pair<int,int>(1, getNumNodes());
93          case DegreesOfFreedom:          case DegreesOfFreedom:
94          case ReducedDegreesOfFreedom:              return pair<int,int>(1, getNumNodes());
             temp=escript::Data( arg, escript::function(*this) );  
             _temp=temp.getDataC();  
             Ripley_Assemble_integrate(m_nodes,mesh->Elements,&_temp,&integrals[0]);  
             break;  
95          case Elements:          case Elements:
96          case ReducedElements:              return pair<int,int>(ptsPerSample, getNumElements());
             Ripley_Assemble_integrate(m_nodes,mesh->Elements,&_arg,&integrals[0]);  
             break;  
97          case FaceElements:          case FaceElements:
98          case ReducedFaceElements:              return pair<int,int>(ptsPerSample/2, getNumFaceElements());
99              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.");  
100          case ReducedNodes:          case ReducedNodes:
101              throw RipleyException("Gradient at reduced nodes is not supported.");              return pair<int,int>(1, getNumReducedNodes());
         case Elements:  
             Ripley_Assemble_gradient(m_nodes,mesh->Elements,&_grad,&nodeDataC);  
             break;  
102          case ReducedElements:          case ReducedElements:
103              Ripley_Assemble_gradient(m_nodes,mesh->Elements,&_grad,&nodeDataC);              return pair<int,int>(1, getNumReducedElements());
             break;  
         case FaceElements:  
             Ripley_Assemble_gradient(m_nodes,mesh->FaceElements,&_grad,&nodeDataC);  
             break;  
104          case ReducedFaceElements:          case ReducedFaceElements:
105              Ripley_Assemble_gradient(m_nodes,mesh->FaceElements,&_grad,&nodeDataC);              return pair<int,int>(1, getNumReducedFaceElements());
             break;  
106          case Points:          case Points:
107              throw RipleyException("Gradient at points is not supported");              return pair<int,int>(1, getNumNodes());
         case DegreesOfFreedom:  
             throw RipleyException("Gradient at degrees of freedom is not supported");  
108          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
109              throw RipleyException("Gradient at reduced degrees of freedom is not supported");              return pair<int,int>(1, getNumReducedNodes());
110                */
111          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);  
112              break;              break;
         case ReducedElements:  
             Ripley_Assemble_getSize(m_nodes, mesh->Elements, &tmp);  
             break;  
         case FaceElements:  
             Ripley_Assemble_getSize(m_nodes, mesh->FaceElements, &tmp);  
             break;  
         case ReducedFaceElements:  
             Ripley_Assemble_getSize(m_nodes, mesh->FaceElements, &tmp);  
             break;  
         case Points:  
             throw RipleyException("Size of point elements is not supported");  
         case DegreesOfFreedom:  
             throw RipleyException("Size of degrees of freedom is not supported");  
         case ReducedDegreesOfFreedom:  
             throw RipleyException("Size of reduced degrees of freedom is not supported");  
         default:  
             stringstream temp;  
             temp << "Element size: Ripley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();  
             throw RipleyException(temp.str());  
113      }      }
 */  
 }  
   
 //  
 // 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");  
114    
115      if (new_x.getFunctionSpace() == continuousFunction(*this)) {      stringstream msg;
116          m_nodes->setCoordinates(new_x);      msg << "Invalid function space type " << fsType << " for "
117      } else {          << getDescription();
118          escript::Data new_x_inter=escript::Data(new_x, continuousFunction(*this) );      throw RipleyException(msg.str());
         m_nodes->setCoordinates(new_x_inter);  
     }  
119  }  }
120    
121  bool RipleyDomain::ownSample(int fs_code, index_t id) const  int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
122  {  {
123  #ifdef ESYS_MPI      throw RipleyException("getTagFromSampleNo() not implemented");
     index_t myFirstNode, myLastNode, k;  
     if (fs_code == ReducedNodes) {  
         myFirstNode = m_nodes->getFirstReducedNode();  
         myLastNode = m_nodes->getLastReducedNode();  
         k = m_nodes->getIndexForGlobalReducedNode(id);  
     } else {  
         myFirstNode = m_nodes->getFirstNode();  
         myLastNode = m_nodes->getLastNode();  
         k = m_nodes->getIndexForGlobalNode(id);  
     }  
     return static_cast<bool>((myFirstNode <= k) && (k < myLastNode));  
 #endif  
     return true;  
124  }  }
125    
126    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  
127  {  {
128      switch(functionSpaceCode) {      stringstream ret;
129          case Nodes:      TagMap::const_iterator it;
130          case DegreesOfFreedom:      for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
131          case ReducedDegreesOfFreedom:          if (it!=m_tagMap.begin()) ret << ", ";
132              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());  
133      }      }
134      return false;      return ret.str();
135  }  }
136    
137  bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const  bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
138  {  {
139     /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]     /*
140        The idea is to use equivalence classes (i.e. types which can be
141        interpolated back and forth):
142      class 1: DOF <-> Nodes      class 1: DOF <-> Nodes
143      class 2: ReducedDOF <-> ReducedNodes      class 2: ReducedDOF <-> ReducedNodes
144      class 3: Points      class 3: Points
# Line 1227  bool RipleyDomain::commonFunctionSpace(c Line 149  bool RipleyDomain::commonFunctionSpace(c
149      class 8: ContactElementZero <-> ContactElementOne      class 8: ContactElementZero <-> ContactElementOne
150      class 9: ReducedContactElementZero <-> ReducedContactElementOne      class 9: ReducedContactElementZero <-> ReducedContactElementOne
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        class 1 and 2 belong to all lines so aren't considered.
155      line 0: class 3      line 0: class 3
156      line 1: class 4,5      line 1: class 4,5
157      line 2: class 6,7      line 2: class 6,7
# Line 1286  bool RipleyDomain::commonFunctionSpace(c Line 209  bool RipleyDomain::commonFunctionSpace(c
209          return false;          return false;
210      else if (totlines==1) {      else if (totlines==1) {
211          // we have points          // we have points
212          if (hasline[0]==1) {          if (hasline[0]==1)
213              resultcode=Points;              resultcode=Points;
214          } else if (hasline[1]==1) {          else if (hasline[1]==1) {
215              if (hasclass[5]==1) {              if (hasclass[5]==1)
216                  resultcode=ReducedElements;                  resultcode=ReducedElements;
217              } else {              else
218                  resultcode=Elements;                  resultcode=Elements;
             }  
219          } else if (hasline[2]==1) {          } else if (hasline[2]==1) {
220              if (hasclass[7]==1) {              if (hasclass[7]==1)
221                  resultcode=ReducedFaceElements;                  resultcode=ReducedFaceElements;
222              } else {              else
223                  resultcode=FaceElements;                  resultcode=FaceElements;
224              }          } else
225          } else { // so we must be in line3              throw RipleyException("Internal Ripley Error!");
             throw RipleyException("Programmer Error - choosing between contact elements - we should never get here");  
         }  
226      } else { // totlines==0      } else { // totlines==0
227          if (hasclass[2]==1) {          if (hasclass[2]==1)
228              // something from class 2              // something from class 2
229              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
230          } else { // something from class 1          else // something from class 1
231              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
         }  
232      }      }
233      return true;      return true;
234  }  }
235    
236  bool RipleyDomain::probeInterpolationOnDomain(int functionSpaceType_source,  escript::Data RipleyDomain::getX() const
237                                               int functionSpaceType_target) const  {
238        return escript::continuousFunction(*this).getX();
239    }
240    
241    escript::Data RipleyDomain::getNormal() const
242    {
243        return escript::functionOnBoundary(*this).getNormal();
244    }
245    
246    escript::Data RipleyDomain::getSize() const
247    {
248        return escript::function(*this).getSize();
249    }
250    
251    void RipleyDomain::setToX(escript::Data& arg) const
252  {  {
253      if (functionSpaceType_target != Nodes &&      const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
254              functionSpaceType_target != ReducedNodes &&              *(arg.getFunctionSpace().getDomain()));
255              functionSpaceType_target != ReducedDegreesOfFreedom &&      if (argDomain != *this)
256              functionSpaceType_target != DegreesOfFreedom &&          throw RipleyException("setToX: Illegal domain of data point locations");
257              functionSpaceType_target != Elements &&      if (!arg.isExpanded())
258              functionSpaceType_target != ReducedElements &&          throw RipleyException("setToX: Expanded Data object expected");
259              functionSpaceType_target != FaceElements &&  
260              functionSpaceType_target != ReducedFaceElements &&      if (arg.getFunctionSpace().getTypeCode()==Nodes) {
261              functionSpaceType_target != Points) {          assembleCoordinates(arg);
262          stringstream temp;      } else {
263          temp << "Interpolation On Domain: Ripley does not know anything about function space type " << functionSpaceType_target;          // interpolate the result
264          throw RipleyException(temp.str());          escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
265            assembleCoordinates(contData);
266            interpolateOnDomain(arg, contData);
267      }      }
268    }
269    
270      switch (functionSpaceType_source) {  bool RipleyDomain::isCellOriented(int fsType) const
271    {
272        switch(fsType) {
273          case Nodes:          case Nodes:
274          case DegreesOfFreedom:          case DegreesOfFreedom:
             return true;  
         case ReducedNodes:  
275          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
276              return (functionSpaceType_target != Nodes &&              return false;
                     functionSpaceType_target != DegreesOfFreedom);  
277          case Elements:          case Elements:
             return (functionSpaceType_target==Elements ||  
                     functionSpaceType_target==ReducedElements);  
         case ReducedElements:  
             return (functionSpaceType_target==ReducedElements);  
278          case FaceElements:          case FaceElements:
             return (functionSpaceType_target==FaceElements ||  
                     functionSpaceType_target==ReducedFaceElements);  
         case ReducedFaceElements:  
             return (functionSpaceType_target==ReducedFaceElements);  
279          case Points:          case Points:
280              return (functionSpaceType_target==Points);          case ReducedElements:
281            case ReducedFaceElements:
282                return true;
283          default:          default:
284              stringstream temp;              break;
             temp << "Interpolation On Domain: Ripley does not know anything about function space type " << functionSpaceType_source;  
             throw RipleyException(temp.str());  
285      }      }
286      return false;      stringstream msg;
287        msg << "Illegal function space type " << fsType << " on " << getDescription();
288        throw RipleyException(msg.str());
289  }  }
290    
291  bool RipleyDomain::probeInterpolationACross(int functionSpaceType_source,  void RipleyDomain::Print_Mesh_Info(const bool full) const
         const AbstractDomain& targetDomain, int functionSpaceType_target) const  
292  {  {
293      return false;      cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
294  }          << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
295        cout << "Number of dimensions: " << m_numDim << endl;
296    
297  bool RipleyDomain::operator==(const AbstractDomain& other) const      // write tags
298  {      if (m_tagMap.size() > 0) {
299      const RipleyDomain* temp=dynamic_cast<const RipleyDomain*>(&other);          cout << "Tags:" << endl;
300      if (temp != NULL) {          TagMap::const_iterator it;
301          return (m_name==temp->m_name && m_nodes==temp->m_nodes &&          for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
302                  m_elements==temp->m_elements &&              cout << "  " << setw(5) << it->second << " "
303                  m_faceElements==temp->m_faceElements &&                  << it->first << endl;
304                  m_points==temp->m_points);          }
305      }      }
     return false;  
 }  
   
 bool RipleyDomain::operator!=(const AbstractDomain& other) const  
 {  
     return !(operator==(other));  
306  }  }
307    
308  int RipleyDomain::getSystemMatrixTypeId(const int solver,  int RipleyDomain::getSystemMatrixTypeId(const int solver,
309          const int preconditioner, const int package, const bool symmetry) const          const int preconditioner, const int package, const bool symmetry) const
310  {  {
311      int out=Paso_SystemMatrix_getSystemMatrixTypeId(      return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
312              SystemMatrixAdapter::mapOptionToPaso(solver),              package, symmetry, m_mpiInfo);
             SystemMatrixAdapter::mapOptionToPaso(preconditioner),  
             SystemMatrixAdapter::mapOptionToPaso(package),  
             symmetry?1:0, m_mpiInfo);  
     checkPasoError();  
     return out;  
313  }  }
314    
315  int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,  int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
316          const int package, const bool symmetry) const          const int package, const bool symmetry) const
317  {  {
318      int out=Paso_TransportProblem_getTypeId(      return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
319              SystemMatrixAdapter::mapOptionToPaso(solver),              package, symmetry, m_mpiInfo);
             SystemMatrixAdapter::mapOptionToPaso(preconditioner),  
             SystemMatrixAdapter::mapOptionToPaso(package),  
             symmetry?1:0, m_mpiInfo);  
     checkPasoError();  
     return out;  
320  }  }
321    
322  escript::Data RipleyDomain::getX() const  escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
323            const escript::FunctionSpace& row_functionspace,
324            const int column_blocksize,
325            const escript::FunctionSpace& column_functionspace,
326            const int type) const
327  {  {
328      return continuousFunction(*this).getX();      bool reduceRowOrder=false;
329        bool reduceColOrder=false;
330        // is the domain right?
331        const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
332        if (row_domain!=*this)
333            throw RipleyException("Domain of row function space does not match the domain of matrix generator");
334        const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
335        if (col_domain!=*this)
336            throw RipleyException("Domain of column function space does not match the domain of matrix generator");
337        // is the function space type right?
338        if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
339            reduceRowOrder=true;
340        else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
341            throw RipleyException("Illegal function space type for system matrix rows");
342        if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
343            reduceColOrder=true;
344        else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
345            throw RipleyException("Illegal function space type for system matrix columns");
346    
347        // generate matrix
348        Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
349        Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
350                row_blocksize, column_blocksize, FALSE);
351        paso::checkPasoError();
352        Paso_SystemMatrixPattern_free(pattern);
353        escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
354                    row_functionspace, column_blocksize, column_functionspace));
355        return sma;
356  }  }
357    
358  escript::Data RipleyDomain::getNormal() const  //
359    // the methods that follow have to be implemented by the subclasses
360    //
361    
362    string RipleyDomain::getDescription() const
363  {  {
364      return functionOnBoundary(*this).getNormal();      throw RipleyException("getDescription() not implemented");
365  }  }
366    
367  escript::Data RipleyDomain::getSize() const  bool RipleyDomain::operator==(const AbstractDomain& other) const
368  {  {
369      return escript::function(*this).getSize();      throw RipleyException("operator==() not implemented");
370  }  }
371    
372  const int* RipleyDomain::borrowSampleReferenceIDs(int functionSpaceType) const  void RipleyDomain::write(const string& filename) const
373  {  {
374      switch (functionSpaceType) {      throw RipleyException("write() not implemented");
         case Nodes:  
             return &m_nodes->getIdVector()[0];  
         case ReducedNodes:  
             return &m_nodes->getReducedNodesId()[0];  
         case Elements:  
         case ReducedElements:  
             return &m_elements->getIdVector()[0];  
         case FaceElements:  
         case ReducedFaceElements:  
             return &m_faceElements->getIdVector()[0];  
         case Points:  
             return &m_points->getIdVector()[0];  
         case DegreesOfFreedom:  
             return &m_nodes->getDegreesOfFreedomId()[0];  
         case ReducedDegreesOfFreedom:  
             return &m_nodes->getReducedDegreesOfFreedomId()[0];  
         default:  
             stringstream msg;  
             msg << "borrowSampleReferenceIDs: Invalid function space type " << functionSpaceType << " for domain: " << getDescription();  
             throw RipleyException(msg.str());  
     }  
375  }  }
376    
377  int RipleyDomain::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  void RipleyDomain::dump(const string& filename) const
378  {  {
379      throw RipleyException("not implemented");      throw RipleyException("dump() not implemented");
380  }  }
381    
382    const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
 void RipleyDomain::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const  
383  {  {
384      switch (functionSpaceType) {      throw RipleyException("borrowSampleReferenceIDs() not implemented");
         case Nodes:  
             m_nodes->setTags(newTag, mask);  
             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");  
         case Elements:  
             m_elements->setTags(newTag, mask);  
             break;  
         case ReducedElements:  
             m_elements->setTags(newTag, mask);  
             break;  
         case FaceElements:  
             m_faceElements->setTags(newTag, mask);  
             break;  
         case ReducedFaceElements:  
             m_faceElements->setTags(newTag, mask);  
             break;  
         case Points:  
             m_points->setTags(newTag, mask);  
             break;  
         default:  
             stringstream msg;  
             msg << "Ripley does not know anything about function space type " << functionSpaceType;  
             throw RipleyException(msg.str());  
     }  
385  }  }
386    
387  void RipleyDomain::setTagMap(const string& name, int tag)  void RipleyDomain::setNewX(const escript::Data& arg)
388  {  {
389      m_tagMap[name] = tag;      throw RipleyException("setNewX(): Operation not supported");
390  }  }
391    
392  int RipleyDomain::getTag(const string& name) const  void RipleyDomain::interpolateOnDomain(escript::Data& target,
393                                           const escript::Data& in) const
394  {  {
395      return m_tagMap.find(name)->second;      throw RipleyException("interpolateOnDomain() not implemented");
396  }  }
397    
398  bool RipleyDomain::isValidTagName(const string& name) const  bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
399                                                 int fsType_target) const
400  {  {
401      return (m_tagMap.find(name)!=m_tagMap.end());      throw RipleyException("probeInterpolationOnDomain() not implemented");
402  }  }
403    
404  string RipleyDomain::showTagNames() const  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
405  {  {
406      stringstream ret;      throw RipleyException("interpolateACross() not implemented");
     TagMap::const_iterator it;  
     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {  
         if (it!=m_tagMap.begin()) ret << ", ";  
         ret << it->first;  
     }  
     return ret.str();  
407  }  }
408    
409  int RipleyDomain::getNumberOfTagsInUse(int functionSpaceCode) const  bool RipleyDomain::probeInterpolationACross(int fsType_source,
410            const escript::AbstractDomain&, int fsType_target) const
411  {  {
412      dim_t numTags=0;      throw RipleyException("probeInterpolationACross() not implemented");
     switch (functionSpaceCode) {  
         case Nodes:  
             numTags=m_nodes->getNumberOfTagsInUse();  
             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");  
         case Elements:  
         case ReducedElements:  
             numTags=m_elements->getNumberOfTagsInUse();  
             break;  
         case FaceElements:  
         case ReducedFaceElements:  
             numTags=m_faceElements->getNumberOfTagsInUse();  
             break;  
         case Points:  
             numTags=m_points->getNumberOfTagsInUse();  
             break;  
         default:  
             stringstream msg;  
             msg << "Ripley does not know anything about function space type " << functionSpaceCode;  
             throw RipleyException(msg.str());  
     }  
     return numTags;  
413  }  }
414    
415  const int* RipleyDomain::borrowListOfTagsInUse(int functionSpaceCode) const  void RipleyDomain::setToNormal(escript::Data& normal) const
416  {  {
417      switch (functionSpaceCode) {      throw RipleyException("setToNormal() not implemented");
         case Nodes:  
             return &m_nodes->getTagsInUse()[0];  
         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");  
         case Elements:  
         case ReducedElements:  
             return &m_elements->getTagsInUse()[0];  
         case FaceElements:  
         case ReducedFaceElements:  
             return &m_faceElements->getTagsInUse()[0];  
         case Points:  
             return &m_points->getTagsInUse()[0];  
         default:  
             stringstream msg;  
             msg << "Ripley does not know anything about function space type " << functionSpaceCode;  
             throw RipleyException(msg.str());  
     }  
418  }  }
419    
420  bool RipleyDomain::canTag(int functionSpaceCode) const  void RipleyDomain::setToSize(escript::Data& size) const
421  {  {
422      switch(functionSpaceCode) {      throw RipleyException("setToSize() not implemented");
         case Nodes:  
         case Elements:  
         case ReducedElements:  
         case FaceElements:  
         case ReducedFaceElements:  
         case Points:  
             return true;  
     }  
     return false;  
423  }  }
424    
425  escript::AbstractDomain::StatusType RipleyDomain::getStatus() const  void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
426  {  {
427      return m_nodes->getStatus();      throw RipleyException("setToGradient() not implemented");
428  }  }
429    
430  void RipleyDomain::prepare(bool optimize)  void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
431  {  {
432      resolveNodeIds();      throw RipleyException("setToIntegrals() not implemented");
   
     // first step is to distribute the elements according to a global  
     // distribution of DOF  
     // - create dense labeling for the DOFs  
     dim_t newGlobalNumDOFs=m_nodes->createDenseDOFLabeling();  
   
     // - 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);  
         }  
   
         // optimize the local labeling of the degrees of freedom  
         optimizeDOFLabeling(distribution);  
     }  
   
     // Rearrange elements in order to try and bring them closer to memory  
     // locations of the nodes (distributed shared memory).  
     optimizeElementOrdering();  
   
 /* useful DEBUG:  
     {  
         printf("prepare: global DOF : %d\n",newGlobalNumDOFs);  
         IndexPair range = m_nodes->getGlobalIdRange();  
         printf("prepare: global node id range = %d :%d\n",range.first,range.second);  
         range = m_nodes->getIdRange();  
         printf("prepare: local node id range = %d :%d\n",range.first,range.second);  
     }  
 */  
   
     // create the global indices  
     IndexVector maskReducedNodes(m_nodes->getNumNodes(), -1);  
     markNodes(maskReducedNodes, 0);  
     IndexVector indexReducedNodes = packMask(maskReducedNodes);  
   
     IndexVector nodeDistribution(m_mpiInfo->size+1);  
     m_nodes->createDenseNodeLabeling(nodeDistribution, distribution);  
     m_nodes->createDenseReducedDOFLabeling(maskReducedNodes);  
     m_nodes->createDenseReducedNodeLabeling(maskReducedNodes);  
   
     // create the missing mappings  
     m_nodes->createNodeFileMappings(indexReducedNodes, distribution,  
                                     nodeDistribution);  
   
     updateTagsInUse();  
433  }  }
434    
435  void RipleyDomain::optimizeElementOrdering()  bool RipleyDomain::ownSample(int fsType, index_t id) const
436  {  {
437      m_elements->optimizeOrdering();      throw RipleyException("ownSample() not implemented");
     m_faceElements->optimizeOrdering();  
     m_points->optimizeOrdering();  
438  }  }
439    
440  void RipleyDomain::updateTagsInUse()  void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& mask) const
441  {  {
442      m_nodes->updateTagsInUse();      throw RipleyException("setTags() not implemented");
     m_elements->updateTagsInUse();  
     m_faceElements->updateTagsInUse();  
     m_points->updateTagsInUse();  
443  }  }
444    
445  void RipleyDomain::createColoring(const IndexVector &node_localDOF_map)  int RipleyDomain::getNumberOfTagsInUse(int fsType) const
446  {  {
447      m_elements->createColoring(node_localDOF_map);      throw RipleyException("getNumberOfTagsInUse() not implemented");
     m_faceElements->createColoring(node_localDOF_map);  
     m_points->createColoring(node_localDOF_map);  
448  }  }
449    
450  void RipleyDomain::distributeByRankOfDOF(const IndexVector &dofDistribution)  const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
451  {  {
452      RankVector mpiRankOfDOF = m_nodes->getOwnerOfDOFs(dofDistribution);      throw RipleyException("borrowListOfTagsInUse() not implemented");
   
     // First the elements are redistributed according to mpiRankOfDOF.  
     // At the input the Node tables refer to the local labeling of  
     // the nodes while at the output they refer to the global labeling  
     // which is rectified in the next step  
     m_elements->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());  
     m_faceElements->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());  
     m_points->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());  
   
     // resolve the node IDs  
     resolveNodeIds();  
   
     // create a local labeling of the DOFs  
     IndexPair dofRange = m_nodes->getDOFRange();  
     dim_t len = dofRange.second - dofRange.first + 1;  
   
     // local mask for used nodes  
     IndexVector tmp_node_localDOF_mask(len, -1);  
     IndexVector tmp_node_localDOF_map(m_nodes->getNumNodes(), -1);  
   
 #pragma omp parallel for schedule(static)  
     for (dim_t n = 0; n < m_nodes->getNumNodes(); n++) {  
 #ifdef BOUNDS_CHECK  
         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;  
     }  
   
     dim_t numDOFs = 0;  
     for (dim_t n = 0; n < len; n++) {  
         if (tmp_node_localDOF_mask[n] >= 0) {  
             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);  
453  }  }
454    
455  /**************************************************************  bool RipleyDomain::canTag(int fsType) const
456     Check whether there is any node which has no vertex. In case  {
457     such a node exists, we don't use ParMetis since it requires      throw RipleyException("canTag() not implemented");
458     that every node has at least 1 vertex (at line 129 of file  }
    "xyzpart.c" in parmetis 3.1.1, variable "nvtxs" would be 0 if  
    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;  
 }  
 #endif  
   
 void RipleyDomain::optimizeDOFDistribution(RankVector &distribution)  
 {  
     const Esys_MPI_rank myRank = m_mpiInfo->rank;  
     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)]);  
             }  
         }  
   
         IndexMatrix indexList(myNumVertices);  
         // ksteube CSR of DOF IDs  
         // 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, 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  
   
         Paso_Pattern_free(pattern);  
   
         // create a new distribution and labeling of the DOF  
         RankVector newDistribution(mpiSize+1, 0);  
 #pragma omp parallel  
         {  
             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];  
             }  
         }  
459    
         // recvbuf will be the concatenation of each CPU's contribution  
         // to newDistribution  
         dim_t *recvbuf = TMPMEMALLOC(mpiSize * mpiSize, dim_t);  
   
 #ifdef ESYS_MPI  
         MPI_Allgather(&newDistribution[0], mpiSize, MPI_INT, recvbuf, mpiSize, MPI_INT, m_mpiInfo->comm);  
 #else  
         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);  
460    
461          // now the overlap needs to be created by sending the partition  void RipleyDomain::addPDEToSystem(
462          // around          escript::AbstractSystemMatrix& mat, escript::Data& rhs,
463          m_nodes->resetGlobalDegreesOfFreedom(newGlobalDOFid, distribution);          const escript::Data& A, const escript::Data& B, const escript::Data& C,
464          for (dim_t i = 0; i < mpiSize+1; ++i)          const escript::Data& D, const escript::Data& X, const escript::Data& Y,
465              distribution[i] = newDistribution[i];          const escript::Data& d, const escript::Data& y,
466      }          const escript::Data& d_contact, const escript::Data& y_contact,
467      TMPMEMFREE(partition);          const escript::Data& d_dirac,const escript::Data& y_dirac) const
468      TMPMEMFREE(xyz);  {
469        throw RipleyException("addPDEToSystem() not implemented");
470  }  }
471    
472  void RipleyDomain::optimizeDOFLabeling(const IndexVector &distribution)  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
473            const escript::Data& D, const escript::Data& d,
474            const escript::Data& d_dirac, const bool useHRZ) const
475  {  {
476      Esys_MPI_rank myRank = m_mpiInfo->rank;      throw RipleyException("addPDEToLumpedSystem() not implemented");
477      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);  
   
     Esys_MPIInfo_noError(m_mpiInfo);  
     checkPasoError();  
   
     /* shift new labeling to create a global id */  
 #pragma omp parallel for  
     for (dim_t i = 0; i < myNumVertices; ++i)  
         newGlobalDOFid[i] += myFirstVertex;  
   
     /* distribute new labeling to other processors */  
     m_nodes->resetGlobalDegreesOfFreedom(newGlobalDOFid, distribution);  
 #if 0  
     for (i = 0; i < m_nodes->getNumNodes(); ++i)  
         printf("%d ", m_nodes->getGlobalDegreesOfFreedom()[i]);  
     printf("\n");  
 #endif  
 }  
   
 void RipleyDomain::resolveNodeIds()  
 {  
     // find the minimum and maximum node ID used by elements  
     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);  
   
 #ifdef Ripley_TRACE  
     index_t globalMinId, globalMaxId;  
 #ifdef ESYS_MPI  
     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;  
     }  
   
     // create a new node file  
     m_nodes = m_nodes->gatherGlobal(newLocalToGlobalNodeLabels);  
478    
479      // relabel nodes of the elements  void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
480      relabelElementNodes(globalToNewLocalNodeLabels, minId);          const escript::Data& Y, const escript::Data& y,
481            const escript::Data& y_contact, const escript::Data& y_dirac) const
482    {
483        throw RipleyException("addPDEToRHS() not implemented");
484  }  }
485    
486  void RipleyDomain::relabelElementNodes(const IndexVector &newNode, index_t offset)  void RipleyDomain::addPDEToTransportProblem(
487            escript::AbstractTransportProblem& tp,
488            escript::Data& source, const escript::Data& M,
489            const escript::Data& A, const escript::Data& B, const escript::Data& C,
490            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
491            const escript::Data& d, const escript::Data& y,
492            const escript::Data& d_contact, const escript::Data& y_contact,
493            const escript::Data& d_dirac, const escript::Data& y_dirac) const
494  {  {
495      m_elements->relabelNodes(newNode, offset);      throw RipleyException("addPDEToTransportProblem() not implemented");
     m_faceElements->relabelNodes(newNode, offset);  
     m_points->relabelNodes(newNode, offset);  
496  }  }
497    
498  void RipleyDomain::markNodes(IndexVector &mask, index_t offset)  escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
499            const int blocksize, const escript::FunctionSpace& functionspace,
500            const int type) const
501  {  {
502      m_elements->markNodes(mask, offset);      throw RipleyException("newTransportProblem() not implemented");
     m_faceElements->markNodes(mask, offset);  
     m_points->markNodes(mask, offset);  
503  }  }
504    
505  void RipleyDomain::createMappings(const IndexVector &dofDistribution, const IndexVector &nodeDistribution)  Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
506                                                       bool reducedColOrder) const
507  {  {
508      IndexVector maskReducedNodes(m_nodes->getNumNodes(), -1);      throw RipleyException("getPattern() not implemented");
509      markNodes(maskReducedNodes, 0);  }
510    
511      IndexVector indexReducedNodes = packMask(maskReducedNodes);  dim_t RipleyDomain::getNumDataPointsGlobal() const
512      m_nodes->createNodeFileMappings(indexReducedNodes, dofDistribution,  {
513                                      nodeDistribution);      throw RipleyException("getNumDataPointsGlobal() not implemented");
514  }  }
515    
516  Paso_SystemMatrixPattern *RipleyDomain::makePattern(bool reduce_row_order, bool reduce_col_order) const  dim_t RipleyDomain::getNumFaceElements() const
517  {  {
518      Paso_Connector *colConnector, *rowConnector;      throw RipleyException("getNumFaceElements() not implemented");
519      NodeMapping_ptr colMap, rowMap;  }
     Paso_Distribution *colDistribution = NULL, *rowDistribution = NULL;  
520    
521      if (reduce_col_order) {  dim_t RipleyDomain::getNumElements() const
522          colMap = m_nodes->getReducedDegreesOfFreedomMapping();  {
523          colDistribution = m_nodes->getReducedDegreesOfFreedomDistribution();      throw RipleyException("getNumElements() not implemented");
524          colConnector = m_nodes->getReducedDegreesOfFreedomConnector();  }
     } else {  
         colMap = m_nodes->getDegreesOfFreedomMapping();  
         colDistribution = m_nodes->getDegreesOfFreedomDistribution();  
         colConnector = m_nodes->getDegreesOfFreedomConnector();  
     }  
525    
526      if (reduce_row_order) {  dim_t RipleyDomain::getNumNodes() const
527          rowMap = m_nodes->getReducedDegreesOfFreedomMapping();  {
528          rowDistribution = m_nodes->getReducedDegreesOfFreedomDistribution();      throw RipleyException("getNumNodes() not implemented");
529          rowConnector = m_nodes->getReducedDegreesOfFreedomConnector();  }
     } else {  
         rowMap = m_nodes->getDegreesOfFreedomMapping();  
         rowDistribution = m_nodes->getDegreesOfFreedomDistribution();  
         rowConnector = m_nodes->getDegreesOfFreedomConnector();  
     }  
530    
531      //  insert contributions from element matrices into columns in indexList  void RipleyDomain::assembleCoordinates(escript::Data& arg) const
532      IndexMatrix indexList(rowMap->numTargets);  {
533      m_elements->insertIntoIndexMatrix(indexList, rowMap->target, colMap->target);      throw RipleyException("assembleCoordinates() not implemented");
     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;  
     */  
534  }  }
535    
536    
537  } // end of namespace ripley  } // end of namespace ripley
538    

Legend:
Removed from v.3690  
changed lines
  Added in v.3691

  ViewVC Help
Powered by ViewVC 1.1.26