diff --git a/dune-solvers/iterationsteps/mmgstep.cc b/dune-solvers/iterationsteps/mmgstep.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c8d38738a1546e29a9179bcc7663bd9eaf122e7a
--- /dev/null
+++ b/dune-solvers/iterationsteps/mmgstep.cc
@@ -0,0 +1,256 @@
+#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;
+
+}
diff --git a/dune-solvers/iterationsteps/mmgstep.hh b/dune-solvers/iterationsteps/mmgstep.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f704fe9a7ee7424cf1895227153977594eac5fe9
--- /dev/null
+++ b/dune-solvers/iterationsteps/mmgstep.hh
@@ -0,0 +1,78 @@
+#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