From d4dd632072981ba32c8bf7aded97ea4675a39ac2 Mon Sep 17 00:00:00 2001
From: Uli Sack <usack@math.fu-berlin.de>
Date: Wed, 25 Jan 2012 09:45:09 +0000
Subject: [PATCH] added specialization for MatrixType=SumOperator<SparseMatrix,
LowRankMatrix>
[[Imported from SVN: r5474]]
---
.../compressedmultigridtransfer.hh | 154 ++++++++++++++++++
1 file changed, 154 insertions(+)
diff --git a/dune/solvers/transferoperators/compressedmultigridtransfer.hh b/dune/solvers/transferoperators/compressedmultigridtransfer.hh
index 37b3ea74..cce655cb 100644
--- a/dune/solvers/transferoperators/compressedmultigridtransfer.hh
+++ b/dune/solvers/transferoperators/compressedmultigridtransfer.hh
@@ -5,6 +5,9 @@
#include <dune/common/fmatrix.hh>
#include <dune/common/bitsetvector.hh>
+#include <dune/solvers/operators/sumoperator.hh>
+#include <dune/solvers/common/staticmatrixtools.hh>
+
#include "genericmultigridtransfer.hh"
#include "multigridtransfer.hh"
@@ -142,6 +145,157 @@ public:
+protected:
+
+ TransferOperatorType* matrix_;
+ bool matrixInternallyAllocated;
+
+};
+
+template<
+ class VectorType,
+ class BitVectorType,
+ class SparseMatrixType,
+ class LowRankMatrixType >
+class CompressedMultigridTransfer<VectorType, BitVectorType, SumOperator<SparseMatrixType, LowRankMatrixType> > :
+ public MultigridTransfer<VectorType, BitVectorType, SumOperator<SparseMatrixType,LowRankMatrixType> >
+{
+
+ enum {blocksize = VectorType::block_type::dimension};
+
+ typedef typename VectorType::field_type field_type;
+
+public:
+
+ typedef SumOperator<SparseMatrixType, LowRankMatrixType> 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.sparseMatrix(), coarseMat.sparseMatrix(), -1);
+
+ const typename LowRankMatrixType::LowRankFactorType& lr_fine_factor(fineMat.lowRankMatrix().lowRankFactor());
+ typename LowRankMatrixType::LowRankFactorType& lr_coarse_factor(coarseMat.lowRankMatrix().lowRankFactor());
+
+ lr_coarse_factor = 0;
+
+ for (size_t i=0; i<matrix_->N(); ++i)
+ {
+ typename TransferOperatorType::ColIterator j_it = (*matrix_)[i].begin();
+ typename TransferOperatorType::ColIterator j_end = (*matrix_)[i].end();
+
+ for (; j_it!=j_end; ++j_it)
+ {
+ size_t j = j_it.index();
+
+ for (size_t k=0; k<lr_fine_factor.N(); ++k)
+ StaticMatrix::axpy(lr_coarse_factor[k][j], *j_it, lr_fine_factor[k][i]);
+ }
+ }
+ }
+
+
+ /** \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.sparseMatrix(), coarseMat.sparseMatrix(), -1);
+
+ coarseMat.lowRankMatrix().lowRankFactor().setSize(fineMat.lowRankMatrix().lowRankFactor().N(), matrix_->M());
+ }
+
+
+ /** \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_;
--
GitLab