A line smoother Gauss-Seidel that respects box constraints. It works by

solving each local system with a TNNMG using a direct solver for the linear
correction problems.

Not completely debugged yet.  Stay tuned.

[[Imported from SVN: r3323]]

dune/solvers/iterationsteps/projectedlinegsstep.cc

+#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];
+
+    }
+
+}

dune/solvers/iterationsteps/projectedlinegsstep.hh

+#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