moved to new module dune-solvers

[[Imported from SVN: r2345]]

dune-solvers/iterationsteps/projectedblockgsstep.cc

+template<class OperatorType, class DiscFuncType>
+inline
+DiscFuncType ProjectedBlockGSStep<OperatorType, DiscFuncType>::getSol()
+{
+    return *(this->x_);
+}
+
+template<class OperatorType, class DiscFuncType>
+inline
+void ProjectedBlockGSStep<OperatorType, DiscFuncType>::iterate()
+{
+    if (hasObstacle_->size()!= (unsigned int)this->x_->size())
+        DUNE_THROW(SolverError, "Size of hasObstacle (" << hasObstacle_->size() 
+                   << ") doesn't match solution vector (" << this->x_->size() << ")");
+
+    const OperatorType& mat = *this->mat_;
+
+    for (int i=0; i<this->x_->size(); i++) {
+
+        /** \todo Handle more general boundary conditions */
+        if ((*this->ignoreNodes_)[i][0])
+            continue;
+
+        bool zeroDiagonal = false;
+        for (int j=0; j<BlockSize; j++) {
+            // When using this solver as part of a truncated multigrid solver,
+            // the diagonal entries of the matrix may get completely truncated
+            // away.  In this case we just do nothing here.
+            zeroDiagonal |= (mat[i][i][j][j] < 1e-10);
+        }
+
+        if (zeroDiagonal)
+            continue;
+
+        VectorBlock r;
+        residual(i, r);
+
+        // Compute x_i += A_{i,i}^{-1} r[i] 
+        VectorBlock v;
+        VectorBlock& x = (*this->x_)[i];
+
+
+        if ((*hasObstacle_)[i] == false) {
+
+            // No obstacle --> solve linear problem
+            mat[i][i].solve(v, r);
+        
+        } else {
+            
+            // Solve the local constraint minimization problem
+            // We use a projected Gauss-Seidel, for lack of anything better
+            BoxConstraint<field_type,BlockSize> defectObstacle = (*obstacles_)[i];
+            defectObstacle -= x;
+
+            // Initial correction
+            v = 0;
+
+            for (int j=0; j<20; j++) {
+
+                for (int k=0; k<BlockSize; k++) {
+
+                    // Compute residual
+                    field_type sr = 0;
+                    for (int l=0; l<BlockSize; l++)
+                        sr += mat[i][i][k][l] * v[l];
+
+                    sr = r[k] - sr;
+                    v[k] += sr / mat[i][i][k][k];
+
+                    // Project
+                    if (v[k] < defectObstacle.lower(k))
+                        v[k] = defectObstacle.lower(k);
+                    else if (v[k] > defectObstacle.upper(k))
+                        v[k] = defectObstacle.upper(k);
+
+                }
+
+            }
+
+        }
+
+        // Add correction
+        x += v;
+
+    }
+
+}

dune-solvers/iterationsteps/projectedblockgsstep.hh

+#ifndef DUNE_PROJECTED_BLOCK_GAUSS_SEIDEL_STEP_HH
+#define DUNE_PROJECTED_BLOCK_GAUSS_SEIDEL_STEP_HH
+
+#include <vector>
+
+#include "blockgsstep.hh"
+#include "boxconstraint.hh"
+
+template<class OperatorType, class DiscFuncType>
+class ProjectedBlockGSStep : public BlockGSStep<OperatorType, DiscFuncType>
+{
+    
+    typedef typename DiscFuncType::block_type VectorBlock;
+    
+    typedef typename DiscFuncType::field_type field_type;
+
+    enum {BlockSize = VectorBlock::dimension};
+    
+public:
+    
+    //! Default constructor.  Doesn't init anything
+    ProjectedBlockGSStep() {}
+    
+    //! Constructor with a linear problem
+    ProjectedBlockGSStep(OperatorType& mat, DiscFuncType& x, DiscFuncType& rhs)
+        : LinearIterationStep<OperatorType,DiscFuncType>(mat, x, rhs)
+    {}
+    
+    //! Perform one iteration
+    virtual void iterate();
+    
+    virtual DiscFuncType getSol();
+    
+    Dune::BitSetVector<1>* hasObstacle_;
+    
+    const std::vector<BoxConstraint<field_type,BlockSize> >* obstacles_;
+};
+
+#include "projectedblockgsstep.cc"
+
+#endif