Transfer operators moved to new module dune-solvers

[[Imported from SVN: r2341]]

dune-solvers/transferoperators/compressedmultigridtransfer.hh

+#ifndef COMPRESSED_MULTIGRID_TRANSFER_HH
+#define COMPRESSED_MULTIGRID_TRANSFER_HH
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include "genericmultigridtransfer.hh"
+#include "multigridtransfer.hh"
+
+
+/** \brief Restriction and prolongation operator for standard multigrid
+ *
+ * This class implements the standard prolongation and restriction
+ * operators for standard multigrid solvers.  Restriction and prolongation
+ * of block vectors is provided.  Internally, the operator is stored
+ * as a scalar BCRSMatrix.  Therefore, the template parameter VectorType
+ * has to comply with the ISTL requirements.
+
+ * \todo Currently only works for first-order Lagrangian elements!
+ */
+template<
+    class VectorType,
+    class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension>,
+    class MatrixType = Dune::BCRSMatrix< typename Dune::FieldMatrix<
+        typename VectorType::field_type, VectorType::block_type::size, VectorType::block_type::size> > >
+class CompressedMultigridTransfer :
+    public MultigridTransfer<VectorType,BitVectorType,MatrixType>
+{
+
+    enum {blocksize = VectorType::block_type::size};
+
+    typedef typename VectorType::field_type field_type;
+
+public:
+
+    typedef MatrixType OperatorType;
+
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type,1,1> > TransferOperatorType;
+
+    CompressedMultigridTransfer()
+    {
+        matrix_ = new TransferOperatorType;
+        matrixInternallyAllocated = true;
+    }
+
+    virtual ~CompressedMultigridTransfer()
+    {
+        if (matrixInternallyAllocated)
+            delete matrix_;
+    }
+
+
+    /** \brief Sets up the operator between two P1 spaces
+     *
+     * \param grid The grid
+     * \param cL The coarse grid level
+     * \param fL The fine grid level
+     */
+    template <class GridType>
+    void setup(const GridType& grid, int cL, int fL)
+    {
+        GenericMultigridTransfer::setup<TransferOperatorType, GridType, field_type>(*matrix_, grid, cL, fL);
+    }
+
+
+    /** \brief Restrict a function from the fine onto the coarse grid
+     */
+    virtual void restrict(const VectorType& f, VectorType &t) const
+    {
+        GenericMultigridTransfer::restrict(*matrix_, f, t, -1);
+    }
+
+
+    /** \brief Restrict a bitfield from the fine onto the coarse grid
+     */
+    virtual void restrictScalarBitField(const Dune::BitSetVector<1>& f, Dune::BitSetVector<1>& t) const
+    {
+        GenericMultigridTransfer::restrictScalarBitField(*matrix_, f, t, -1);
+    }
+
+
+    /** \brief Restrict a vector valued bitfield from the fine onto the coarse grid
+     */
+    virtual void restrict(const BitVectorType& f, BitVectorType& t) const
+    {
+        GenericMultigridTransfer::restrictBitField(*matrix_, f, t, -1);
+    }
+
+
+    /** \brief Restrict a vector valued bitfield from the fine onto the coarse grid
+     * Fine bits only influence their father bits
+     */
+    virtual void restrictToFathers(const BitVectorType& f, BitVectorType& t) const
+    {
+        GenericMultigridTransfer::restrictBitFieldToFathers(*matrix_, f, t, -1);
+    }
+
+
+    /** \brief Prolong a function from the coarse onto the fine grid
+     */
+    virtual void prolong(const VectorType& f, VectorType &t) const
+    {
+        GenericMultigridTransfer::prolong(*matrix_, f, t, -1);
+    }
+
+
+    /** \brief Galerkin assemble a coarse stiffness matrix
+     */
+    virtual void galerkinRestrict(const OperatorType& fineMat, OperatorType& coarseMat) const
+    {
+        GenericMultigridTransfer::galerkinRestrict(*matrix_, fineMat, coarseMat, -1);
+    }
+
+
+    /** \brief Set Occupation of Galerkin restricted coarse stiffness matrix   
+    *
+    * Set occupation of Galerkin restricted coarse matrix. Call this one before
+    * galerkinRestrict to ensure all non-zeroes are present
+    * \param fineMat The fine level matrix
+    * \param coarseMat The coarse level matrix
+    */
+    virtual void galerkinRestrictSetOccupation(const OperatorType& fineMat, OperatorType& coarseMat) const
+    {
+        GenericMultigridTransfer::galerkinRestrictSetOccupation(*matrix_, fineMat, coarseMat, -1);
+    }
+
+
+    /** \brief Direct access to the operator matrix, if you absolutely want it! */
+    const TransferOperatorType& getMatrix() const {
+        return *matrix_;
+    }
+
+    /** \brief Set matrix! */
+    void setMatrix(TransferOperatorType& matrix)
+    {
+        if (matrixInternallyAllocated)
+            delete matrix_;
+        matrixInternallyAllocated = false;
+        matrix_ = &matrix;
+    }
+
+
+
+protected:
+
+    TransferOperatorType* matrix_;
+    bool matrixInternallyAllocated;
+
+};
+
+#endif

dune-solvers/transferoperators/densemultigridtransfer.hh

+#ifndef DENSE_MULTIGRID_TRANSFER_HH
+#define DENSE_MULTIGRID_TRANSFER_HH
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include "multigridtransfer.hh"
+#include "genericmultigridtransfer.hh"
+
+
+/** \brief Restriction and prolongation operator for standard multigrid
+ *
+ * This class implements the standard prolongation and restriction
+ * operators for standard multigrid solvers.  Restriction and prolongation
+ * of block vectors is provided.  Internally, the operator is stored
+ * as a BCRSMatrix of small dense matrices.
+ * Therefore, the template parameter VectorType has to comply with the ISTL requirements.
+
+ * \todo Currently only works for first-order Lagrangian elements!
+ */
+template<
+    class VectorType,
+    class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension>,
+    class MatrixType = Dune::BCRSMatrix< typename Dune::FieldMatrix<
+        typename VectorType::field_type, VectorType::block_type::size, VectorType::block_type::size> > >
+class DenseMultigridTransfer :
+    public MultigridTransfer<VectorType,BitVectorType,MatrixType>
+{
+
+    enum {blocksize = VectorType::block_type::size};
+
+    typedef typename VectorType::field_type field_type;
+
+public:
+
+    typedef Dune::BCRSMatrix< Dune::FieldMatrix< field_type, blocksize, blocksize > > TransferOperatorType;
+
+    virtual ~DenseMultigridTransfer() {}
+
+
+    /** \brief Sets up the operator between two P1 spaces
+     *
+     * \param grid The grid
+     * \param cL The coarse grid level
+     * \param fL The fine grid level
+     */
+    template <class GridType>
+    void setup(const GridType& grid, int cL, int fL)
+    {
+        GenericMultigridTransfer::setup<TransferOperatorType, GridType, field_type>(matrix_, grid, cL, fL);
+    }
+
+
+    /** \brief Restrict a function from the fine onto the coarse grid
+     */
+    virtual void restrict(const VectorType& f, VectorType &t) const
+    {
+        GenericMultigridTransfer::restrict(matrix_, f, t);
+    }
+
+
+    /** \brief Restrict a bitfield from the fine onto the coarse grid
+     */
+    virtual void restrictScalarBitField(const Dune::BitSetVector<1>& f, Dune::BitSetVector<1>& t) const
+    {
+        GenericMultigridTransfer::restrictScalarBitField(matrix_, f, t);
+    }
+
+
+    /** \brief Restrict a vector valued bitfield from the fine onto the coarse grid
+     */
+    virtual void restrict(const BitVectorType& f, BitVectorType& t) const
+    {
+        GenericMultigridTransfer::restrictBitField(matrix_, f, t);
+    }
+
+
+    /** \brief Restrict a vector valued bitfield from the fine onto the coarse grid
+     * Fine bits only influence their father bits
+     */
+    virtual void restrictToFathers(const BitVectorType& f, BitVectorType& t) const
+    {
+        GenericMultigridTransfer::restrictBitFieldToFathers(matrix_, f, t);
+    }
+
+
+    /** \brief Prolong a function from the coarse onto the fine grid
+     */
+    virtual void prolong(const VectorType& f, VectorType &t) const
+    {
+        GenericMultigridTransfer::prolong(matrix_, f, t);
+    }
+
+
+    /** \brief Galerkin assemble a coarse stiffness matrix
+     */
+    virtual void galerkinRestrict(const MatrixType& fineMat, MatrixType& coarseMat) const
+    {
+        GenericMultigridTransfer::galerkinRestrict(matrix_, fineMat, coarseMat);
+    }
+
+
+    /** \brief Set Occupation of Galerkin restricted coarse stiffness matrix   
+    *
+    * Set occupation of Galerkin restricted coarse matrix. Call this one before
+    * galerkinRestrict to ensure all non-zeroes are present
+    * \param fineMat The fine level matrix
+    * \param coarseMat The coarse level matrix
+    */
+    virtual void galerkinRestrictSetOccupation(const MatrixType& fineMat, MatrixType& coarseMat) const
+    {
+        GenericMultigridTransfer::galerkinRestrictSetOccupation(matrix_, fineMat, coarseMat);
+    }
+
+
+    /** \brief Direct access to the operator matrix, if you absolutely want it! */
+    const TransferOperatorType& getMatrix() const {
+        return matrix_;
+    }
+
+
+protected:
+
+    TransferOperatorType matrix_;
+
+};
+
+#endif

dune-solvers/transferoperators/multigridtransfer.hh

+#ifndef MULTIGRID_TRANSFER_HH
+#define MULTIGRID_TRANSFER_HH
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/bitsetvector.hh>
+
+
+/** \brief Restriction and prolongation operator for standard multigrid
+ *
+ * This class is the base class for prolongation and restriction
+ * operators for standard multigrid solvers.  Restriction and prolongation
+ * of block vectors is provided.
+
+ * \todo Currently only works for first-order Lagrangian elements!
+ */
+template<
+    class VectorType,
+    class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension>,
+    class MatrixType = Dune::BCRSMatrix< typename Dune::FieldMatrix<
+        typename VectorType::field_type, VectorType::block_type::size, VectorType::block_type::size> > >
+class MultigridTransfer {
+
+public:
+
+    typedef MatrixType OperatorType;
+
+    virtual ~MultigridTransfer() {}
+
+
+    /** \brief Restrict a function from the fine onto the coarse grid
+     */
+    virtual void restrict(const VectorType& f, VectorType &t) const = 0;
+
+
+    /** \brief Restrict a bitfield from the fine onto the coarse grid
+     */
+    virtual void restrictScalarBitField(const Dune::BitSetVector<1>& f, Dune::BitSetVector<1>& t) const = 0;
+
+
+    /** \brief Restrict a vector valued bitfield from the fine onto the coarse grid
+     */
+    virtual void restrict(const BitVectorType& f, BitVectorType& t) const = 0;
+
+
+    /** \brief Restrict a vector valued bitfield from the fine onto the coarse grid
+     * Fine bits only influence their father bits
+     */
+    virtual void restrictToFathers(const BitVectorType& f, BitVectorType& t) const = 0;
+
+
+    /** \brief Prolong a function from the coarse onto the fine grid
+     */
+    virtual void prolong(const VectorType& f, VectorType &t) const = 0;
+
+
+    /** \brief Galerkin assemble a coarse stiffness matrix
+     */
+    virtual void galerkinRestrict(const MatrixType& fineMat, MatrixType& coarseMat) const = 0;
+
+
+    /** \brief Set Occupation of Galerkin restricted coarse stiffness matrix   
+    *
+    * Set occupation of Galerkin restricted coarse matrix. Call this one before
+    * galerkinRestrict to ensure all non-zeroes are present
+    * \param fineMat The fine level matrix
+    * \param coarseMat The coarse level matrix
+    */
+    virtual void galerkinRestrictSetOccupation(const MatrixType& fineMat, MatrixType& coarseMat) const = 0;
+
+};
+
+#endif

dune-solvers/transferoperators/truncatedcompressedmgtransfer.cc

+
+
+template<class DiscFuncType, class BitVectorType, class OperatorType>
+void TruncatedCompressedMGTransfer<DiscFuncType, BitVectorType, OperatorType>::prolong(const DiscFuncType& f, DiscFuncType& t,
+                                              const BitVectorType& critical) const
+{
+    if (f.size() != this->matrix_->M())
+        DUNE_THROW(Dune::Exception, "Number of entries in the coarse grid vector is not equal "
+                   << "to the number of columns of the prolongation matrix!");
+
+    if (((unsigned int)this->matrix_->N()) != critical.size())
+        DUNE_THROW(Dune::Exception, "Number of entries in the critical is not equal "
+                   << "to the number of rows of the prolongation matrix!");
+
+    t.resize(this->matrix_->N());
+
+    typedef typename DiscFuncType::Iterator      Iterator;
+    typedef typename DiscFuncType::ConstIterator ConstIterator;
+
+    typedef typename TransferOperatorType::row_type RowType;
+    typedef typename RowType::ConstIterator ColumnIterator;
+    
+    Iterator tIt      = t.begin();
+    ConstIterator fIt = f.begin(); 
+
+    for(int rowIdx=0; rowIdx<this->matrix_->N(); rowIdx++) {
+
+        const RowType& row = (*this->matrix_)[rowIdx];
+
+        *tIt = 0.0;
+        ColumnIterator cIt    = row.begin();
+        ColumnIterator cEndIt = row.end();
+
+        for(; cIt!=cEndIt; ++cIt) {
+
+            // The following lines are a matrix-vector loop with a multiple of the
+            // identity matrix, but rows belonging to critical dofs are left out
+            for (int i=0; i<blocksize; i++) {
+
+                if (!critical[rowIdx][i])
+                    (*tIt)[i] += (*cIt)[0][0] * f[cIt.index()][i];
+
+            }
+
+        }
+        
+        ++tIt;
+    }
+
+}
+
+template<class DiscFuncType, class BitVectorType, class OperatorType>
+void TruncatedCompressedMGTransfer<DiscFuncType, BitVectorType, OperatorType>::restrict(const DiscFuncType& f, DiscFuncType& t,
+                                              const BitVectorType& critical) const
+{
+    if (f.size() != this->matrix_->N())
+        DUNE_THROW(Dune::Exception, "Fine grid vector has " << f.size() << " entries "
+                   << "but the interpolation matrix has " << this->matrix_->N() << " rows!");
+
+
+    if (((unsigned int)this->matrix_->N()) != critical.size())
+        DUNE_THROW(Dune::Exception, "Number of entries in the critical is not equal "
+                   << "to the number of rows of the prolongation matrix!");
+
+    t.resize(this->matrix_->M());
+    t = 0;
+
+    typedef typename DiscFuncType::Iterator      Iterator;
+    typedef typename DiscFuncType::ConstIterator ConstIterator;
+
+    typedef typename TransferOperatorType::row_type RowType;
+    typedef typename RowType::ConstIterator ColumnIterator;
+
+    Iterator tIt      = t.begin();
+    ConstIterator fIt = f.begin();
+    
+    for (int rowIdx=0; rowIdx<this->matrix_->N(); rowIdx++) {
+        
+        const RowType& row = (*this->matrix_)[rowIdx];
+
+        ColumnIterator cIt    = row.begin();
+        ColumnIterator cEndIt = row.end();
+        
+        for(; cIt!=cEndIt; ++cIt) {
+            
+            // The following lines are a matrix-vector loop, but rows belonging
+            // to critical dofs are left out
+            typename DiscFuncType::block_type& tEntry = t[cIt.index()];
+            for (int i=0; i<blocksize; i++) {
+                
+                if (!critical[rowIdx][i])
+                    tEntry[i] += (*cIt)[0][0] * f[rowIdx][i];
+
+            }
+
+        }
+
+    }
+ 
+}
+
+
+template<class DiscFuncType, class BitVectorType, class OperatorType>
+void TruncatedCompressedMGTransfer<DiscFuncType, BitVectorType, OperatorType>::
+galerkinRestrict(const OperatorType& fineMat, OperatorType& coarseMat,
+                 const BitVectorType& critical) const
+{
+
+    if (this->recompute_ != NULL && this->recompute_->size() != (unsigned int)this->matrix_->M())
+        DUNE_THROW(Dune::Exception, "The size of the recompute_-bitfield doesn't match the "
+                   << "size of the coarse grid space!");
+
+  // ////////////////////////
+  // Nonsymmetric case
+  // ////////////////////////
+  typedef typename OperatorType::row_type RowType;
+  typedef typename RowType::Iterator ColumnIterator;
+  typedef typename RowType::ConstIterator ConstColumnIterator;
+
+  typedef typename TransferOperatorType::row_type::ConstIterator TransferConstColumnIterator;
+
+  // Clear coarse matrix
+  if (this->recompute_ == NULL)
+      coarseMat = 0;
+  else {
+      for (int i=0; i<coarseMat.N(); i++) {
+
+          RowType& row = coarseMat[i];
+          
+          // Loop over all columns of the stiffness matrix
+          ColumnIterator m    = row.begin();
+          ColumnIterator mEnd = row.end();
+          
+          for (; m!=mEnd; ++m) {
+              
+              if ((*this->recompute_)[i][0] || (*this->recompute_)[m.index()][0])
+                  *m = 0;
+
+          }
+
+      }
+
+  }
+
+  // Loop over all rows of the stiffness matrix
+  for (int v=0; v<fineMat.N(); v++) {
+      
+      const RowType& row = fineMat[v];
+
+      // Loop over all columns of the stiffness matrix
+      ConstColumnIterator m    = row.begin();
+      ConstColumnIterator mEnd = row.end();
+        
+      for (; m!=mEnd; ++m) {
+
+          int w = m.index();
+
+          // Loop over all coarse grid vectors iv that have v in their support
+          TransferConstColumnIterator im    = (*this->matrix_)[v].begin();
+          TransferConstColumnIterator imEnd = (*this->matrix_)[v].end();
+          for (; im!=imEnd; ++im) {
+
+              int iv = im.index();
+
+              // Loop over all coarse grid vectors jv that have w in their support
+              TransferConstColumnIterator jm    = (*this->matrix_)[w].begin();
+              TransferConstColumnIterator jmEnd = (*this->matrix_)[w].end();
+              
+              for (; jm!=jmEnd; ++jm) {
+                  
+                  int jv = jm.index();
+
+                  if (this->recompute_ && !((*this->recompute_)[iv][0]) && !((*this->recompute_)[jv][0]))
+                      continue;
+                  
+                  typename OperatorType::block_type& cm = coarseMat[iv][jv];
+                  
+                  // Compute im * m * jm, but omitting the critical entries
+                  for (int i=0; i<blocksize; i++) {
+                      
+                      for (int j=0; j<blocksize; j++) {
+                          
+                          // Truncated Multigrid:  Omit coupling if at least
+                          // one of the two vectors is critical
+                          if (!critical[v][i] && !critical[w][j])
+                              cm[i][j] += (*im)[0][0] * (*m)[i][j] * (*jm)[0][0];
+                          
+                      }
+                      
+                  }
+
+              }
+          }
+      }
+  }
+
+}

dune-solvers/transferoperators/truncatedcompressedmgtransfer.hh

+#ifndef TRUNCATED_COMPRESSED_MG_TRANSFER_HH
+#define TRUNCATED_COMPRESSED_MG_TRANSFER_HH
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include "truncatedmgtransfer.hh"
+#include "compressedmultigridtransfer.hh"
+
+
+/** \brief Restriction and prolongation operator for truncated multigrid
+ *
+ * This class provides prolongation and restriction operators for truncated
+ * multigrid.  That means that you can explicitly switch off certain
+ * fine grid degrees-of-freedom.  This often leads to better coarse
+ * grid corrections when treating obstacle problems.  For an introduction
+ * to the theory of truncated multigrid see 'Adaptive Monotone Multigrid
+ * Methods for Nonlinear Variational Problems' by R. Kornhuber.
+ *
+ * Currently only works for first-order Lagrangian elements!
+ */
+
+template<
+    class VectorType,
+    class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension>,
+    class MatrixType = Dune::BCRSMatrix< typename Dune::FieldMatrix<
+        typename VectorType::field_type, VectorType::block_type::size, VectorType::block_type::size> > >
+class TruncatedCompressedMGTransfer :
+    public CompressedMultigridTransfer<VectorType, BitVectorType, MatrixType>,
+    public TruncatedMGTransfer<VectorType, BitVectorType, MatrixType>
+{
+
+    enum {blocksize = VectorType::block_type::size};
+
+    typedef typename VectorType::field_type field_type;
+
+public:
+
+    typedef typename CompressedMultigridTransfer<VectorType, BitVectorType, MatrixType>::TransferOperatorType TransferOperatorType;
+    
+
+    /** \brief Default constructor */
+    TruncatedCompressedMGTransfer()
+    {}
+
+    /** \brief Restrict level fL of f and store the result in level cL of t
+     *
+     * \param critical Has to contain an entry for each degree of freedom.
+     *        Those dofs with a set bit are treated as critical.
+     */
+    void restrict(const VectorType& f, VectorType &t, const BitVectorType& critical) const;
+
+    /** \brief Restriction of  MultiGridTransfer*/
+    using CompressedMultigridTransfer< VectorType, BitVectorType, MatrixType >::restrict;
+
+    /** \brief Prolong level cL of f and store the result in level fL of t
+     * 
+     * \param critical Has to contain an entry for each degree of freedom.
+     *        Those dofs with a set bit are treated as critical.
+     */
+    void prolong(const VectorType& f, VectorType &t, const BitVectorType& critical) const;
+
+    /** \brief Prolongation of  MultiGridTransfer*/
+    using CompressedMultigridTransfer< VectorType, BitVectorType, MatrixType >::prolong;
+
+    /** \brief Galerkin assemble a coarse stiffness matrix
+     *
+     * \param critical Has to contain an entry for each degree of freedom.
+     *        Those dofs with a set bit are treated as critical.
+     */
+    void galerkinRestrict(const MatrixType& fineMat, MatrixType& coarseMat, const BitVectorType& critical) const;
+
+    /** \brief Galerkin restriction of  MultiGridTransfer*/
+    using CompressedMultigridTransfer< VectorType, BitVectorType, MatrixType >::galerkinRestrict;
+
+};
+
+
+#include "truncatedcompressedmgtransfer.cc"
+
+#endif

dune-solvers/transferoperators/truncateddensemgtransfer.cc

+
+
+template<class DiscFuncType, class BitVectorType, class OperatorType>
+void TruncatedDenseMGTransfer<DiscFuncType, BitVectorType, OperatorType>::prolong(const DiscFuncType& f, DiscFuncType& t,
+                                              const BitVectorType& critical) const
+{
+    if (f.size() != this->matrix_.M())
+        DUNE_THROW(Dune::Exception, "Number of entries in the coarse grid vector is not equal "
+                   << "to the number of columns of the prolongation matrix!");
+
+    if (((unsigned int)this->matrix_.N()) != critical.size())
+        DUNE_THROW(Dune::Exception, "Number of entries in the critical is not equal "
+                   << "to the number of rows of the prolongation matrix!");
+
+    t.resize(this->matrix_.N());
+
+    typedef typename DiscFuncType::Iterator      Iterator;
+    typedef typename DiscFuncType::ConstIterator ConstIterator;
+
+    typedef typename OperatorType::row_type RowType;
+    typedef typename RowType::ConstIterator ColumnIterator;
+    
+    Iterator tIt      = t.begin();
+    ConstIterator fIt = f.begin(); 
+
+
+
+    for(int rowIdx=0; rowIdx<this->matrix_.N(); rowIdx++) {
+
+        const RowType& row = this->matrix_[rowIdx];
+
+        *tIt = 0.0;
+        ColumnIterator cIt    = row.begin();
+        ColumnIterator cEndIt = row.end();
+
+        for(; cIt!=cEndIt; ++cIt) {
+
+            // The following lines are a matrix-vector loop, but rows belonging
+            // to critical dofs are left out
+            for (int i=0; i<blocksize; i++) {
+
+                for (int j=0; j<blocksize; j++) {
+
+                    if (!critical[rowIdx][i])
+                        (*tIt)[i] += (*cIt)[i][j] * f[cIt.index()][j];
+
+                }
+
+            }
+
+        }
+        
+        ++tIt;
+    }
+
+}
+
+template<class DiscFuncType, class BitVectorType, class OperatorType>
+void TruncatedDenseMGTransfer<DiscFuncType, BitVectorType, OperatorType>::restrict(const DiscFuncType& f, DiscFuncType& t,
+                                              const BitVectorType& critical) const
+{
+    if (f.size() != this->matrix_.N())
+        DUNE_THROW(Dune::Exception, "Fine grid vector has " << f.size() << " entries "
+                   << "but the interpolation matrix has " << this->matrix_.N() << " rows!");
+
+
+    if (((unsigned int)this->matrix_.N()) != critical.size())
+        DUNE_THROW(Dune::Exception, "Number of entries in the critical is not equal "
+                   << "to the number of rows of the prolongation matrix!");
+
+    t.resize(this->matrix_.M());
+    t = 0;
+
+    typedef typename DiscFuncType::Iterator      Iterator;
+    typedef typename DiscFuncType::ConstIterator ConstIterator;
+
+    typedef typename OperatorType::row_type RowType;
+    typedef typename RowType::ConstIterator ColumnIterator;
+
+    Iterator tIt      = t.begin();
+    ConstIterator fIt = f.begin();
+    
+    for (int rowIdx=0; rowIdx<this->matrix_.N(); rowIdx++) {
+        
+        const RowType& row = this->matrix_[rowIdx];
+
+        ColumnIterator cIt    = row.begin();
+        ColumnIterator cEndIt = row.end();
+        
+        for(; cIt!=cEndIt; ++cIt) {
+            
+            // The following lines are a matrix-vector loop, but rows belonging
+            // to critical dofs are left out
+            typename DiscFuncType::block_type& tEntry = t[cIt.index()];
+            for (int i=0; i<blocksize; i++) {
+                
+                for (int j=0; j<blocksize; j++) {
+                    
+                    if (!critical[rowIdx][j])
+                        tEntry[i] += (*cIt)[j][i] * f[rowIdx][j];
+                    
+                }
+
+            }
+
+        }
+
+    }
+ 
+}
+
+
+template<class DiscFuncType, class BitVectorType, class OperatorType>
+void TruncatedDenseMGTransfer<DiscFuncType, BitVectorType, OperatorType>::
+galerkinRestrict(const OperatorType& fineMat, OperatorType& coarseMat,
+                 const BitVectorType& critical) const
+{
+
+    if (this->recompute_ != NULL && this->recompute_->size() != (unsigned int)this->matrix_.M())
+        DUNE_THROW(Dune::Exception, "The size of the recompute_-bitfield doesn't match the "
+                   << "size of the coarse grid space!");
+
+  // ////////////////////////
+  // Nonsymmetric case
+  // ////////////////////////
+  typedef typename OperatorType::row_type RowType;
+  typedef typename RowType::Iterator ColumnIterator;
+  typedef typename RowType::ConstIterator ConstColumnIterator;
+
+  // Clear coarse matrix
+  if (this->recompute_ == NULL)
+      coarseMat = 0;
+  else {
+      for (int i=0; i<coarseMat.N(); i++) {
+
+          RowType& row = coarseMat[i];
+          
+          // Loop over all columns of the stiffness matrix
+          ColumnIterator m    = row.begin();
+          ColumnIterator mEnd = row.end();
+          
+          for (; m!=mEnd; ++m) {
+              
+              if ((*this->recompute_)[i][0] || (*this->recompute_)[m.index()][0])
+                  *m = 0;
+
+          }
+
+      }
+
+  }
+
+  // Loop over all rows of the stiffness matrix
+  for (int v=0; v<fineMat.N(); v++) {
+      
+      const RowType& row = fineMat[v];
+
+      // Loop over all columns of the stiffness matrix
+      ConstColumnIterator m    = row.begin();
+      ConstColumnIterator mEnd = row.end();
+        
+      for (; m!=mEnd; ++m) {
+
+          int w = m.index();
+
+          // Loop over all coarse grid vectors iv that have v in their support
+          ConstColumnIterator im    = this->matrix_[v].begin();
+          ConstColumnIterator imEnd = this->matrix_[v].end();
+          for (; im!=imEnd; ++im) {
+
+              int iv = im.index();
+
+              // Loop over all coarse grid vectors jv that have w in their support
+              ConstColumnIterator jm = this->matrix_[w].begin();
+              ConstColumnIterator jmEnd = this->matrix_[w].end();
+              
+              for (; jm!=jmEnd; ++jm) {
+                  
+                  int jv = jm.index();
+
+                  if (this->recompute_ && !((*this->recompute_)[iv][0]) && !((*this->recompute_)[jv][0]))
+                      continue;
+                  
+                  typename OperatorType::block_type& cm = coarseMat[iv][jv];
+                  
+                  // Compute im * m * jm, but omitting the critical entries
+                  for (int i=0; i<blocksize; i++) {
+                      
+                      for (int j=0; j<blocksize; j++) {
+                          
+                          double sum = 0.0;
+                          
+                          for (int k=0; k<blocksize; k++) {
+                              
+                              for (int l=0; l<blocksize; l++) {
+                                  // Truncated Multigrid:  Omit coupling if at least
+                                  // one of the two vectors is critical
+                                  if (!critical[v][k] && !critical[w][l]) {
+                                      sum += (*im)[k][i] * (*m)[k][l] * (*jm)[l][j];
+                                      
+                                  } 
+                                  
+                              }
+                              
+                          }
+                          
+                          cm[i][j] += sum;
+                          
+                      }
+                      
+                  }
+
+              }
+          }
+      }
+  }
+
+}

dune-solvers/transferoperators/truncateddensemgtransfer.hh

+#ifndef TRUNCATED_DENSE_MG_TRANSFER_HH
+#define TRUNCATED_DENSE_MG_TRANSFER_HH
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include "densemultigridtransfer.hh"
+#include "truncatedmgtransfer.hh"
+
+/** \brief Restriction and prolongation operator for truncated multigrid
+ *
+ * This class provides prolongation and restriction operators for truncated
+ * multigrid.  That means that you can explicitly switch off certain
+ * fine grid degrees-of-freedom.  This often leads to better coarse
+ * grid corrections when treating obstacle problems.  For an introduction
+ * to the theory of truncated multigrid see 'Adaptive Monotone Multigrid
+ * Methods for Nonlinear Variational Problems' by R. Kornhuber.
+ *
+ * Currently only works for first-order Lagrangian elements!
+ */
+
+template<
+    class VectorType,
+    class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension>,
+    class MatrixType = Dune::BCRSMatrix< typename Dune::FieldMatrix<
+        typename VectorType::field_type, VectorType::block_type::size, VectorType::block_type::size> > >
+class TruncatedDenseMGTransfer :
+    public DenseMultigridTransfer<VectorType, BitVectorType, MatrixType>,
+    public TruncatedMGTransfer<VectorType, BitVectorType, MatrixType>
+{
+
+    enum {blocksize = VectorType::block_type::size};
+
+    typedef typename VectorType::field_type field_type;
+
+public:
+
+    typedef MatrixType OperatorType;
+
+    typedef Dune::BCRSMatrix< Dune::FieldMatrix< field_type, blocksize, blocksize > > TransferOperatorType;
+
+    /** \brief Default constructor */
+    TruncatedDenseMGTransfer()
+    {}
+
+    /** \brief Restrict level fL of f and store the result in level cL of t
+     *
+     * \param critical Has to contain an entry for each degree of freedom.
+     *        Those dofs with a set bit are treated as critical.
+     */
+    void restrict(const VectorType& f, VectorType &t, const BitVectorType& critical) const;
+
+    /** \brief Restriction of  MultiGridTransfer*/
+    using DenseMultigridTransfer< VectorType, BitVectorType, MatrixType >::restrict;
+
+    /** \brief Prolong level cL of f and store the result in level fL of t
+     * 
+     * \param critical Has to contain an entry for each degree of freedom.
+     *        Those dofs with a set bit are treated as critical.
+     */
+    void prolong(const VectorType& f, VectorType &t, const BitVectorType& critical) const;
+
+    /** \brief Prolongation of  MultiGridTransfer*/
+    using DenseMultigridTransfer< VectorType, BitVectorType, MatrixType >::prolong;
+
+    /** \brief Galerkin assemble a coarse stiffness matrix
+     *
+     * \param critical Has to contain an entry for each degree of freedom.
+     *        Those dofs with a set bit are treated as critical.
+     */
+    void galerkinRestrict(const MatrixType& fineMat, MatrixType& coarseMat, const BitVectorType& critical) const;
+
+    /** \brief Galerkin restriction of  MultiGridTransfer*/
+    using DenseMultigridTransfer< VectorType, BitVectorType, MatrixType >::galerkinRestrict;
+
+};
+
+
+#include "truncateddensemgtransfer.cc"
+
+#endif

dune-solvers/transferoperators/truncatedmgtransfer.hh

+#ifndef TRUNCATED_MG_TRANSFER_HH
+#define TRUNCATED_MG_TRANSFER_HH
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/bitsetvector.hh>
+
+/** \brief Restriction and prolongation interface class for truncated multigrid
+ *
+ * This class provides prolongation and restriction operators for truncated
+ * multigrid.  That means that you can explicitly switch off certain
+ * fine grid degrees-of-freedom.  This often leads to better coarse
+ * grid corrections when treating obstacle problems.  For an introduction
+ * to the theory of truncated multigrid see 'Adaptive Monotone Multigrid
+ * Methods for Nonlinear Variational Problems' by R. Kornhuber.
+ *
+ * Currently only works for first-order Lagrangian elements!
+ */
+
+template<
+    class VectorType,
+    class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension>,
+    class MatrixType = Dune::BCRSMatrix< typename Dune::FieldMatrix<
+        typename VectorType::field_type, VectorType::block_type::size, VectorType::block_type::size> > >
+class TruncatedMGTransfer
+{
+
+    enum {blocksize = VectorType::block_type::size};
+
+    typedef typename VectorType::field_type field_type;
+
+public:
+
+    /** \brief Default constructor */
+    TruncatedMGTransfer() : recompute_(NULL)
+    {}
+
+    /** \brief Restrict level fL of f and store the result in level cL of t
+     *
+     * \param critical Has to contain an entry for each degree of freedom.
+     *        Those dofs with a set bit are treated as critical.
+     */
+    virtual void restrict(const VectorType& f, VectorType &t, const BitVectorType& critical) const = 0;
+
+    /** \brief Prolong level cL of f and store the result in level fL of t
+     * 
+     * \param critical Has to contain an entry for each degree of freedom.
+     *        Those dofs with a set bit are treated as critical.
+     */
+    virtual void prolong(const VectorType& f, VectorType &t, const BitVectorType& critical) const = 0;
+
+    /** \brief Galerkin assemble a coarse stiffness matrix
+     *
+     * \param critical Has to contain an entry for each degree of freedom.
+     *        Those dofs with a set bit are treated as critical.
+     */
+    virtual void galerkinRestrict(const MatrixType& fineMat, 
+                                  MatrixType& coarseMat, 
+                                  const BitVectorType& critical) const = 0;
+
+    /** \brief Bitfield specifying a subsets of dofs which need to be recomputed
+     * when doing Galerkin restriction
+     *
+     * If this pointer points to NULL it has no effect.  If not it is expected
+     * to point to a bitfield with the size of the coarse function space.
+     * Then, when calling galerkinRestrict(), only those matrix entries which
+     * involve at least one dof which has a set bit get recomputed.  Depending
+     * on the problem this can lead to considerable time savings.
+     */
+    const Dune::BitSetVector<1>* recompute_;
+};
+
+#endif