diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5746cb97dd0abb5a5f6084cbac13fc3b3930cabe..87bb694078c4ada30103e76ec4f04a3fd8a74126 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,4 +4,5 @@
 
 ## Deprecations and removals
 
-- ...
+- The file `blockgsstep.hh` has been removed after four years of deprecation.  Use the replacement
+  in `blockgssteps.hh` instead.
diff --git a/dune/solvers/iterationsteps/CMakeLists.txt b/dune/solvers/iterationsteps/CMakeLists.txt
index 5fde9a8fc16ff62b4d9efe5f674b834eed3e8750..935f4cda11435d061f2c86cbd322eef3dcc5a600 100644
--- a/dune/solvers/iterationsteps/CMakeLists.txt
+++ b/dune/solvers/iterationsteps/CMakeLists.txt
@@ -1,7 +1,5 @@
 install(FILES
     amgstep.hh
-    blockgsstep.cc
-    blockgsstep.hh
     blockgssteps.hh
     cgstep.cc
     cgstep.hh
diff --git a/dune/solvers/iterationsteps/blockgsstep.cc b/dune/solvers/iterationsteps/blockgsstep.cc
deleted file mode 100644
index 9092e4d07649e47491ec8386878b3a7342df6f57..0000000000000000000000000000000000000000
--- a/dune/solvers/iterationsteps/blockgsstep.cc
+++ /dev/null
@@ -1,75 +0,0 @@
-// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=8 sw=4 sts=4:
-
-#include <cassert>
-
-#include <dune/matrix-vector/axpy.hh>
-
-template<class MatrixType, class DiscFuncType, class BitVectorType>
-inline
-void BlockGSStep<MatrixType, DiscFuncType, BitVectorType>::
-residual(int index, VectorBlock& r) const
-{
-    r = (*this->rhs_)[index];
-
-    const auto& row = (*this->mat_)[index];
-    for (auto cIt = row.begin(); cIt!=row.end(); ++cIt) {
-        // r_i -= A_ij x_j
-        Dune::MatrixVector::subtractProduct(r, *cIt, (*this->x_)[cIt.index()]);
-    }
-
-}
-
-template<class MatrixType, class DiscFuncType, class BitVectorType>
-inline
-void BlockGSStep<MatrixType, DiscFuncType, BitVectorType>::iterate()
-{
-    assert(this->hasIgnore());
-
-    if (gs_type_ != Direction::BACKWARD)
-        for (std::size_t i=0; i<this->x_->size(); i++)
-            iterate_step(i);
-
-    if (gs_type_ != Direction::FORWARD)
-        for (std::size_t i=this->x_->size()-1; i>=0 && i<this->x_->size(); i--)
-            iterate_step(i);
-}
-
-template<class MatrixType, class DiscFuncType, class BitVectorType>
-inline
-void BlockGSStep<MatrixType, DiscFuncType, BitVectorType>::iterate_step(int i)
-{
-    auto const &ignore_i = this->ignore()[i];
-    if (ignore_i.all())
-        return;
-
-    VectorBlock r;
-    residual(i, r);
-    const auto& mat_ii = (*this->mat_)[i][i];
-
-    // Compute correction v = A_{i,i}^{-1} r[i]
-    VectorBlock v;
-
-    if (ignore_i.none()) {
-        // No degree of freedom shall be ignored --> solve linear problem
-        mat_ii.solve(v, r);
-    } else {
-        // Copy the matrix and adjust rhs and matrix so that dofs given in ignoreNodes[i]
-        // are not touched
-        typename MatrixType::block_type matRes;
-        for (int j = 0; j < BlockSize; ++j) {
-            if (ignore_i[j]) {
-                r[j] = 0;
-
-                for (int k = 0; k < BlockSize; ++k)
-                    matRes[j][k] = (k == j);
-            } else
-                matRes[j] = mat_ii[j];
-
-        }
-        matRes.solve(v, r);
-    }
-    // Add correction
-    VectorBlock& x = (*this->x_)[i];
-    x += v;
-}
diff --git a/dune/solvers/iterationsteps/blockgsstep.hh b/dune/solvers/iterationsteps/blockgsstep.hh
deleted file mode 100644
index ecf89f21609fd0fc89ae034c86aafbe2bc01efca..0000000000000000000000000000000000000000
--- a/dune/solvers/iterationsteps/blockgsstep.hh
+++ /dev/null
@@ -1,61 +0,0 @@
-// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=8 sw=4 sts=4:
-#ifndef DUNE_BLOCK_GAUSS_SEIDEL_STEP_HH
-#define DUNE_BLOCK_GAUSS_SEIDEL_STEP_HH
-
-#warning This header is deprecated. Please use Dune::Solvers::BlockGSStep from blockgssteps.hh instead of ::BlockGSStep. See also gssteptest.cc for usage examples.
-
-#include <dune/common/bitsetvector.hh>
-
-#include "lineariterationstep.hh"
-
-/** \brief A linear Gauß-Seidel iteration step for convex problems which have a block-structure.
- *
- *  \tparam MatrixType The linear operator type
- *  \tparam DiscFuncType The block vector type of the right hand side and the iterates
- *  \tparam BitVectorType The type of the bit-vector specifying degrees of freedom that shall be ignored.
- *
- */
-template<class MatrixType,
-         class DiscFuncType,
-         class BitVectorType = Dune::BitSetVector<DiscFuncType::block_type::dimension> >
-         class BlockGSStep : public LinearIterationStep<MatrixType, DiscFuncType, BitVectorType>
-    {
-
-        using VectorBlock = typename DiscFuncType::block_type;
-
-        enum {BlockSize = VectorBlock::dimension};
-
-    public:
-        enum Direction { FORWARD, BACKWARD, SYMMETRIC };
-
-        //! Default constructor.
-        BlockGSStep(Direction gs_type = FORWARD)
-            : gs_type_(gs_type)
-        {}
-
-        //! Constructor with a linear problem
-        BlockGSStep(const MatrixType& mat, DiscFuncType& x, const DiscFuncType& rhs, Direction gs_type = FORWARD)
-            : LinearIterationStep<MatrixType,DiscFuncType>(mat, x, rhs), gs_type_(gs_type)
-        {}
-
-        //! Perform one iteration
-        virtual void iterate();
-
-        /** \brief Compute the residual of the current iterate of a (block) degree of freedom
-         *
-         *  \param index Global index of the (block) degree of freedom
-         *  \param r Write residual in this vector
-         */
-        virtual void residual(int index, VectorBlock& r) const;
-
-    protected:
-        const Direction gs_type_;
-
-        void iterate_step(int i);
-    };
-
-
-#include "blockgsstep.cc"
-
-#endif
diff --git a/dune/solvers/iterationsteps/truncatedblockgsstep.hh b/dune/solvers/iterationsteps/truncatedblockgsstep.hh
index 86fa738fdce830c81a9e8bf8c6479ee0f46456f8..57e2863b84935620d08d759631cd1253ff031c62 100644
--- a/dune/solvers/iterationsteps/truncatedblockgsstep.hh
+++ b/dune/solvers/iterationsteps/truncatedblockgsstep.hh
@@ -42,7 +42,7 @@ public:
     //! Perform one iteration
     virtual void iterate()
     {
-        TruncatedBlockGSStepNS::RecursiveGSStep<MatrixType::blocklevel>::apply(*this->mat_, *this->rhs_, this->ignore(), *this->x_, innerLoops_);
+        TruncatedBlockGSStepNS::RecursiveGSStep<Dune::blockLevel<MatrixType>()>::apply(*this->mat_, *this->rhs_, this->ignore(), *this->x_, innerLoops_);
     }
 
 protected:
diff --git a/dune/solvers/test/common.hh b/dune/solvers/test/common.hh
index 028dad406088381cb1b771594812bcf02df47526..7015536debd45974e770bb404e2b8fed4f9c5275 100644
--- a/dune/solvers/test/common.hh
+++ b/dune/solvers/test/common.hh
@@ -7,7 +7,6 @@
 
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
-#include <dune/common/function.hh>
 
 #include <dune/istl/matrixindexset.hh>
 
@@ -30,23 +29,6 @@
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/norms/twonorm.hh>
 
-template <class DomainType, class RangeType>
-class ConstantFunction :
-    public Dune::VirtualFunction<DomainType, RangeType>
-{
-    public:
-        ConstantFunction(double c):
-            c_(c)
-        {}
-
-        void evaluate(const DomainType& x, RangeType& y) const
-        {
-            y = c_;
-        }
-
-    private:
-        double c_;
-};
 
 template<class GridView, class Matrix>
 void constructPQ1Pattern(const GridView& gridView, Matrix& matrix)
@@ -231,7 +213,6 @@ void assemblePQ1RHS(const GridView& gridView, Vector& r, const Function& f)
     typedef typename FiniteElementCache::FiniteElementType FiniteElement;
 
     typedef typename FiniteElement::Traits::LocalBasisType::Traits::RangeType RangeType;
-    typedef typename Function::RangeType FunctionRangeType;
 
     const auto& indexSet = gridView.indexSet();
     FiniteElementCache cache;
@@ -265,8 +246,7 @@ void assemblePQ1RHS(const GridView& gridView, Vector& r, const Function& f)
             fe.localBasis().evaluateFunction(quadPos, values);
 
             // evaluate function
-            FunctionRangeType fAtPos;
-            f.evaluate(geometry.global(quadPos), fAtPos);
+            auto fAtPos = f(geometry.global(quadPos));
 
             // add vector entries
             double z = pt.weight() * integrationElement;
diff --git a/dune/solvers/test/gssteptest.cc b/dune/solvers/test/gssteptest.cc
index 5538592858ac73e4c249c2d5bc00cc0ad2147936..424ea189a240bf2e010a980a122aaa220fb285b4 100644
--- a/dune/solvers/test/gssteptest.cc
+++ b/dune/solvers/test/gssteptest.cc
@@ -10,7 +10,6 @@
 #include <dune/common/stringutility.hh>
 
 #include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/iterationsteps/blockgsstep.hh>
 #include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
 #include <dune/solvers/iterationsteps/blockgssteps.hh>
 #include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
@@ -102,14 +101,9 @@ struct GSTestSuite {
     };
 
     {
-      BlockGSStep<Matrix, Vector> gsStep;
-      test(&gsStep, "BlockGS");
-    }
-    {
-      BlockGSStep<Matrix, Vector> gsStep1;
       auto gsStep2 = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(
           Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
-      diffTest(&gsStep1, "BlockGS", &gsStep2, "Dune::Solvers::BlockGS(Direct)");
+      test(&gsStep2, "Dune::Solvers::BlockGS(Direct)");
     }
     {
       TruncatedBlockGSStep<Matrix, Vector> gsStep1(1);
diff --git a/dune/solvers/test/mmgtest.cc b/dune/solvers/test/mmgtest.cc
index ae5268f0dd7872aa2a3eebb326bee4ef635f23c4..80318b5d6927bf6dfe203a039d141cb4fb8c255a 100644
--- a/dune/solvers/test/mmgtest.cc
+++ b/dune/solvers/test/mmgtest.cc
@@ -123,7 +123,9 @@ bool checkWithGrid(const GridType& grid, const std::string fileName="")
 
     Vector rhs(A.N());
     rhs = 0;
-    ConstantFunction<DomainType, RangeType> f(50);
+    auto f = [](const DomainType& x) -> RangeType
+      {return 50.0; };
+
     assemblePQ1RHS(gridView, rhs, f);
 
     Vector u(A.N());
diff --git a/dune/solvers/test/obstacletnnmgtest.cc b/dune/solvers/test/obstacletnnmgtest.cc
index 23858a85b5d310297d36d13fc1a9a1d043c7d16c..f6b938b2dd2eab6e8bc76b4020b24dcd45fbaa9b 100644
--- a/dune/solvers/test/obstacletnnmgtest.cc
+++ b/dune/solvers/test/obstacletnnmgtest.cc
@@ -114,7 +114,9 @@ bool checkWithGrid(const GridType& grid, const std::string fileName="")
 
     Vector rhs(A.N());
     rhs = 0;
-    ConstantFunction<DomainType, RangeType> f(50);
+    auto f = [](const DomainType& x) -> RangeType
+      {return 50.0; };
+
     assemblePQ1RHS(gridView, rhs, f);
 
     Vector u(A.N());
diff --git a/dune/solvers/test/projectedgradienttest.cc b/dune/solvers/test/projectedgradienttest.cc
index 131c73b0232efaa0804f7ded640769eddf540f4f..4e76bd441d923cc501452b13e71f5dc574ff55e0 100644
--- a/dune/solvers/test/projectedgradienttest.cc
+++ b/dune/solvers/test/projectedgradienttest.cc
@@ -84,7 +84,9 @@ bool checkWithGrid(const GridType& grid, const std::string fileName, int maxIter
 
     Vector rhs(A.N());
     rhs = 0;
-    ConstantFunction<DomainType, RangeType> f(50);
+    auto f = [](const DomainType& x) -> RangeType
+      {return 50.0; };
+
     assemblePQ1RHS(gridView, rhs, f);
 
     Vector u(A.N());
diff --git a/dune/solvers/test/quadraticipoptsolvertest.cc b/dune/solvers/test/quadraticipoptsolvertest.cc
index 52c177aab04d3484a8d47898b9fd7702259aa6b8..4ea4f561f480438188f1fa1e0c0a9ca8b292e1c1 100644
--- a/dune/solvers/test/quadraticipoptsolvertest.cc
+++ b/dune/solvers/test/quadraticipoptsolvertest.cc
@@ -77,7 +77,9 @@ bool checkWithGrid(const GridType& grid, const std::string fileName="")
 
     Vector rhs(A.N());
     rhs = 0;
-    ConstantFunction<DomainType, RangeType> f(50);
+    auto f = [](const DomainType& x) -> RangeType
+      {return 50.0; };
+
     assemblePQ1RHS(gridView, rhs, f);
 
     Vector u(A.N());