moved to new module dune-solvers

[[Imported from SVN: r2359]]

dune-solvers/iterationsteps/mmgstep.cc

+#include <dune/ag-common/transferoperators/truncatedmgtransfer.hh>
+#include <dune/ag-common/projectedblockgsstep.hh>
+#include <dune/ag-common/computeenergy.hh>
+
+
+#ifdef HAVE_IPOPT
+#include <dune/ag-common/quadraticipopt.hh>
+#endif
+
+template <class OperatorType, class DiscFuncType>
+void MonotoneMGStep<OperatorType, DiscFuncType>::
+preprocess()
+{
+    // Call preprocess of the base class
+    MultigridStep<OperatorType,DiscFuncType>::preprocess();
+
+    // Then specify the subset of entries to be reassembled at each iteration
+    recompute_.resize(this->numLevels_);
+    recompute_[recompute_.size()-1] = (*hasObstacle_)[recompute_.size()-1];
+    for (int i=this->mgTransfer_.size()-1; i>=0; i--)
+        this->mgTransfer_[i]->restrictScalarBitField(recompute_[i+1], recompute_[i]);
+        
+
+    for (size_t i=0; i<this->mgTransfer_.size(); i++)
+        dynamic_cast<TruncatedMGTransfer<DiscFuncType>*>(this->mgTransfer_[i])->recompute_ = &recompute_[i];
+    
+    // /////////////////////////////////////////////
+    //   Set up base solver
+    // /////////////////////////////////////////////
+    if (typeid(*this->basesolver_) == typeid(LoopSolver<DiscFuncType>)) {
+
+        LoopSolver<DiscFuncType>* loopBaseSolver = dynamic_cast<LoopSolver<DiscFuncType>* > (this->basesolver_);
+     
+        typedef ProjectedBlockGSStep<OperatorType, DiscFuncType> SmootherType;
+        
+        dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->hasObstacle_ = &(*hasObstacle_)[0];
+        dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->obstacles_   = &(*obstacles_)[0];
+        
+#ifdef HAVE_IPOPT
+    } else if (typeid(*this->basesolver_) == typeid(QuadraticIPOptSolver<OperatorType,DiscFuncType>)) {
+
+        QuadraticIPOptSolver<OperatorType,DiscFuncType>* ipoptBaseSolver = dynamic_cast<QuadraticIPOptSolver<OperatorType,DiscFuncType>*> (this->basesolver_);
+     
+        ipoptBaseSolver->obstacles_   = (this->numLevels_ == 1) ? &(*obstacles_)[0]   : NULL;
+#endif
+    } else {
+        DUNE_THROW(SolverError, "You can't use " << typeid(*this->basesolver_).name()
+                   << " as a base solver for contact problems!");
+    }
+}
+
+template<class OperatorType, class DiscFuncType>
+void MonotoneMGStep<OperatorType, DiscFuncType>::nestedIteration()
+{
+    for (int i=this->numLevels_-1; i>0; i--) {
+
+        Dune::BitSetVector<dim> dummy(this->rhs_[i].size(), false);
+
+        // /////////////////////////////
+        //   Restrict obstacles
+        // //////////////////////////////
+        obstacleRestrictor_->restrict((*obstacles_)[i], 
+                                    (*obstacles_)[i-1],
+                                    (*hasObstacle_)[i], 
+                                    (*hasObstacle_)[i-1],
+                                    *this->mgTransfer_[i-1],
+                                    dummy);
+
+        // Restrict right hand side
+        this->mgTransfer_[i-1]->restrict(this->rhs_[i], this->rhs_[i-1]);
+    }
+
+    for (this->level_ = 0; this->level_<this->numLevels_-1; this->level_++) {
+
+        // If we start from an infeasible configuration, the restricted
+        // obstacles may be inconsistent.  We do an ad hoc correction here.
+        // The true obstacles on the maxlevel are not touched.
+        for (int i=0; i<(*obstacles_)[this->level_].size(); i++) {
+            BoxConstraint<field_type,dim>& cO = (*obstacles_)[this->level_][i];
+            for (int j=0; j<dim; j++)
+                if (cO.lower(j) > cO.upper(j))
+                    cO.lower(j) = cO.upper(j) = (cO.lower(j) + cO.upper(j)) / 2;
+        }
+        
+        std::cout << "Nested iteration on level " << this->level_ << std::endl;
+        iterate();
+        iterate();
+        //std::cout << this->x_[this->level_] << std::endl;
+        this->mgTransfer_[this->level_]->prolong(this->x_[this->level_], this->x_[this->level_+1]);
+        
+    }
+     
+}
+
+
+template <class OperatorType, class DiscFuncType>
+void MonotoneMGStep<OperatorType, DiscFuncType>::iterate()
+{
+    this->preprocessCalled = false;
+
+    int& level = this->level_;
+
+    // Define references just for ease of notation
+    std::vector<const OperatorType*>& mat = this->mat_;
+    std::vector<DiscFuncType>& x   = this->x_;
+    std::vector<DiscFuncType>& rhs = this->rhs_;
+    std::vector<std::vector<BoxConstraint<field_type,dim> > >& obstacles = *this->obstacles_;
+
+    Dune::BitSetVector<dim> critical(x[level].size(), false);
+
+    // Solve directly if we're looking at the coarse problem
+    if (level == 0) {
+
+        this->basesolver_->solve();
+
+        // Determine critical nodes  (only for debugging)
+        const double eps = 1e-12;
+        for (int i=0; i<obstacles[level].size(); i++) {
+            
+            for (int j=0; j<dim; j++) {
+                
+                if (obstacles[level][i].lower(j) >= x[level][i][j] - eps
+                    || obstacles[level][i].upper(j) <= x[level][i][j] + eps) {
+                    critical[i][j] = true;
+                }
+                
+            }
+        }
+
+    } else {
+    
+        // Presmoothing
+        this->presmoother_->setProblem(*(this->mat_[level]), x[level], rhs[level]);
+        dynamic_cast<ProjectedBlockGSStep<OperatorType, DiscFuncType>*>(this->presmoother_)->obstacles_ = &(obstacles[level]);
+        dynamic_cast<ProjectedBlockGSStep<OperatorType, DiscFuncType>*>(this->presmoother_)->hasObstacle_ = &((*hasObstacle_)[level]);
+        this->presmoother_->ignoreNodes_ = this->ignoreNodesHierarchy_[level];
+        
+        for (int i=0; i<this->nu1_; i++)
+            this->presmoother_->iterate();
+        
+        // First, a backup of the obstacles for postsmoothing
+        std::vector<BoxConstraint<field_type,dim> > obstacleBackup = obstacles[level];
+        
+        // Compute defect obstacles
+        for (int i=0; i<(*obstacles_)[level].size(); i++)
+            obstacles[level][i] -= x[level][i];
+        
+        // ///////////////////////
+        //    Truncation
+        // ///////////////////////
+        
+        // Determine critical nodes
+        const double eps = 1e-12;
+        for (int i=0; i<obstacles[level].size(); i++) {
+            
+            for (int j=0; j<dim; j++) {
+                
+                if (obstacles[level][i].lower(j) >= -eps
+                    || obstacles[level][i].upper(j) <= eps) {
+                    critical[i][j] = true;
+                }
+                
+            }
+        }
+        
+     
+        // Restrict stiffness matrix
+        dynamic_cast<TruncatedMGTransfer<DiscFuncType>*>(this->mgTransfer_[level-1])
+            ->galerkinRestrict(*(mat[level]), *(const_cast<OperatorType*>(mat[level-1])), critical);
+        
+        // Restrict obstacles
+        obstacleRestrictor_->restrict(obstacles[level], obstacles[level-1],
+                                    (*hasObstacle_)[level], (*hasObstacle_)[level-1],
+                                    *this->mgTransfer_[level-1],
+                                    critical);
+        
+        // //////////////////////////////////////////////////////////////////////
+        // Restriction:  fineResiduum = rhs[level] - mat[level] * x[level];
+        // //////////////////////////////////////////////////////////////////////
+        DiscFuncType fineResidual = rhs[level];
+        mat[level]->mmv(x[level], fineResidual);
+
+        // restrict residual
+        dynamic_cast<TruncatedMGTransfer<DiscFuncType>*>(this->mgTransfer_[level-1])->restrict(fineResidual, rhs[level-1], critical);
+        
+        // Choose all zeros as the initial correction
+        x[level-1] = 0;
+        
+        // ///////////////////////////////////////
+        // Recursively solve the coarser system
+        // ///////////////////////////////////////
+        level--;
+        for (int i=0; i<this->mu_; i++)
+            iterate();
+        level++;
+        
+        // ////////////////////////////////////////
+        // Prolong
+        // ////////////////////////////////////////
+        
+        // add correction to the presmoothed solution
+        DiscFuncType tmp;
+        dynamic_cast<TruncatedMGTransfer<DiscFuncType>*>(this->mgTransfer_[level-1])->prolong(x[level-1], tmp, critical);
+        x[level] += tmp;
+        
+        // restore the true (non-defect) obstacles
+        obstacles[level] = obstacleBackup;
+        
+#ifndef NDEBUG
+        // Debug: is the current iterate really admissible?
+        for (int i=0; i<obstacles[level].size(); i++)
+            for (int j=0; j<DiscFuncType::block_type::size; j++)
+                if (x[level][i][j] < obstacles[level][i].lower(j)
+                    || x[level][i][j] > obstacles[level][i].upper(j)) {
+                    
+                    std::cout << "Obstacle disregarded!\n";
+                    std::cout << x[level][i] << std::endl << obstacles[level][i] << std::endl;
+                    printf("level: %d   index: %d   komponent: %d\n", level, i, j);
+                    printf("is %s critical\n", (critical[i][j]) ? "" : "not");
+                }
+#endif
+        
+        // Postsmoothing
+        this->postsmoother_->setProblem(*(mat[level]), x[level], rhs[level]);
+        dynamic_cast<ProjectedBlockGSStep<OperatorType, DiscFuncType>*>(this->postsmoother_)->obstacles_ = &obstacles[level];
+        dynamic_cast<ProjectedBlockGSStep<OperatorType, DiscFuncType>*>(this->postsmoother_)->hasObstacle_ = &((*hasObstacle_)[level]);
+        this->postsmoother_->ignoreNodes_ = this->ignoreNodesHierarchy_[level];
+        
+        for (int i=0; i<this->nu2_; i++)
+            this->postsmoother_->iterate();
+        
+    }
+    
+    // ////////////////////////////////////////////////////////////////////
+    //   Track the number of critical nodes found during this iteration
+    // ////////////////////////////////////////////////////////////////////
+    if (level==this->numLevels_-1 && this->verbosity_==NumProc::FULL) {
+            
+        std::cout << critical.count() << " critical nodes found on level " << level;
+            
+        int changes = 0;
+        for (unsigned int i=0; i<oldCritical.size(); i++)
+            for (int j=0; j<dim; j++)
+                if (oldCritical[i][j] !=critical[i][j])
+                    changes++;
+        
+        std::cout << ", and " << changes << " changes." << std::endl;
+        oldCritical = critical;
+    }
+
+    // Debug: output energy
+    if (level==this->numLevels_-1 && this->verbosity_==NumProc::FULL)
+        std::cout << "Total energy: " 
+                  << std::setprecision(10) << computeEnergy(*mat[level], x[level], rhs[level]) << std::endl;
+
+}

dune-solvers/iterationsteps/mmgstep.hh

+#ifndef DUNE_MONOTONE_MULTIGRID_STEP_HH
+#define DUNE_MONOTONE_MULTIGRID_STEP_HH
+
+#include <dune/ag-common/multigridstep.hh>
+#include <dune/ag-common/transferoperators/multigridtransfer.hh>
+#include <dune/ag-common/projectedblockgsstep.hh>
+#include <dune/ag-common/boxconstraint.hh>
+#include <dune/ag-common/obstaclerestrictor.hh>
+
+/** \brief The general monotone multigrid solver
+*/
+template<class OperatorType, class DiscFuncType>
+class MonotoneMGStep : public MultigridStep<OperatorType, DiscFuncType>
+{
+    
+    static const int dim = DiscFuncType::block_type::dimension;
+    
+    typedef typename DiscFuncType::field_type field_type;
+    
+public:
+    
+    MonotoneMGStep() {} 
+    
+    MonotoneMGStep(int numLevels) {
+        this->numLevels_ = numLevels;
+        this->mat_.resize(numLevels);
+        this->x_.resize(numLevels);
+        this->rhs_.resize(numLevels);
+    }
+    
+    
+    MonotoneMGStep(const OperatorType& mat, 
+                   DiscFuncType& x, 
+                   DiscFuncType& rhs, int numLevels)
+        : MultigridStep<OperatorType, DiscFuncType>(mat, x, rhs, numLevels)
+    {
+        oldCritical.resize(x.size(), false);
+    }
+    
+    virtual ~MonotoneMGStep() {
+        // Delete the coarse grid matrices we created
+        for (int i=0; i<this->numLevels_-1; i++) {
+            delete(this->mat_[i]);
+            this->mat_[i] = NULL;
+        }
+    }
+    
+    virtual void setProblem(const OperatorType& mat, 
+                            DiscFuncType& x, 
+                            DiscFuncType& rhs,
+                            int numLevels) 
+    {
+        MultigridStep<OperatorType, DiscFuncType>::setProblem(mat,x,rhs, numLevels);
+        oldCritical.resize(x.size()*dim, false);
+    }
+    
+    virtual void iterate();
+    
+    virtual void preprocess();
+    
+    virtual void nestedIteration();
+    
+    ObstacleRestrictor<DiscFuncType>* obstacleRestrictor_;
+    
+    std::vector<Dune::BitSetVector<1> >* hasObstacle_;
+    
+    std::vector<Dune::BitSetVector<1> > recompute_;
+    
+    std::vector<std::vector<BoxConstraint<field_type,dim> > >* obstacles_;
+    
+    // Needed to track changes in the set of critical bits, and allows
+    // to check which dofs where critical after the last iteration
+    Dune::BitSetVector<dim> oldCritical;
+};
+
+#include "mmgstep.cc"
+
+#endif