diff --git a/dune/solvers/iterationsteps/projectedlinegsstep.cc b/dune/solvers/iterationsteps/projectedlinegsstep.cc
new file mode 100755
index 0000000000000000000000000000000000000000..91c2e930006e4ba35724cb8db3895da869ecf67c
--- /dev/null
+++ b/dune/solvers/iterationsteps/projectedlinegsstep.cc
@@ -0,0 +1,231 @@
+#include <dune/istl/scaledidmatrix.hh>
+
+#include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
+
+template<class MatrixType, class VectorType, class BitVectorType>
+void ProjectedLineGSStep<MatrixType, VectorType, BitVectorType >::
+solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
+ VectorType& x,
+ const VectorType& rhs,
+ const std::vector<BoxConstraint<field_type,BlockSize> >& localObstacles) const
+{
+ typedef Dune::BTDMatrix<typename MatrixType::block_type> LocalMatrixType;
+
+ // ///////////////////////////
+ // Set initial iterate
+ // Make rhs the initial iterate: it already contains the correct Dirichlet values.
+ // The rest doesn't matter
+ // ///////////////////////////
+ x = rhs;
+
+
+ // //////////////////////////////////////
+ // The TNNMG loop
+ // //////////////////////////////////////
+
+ // stupid:
+ Dune::BitSetVector<1> hasObstacle(rhs.size(), true);
+ // The nodes to ignore are already written into the matrix/rhs
+ Dune::BitSetVector<1> ignoreNodes(rhs.size(), false);
+
+ do {
+
+ // ///////////////////////////
+ // Nonlinear presmoothing
+ // ///////////////////////////
+ VectorType rhsCopy = rhs; // need a non-const version of rhs
+
+ ProjectedBlockGSStep<LocalMatrixType, VectorType> nonlinearSmootherStep(matrix, x, rhsCopy);
+
+ nonlinearSmootherStep.ignoreNodes_ = &ignoreNodes;
+ nonlinearSmootherStep.hasObstacle_ = &hasObstacle;
+ nonlinearSmootherStep.obstacles_ = &localObstacles;
+
+ for (int i=0; i<3; i++)
+ nonlinearSmootherStep.iterate();
+
+ // Compute residual
+ VectorType residual = rhs;
+ matrix.mmv(x,residual);
+
+ // ///////////////////////////
+ // Truncation
+ // ///////////////////////////
+
+ // ///////////////////////////
+ // Linear correction
+ // ///////////////////////////
+
+ VectorType corr(residual.size());
+ matrix.solve(corr, residual); // the direct linear solution algorithm
+
+ // //////////////////////////////////////////////////////////
+ // Line search in the direction of the projected coarse
+ // grid correction to ensure monotonicity.
+ // //////////////////////////////////////////////////////////
+
+ // L2-projection of the correction onto the defect obstacle
+ for (int i=0; i<corr.size(); i++) {
+
+ for (int j=0; j<BlockSize; j++) {
+
+ corr[i][j] = std::max(corr[i][j], localObstacles[i].lower(j) - x[i][j]);
+ corr[i][j] = std::min(corr[i][j], localObstacles[i].upper(j) - x[i][j]);
+
+ }
+
+ }
+
+ // Construct obstacles in the direction of the projected correction
+ BoxConstraint<field_type,1> lineSearchObs;
+ for (int i=0; i<corr.size(); i++) {
+
+ for (int j=0; j<BlockSize; j++) {
+
+ if (corr[i][j] > 0) {
+
+ // This division can cause nan on some platforms...
+ if (!isnan( (localObstacles[i].lower(j)-x[i][j]) / corr[i][j]) )
+ lineSearchObs.lower(0) = std::max(lineSearchObs.lower(0),
+ (localObstacles[i].lower(j)-x[i][j]) / corr[i][j]);
+
+ if (!isnan( (localObstacles[i].upper(j)-x[i][j]) / corr[i][j]) )
+ lineSearchObs.upper(0) = std::min(lineSearchObs.upper(0),
+ (localObstacles[i].upper(j)-x[i][j]) / corr[i][j]);
+ }
+
+ if (corr[i][j] < 0) {
+ if (!isnan( (localObstacles[i].upper(j)-x[i][j]) / corr[i][j]) )
+ lineSearchObs.lower(0) = std::max(lineSearchObs.lower(0),
+ (localObstacles[i].upper(j)-x[i][j]) / corr[i][j]);
+ if (!isnan( (localObstacles[i].lower(j)-x[i][j]) / corr[i][j]) )
+ lineSearchObs.upper(0) = std::min(lineSearchObs.upper(0),
+ (localObstacles[i].lower(j)-x[i][j]) / corr[i][j]);
+ }
+
+ }
+
+ }
+
+ // abort when the linear correction has a small norm
+ field_type correctionNormSquared = EnergyNorm<LocalMatrixType,VectorType>::normSquared(corr,matrix);
+
+ // If the correction is very small we may leave the loop here
+ // A relative criterion may be more appropriate here
+ if (correctionNormSquared < 1e-10)
+ break;
+
+ // Line search
+ field_type alpha = (residual*corr) / correctionNormSquared;
+ alpha = std::max(std::min(alpha, lineSearchObs.upper(0)), lineSearchObs.lower(0));
+
+ // Even if the correction has a large norm the scaled norm may be small.
+ // A relative criterion may be more appropriate here
+ if (alpha*correctionNormSquared < 1e-10)
+ break;
+
+ // add scaled correction
+ x.axpy(alpha,corr);
+
+ } while (true);
+
+}
+
+
+
+template<class MatrixType, class VectorType, class BitVectorType>
+void ProjectedLineGSStep<MatrixType, VectorType, BitVectorType>::iterate()
+{
+
+ // input of this method: x^(k) (not the permuted version of x^(k)!)
+
+ const MatrixType& matrix = *this->mat_;
+
+ int number_of_blocks = this->blockStructure_.size();
+
+ // iterate over the blocks
+ for (int b_num=0; b_num < number_of_blocks; b_num++ ) {
+
+ const int current_block_size = this->blockStructure_[b_num].size();
+
+ //! compute and save the residuals for the curent block:
+ // determine the (permuted) residuals r[p(i)],..., r[p(i+current_block_size-1)]
+ // this means that we determine the residuals for the current block
+ VectorType permuted_r_i(current_block_size);
+ for (int k=0; k<current_block_size; k++)
+ {
+ if ( (*this->ignoreNodes_)[this->blockStructure_[b_num][k]][0])
+ permuted_r_i[k] = 0.0;
+ else
+ residual( this->blockStructure_[b_num][k], permuted_r_i[k]); // get r[p(i+k)]
+ }
+ // permuted_r_i[k] is the residual for r[p(i+k)], which means the residual \tilde{r}[i+k]
+ // note: calling the method 'residual' with the permuted index implies that there is no additional permutation required within 'residual'
+
+ // permuted_v_i is for saving the correction for x[p(i)],..., x[p(i+current_block_size-1)]:
+ VectorType permuted_v_i(current_block_size); // permuted_v_i[k] = v[p(i+k)]
+ // (Later: x_now[p(i+k)] = x_now[p(i+k)] + v[p(i+k)] )
+
+ // //////////////////////////////////////////////////////////////////////////////////////////////////
+ // Copy the linear system for the current line/block into a tridiagonal matrix
+ // //////////////////////////////////////////////////////////////////////////////////////////////////
+
+ Dune::BTDMatrix<typename MatrixType::block_type> tridiagonalMatrix(current_block_size);
+
+ for (int j=0; j<current_block_size; j++) {
+
+ if ( (*this->ignoreNodes_)[this->blockStructure_[b_num][j]][0] ) {
+
+ // left off-diagonal:
+ if (j>0)
+ tridiagonalMatrix[j][j-1] = 0;
+
+ // diagonal
+ tridiagonalMatrix[j][j] = Dune::ScaledIdentityMatrix<field_type,BlockSize>(1);
+
+ // left off-diagonal:
+ if (j<current_block_size-1)
+ tridiagonalMatrix[j][j+1] = 0;
+
+ } else {
+
+ // left off-diagonal:
+ if (j>0)
+ tridiagonalMatrix[j][j-1] = matrix[this->blockStructure_[b_num][j]][this->blockStructure_[b_num][j-1]];
+
+ // diagonal
+ tridiagonalMatrix[j][j] = matrix[this->blockStructure_[b_num][j]][this->blockStructure_[b_num][j]];
+
+ // left off-diagonal:
+ if (j<current_block_size-1)
+ tridiagonalMatrix[j][j+1] = matrix[this->blockStructure_[b_num][j]][this->blockStructure_[b_num][j+1]];
+
+ }
+
+ }
+
+ // ///////////////////////////////////////////////////////////////////
+ // Get the local set of obstacles for this block
+ // ///////////////////////////////////////////////////////////////////
+
+ std::vector<BoxConstraint<field_type,BlockSize> > localDefectConstraints(current_block_size);
+
+ for (int j=0; j<current_block_size; j++) {
+ localDefectConstraints[j] = (*obstacles_)[this->blockStructure_[b_num][j]];
+ localDefectConstraints[j] -= (*this->x_)[this->blockStructure_[b_num][j]];
+ }
+
+ // ////////////////////////////////////////
+ // Solve the tridiagonal system
+ // ////////////////////////////////////////
+ this->solveLocalSystem(tridiagonalMatrix, permuted_v_i, permuted_r_i, localDefectConstraints);
+
+ // //////////////////////////////////////////////////////////////////////
+ // Add the correction to the current iterate
+ // //////////////////////////////////////////////////////////////////////
+ for (int k=0; k<current_block_size; k++)
+ (*this->x_)[this->blockStructure_[b_num][k]] += permuted_v_i[k];
+
+ }
+
+}
diff --git a/dune/solvers/iterationsteps/projectedlinegsstep.hh b/dune/solvers/iterationsteps/projectedlinegsstep.hh
new file mode 100755
index 0000000000000000000000000000000000000000..2cd0c74bd8f91e3a00586ae0d9947c3ecf229f68
--- /dev/null
+++ b/dune/solvers/iterationsteps/projectedlinegsstep.hh
@@ -0,0 +1,49 @@
+#ifndef DUNE_PROJECTED_LINE_GAUSS_SEIDEL_STEP_HH
+#define DUNE_PROJECTED_LINE_GAUSS_SEIDEL_STEP_HH
+
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/istl/btdmatrix.hh>
+
+#include <dune/solvers/boxconstraint.hh>
+#include <dune/solvers/iterationsteps/linegsstep.hh>
+
+template<class MatrixType,
+ class VectorType,
+ class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> >
+class ProjectedLineGSStep : public LineGSStep<MatrixType, VectorType, BitVectorType>
+{
+
+ typedef typename VectorType::block_type VectorBlock;
+
+ typedef typename VectorType::field_type field_type;
+
+ enum {BlockSize = VectorBlock::dimension};
+
+public:
+
+ //! Default constructor. Doesn't init anything
+ ProjectedLineGSStep() {}
+
+ //! Constructor for usage of Permutation Manager
+ ProjectedLineGSStep( const std::vector<std::vector<unsigned int> >& blockStructure)
+ : LineGSStep<MatrixType, VectorType, BitVectorType>( blockStructure )
+ {}
+
+ /** \brief Solve one tridiagonal system */
+ void solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
+ VectorType& x,
+ const VectorType& rhs,
+ const std::vector<BoxConstraint<field_type,BlockSize> >& localObstacles) const;
+
+ //! Perform one iteration
+ virtual void iterate();
+
+ const std::vector<BoxConstraint<field_type,BlockSize> >* obstacles_;
+
+};
+
+
+#include "projectedlinegsstep.cc"
+
+#endif