diff --git a/dune-solvers/iterationsteps/multigridstep.cc b/dune-solvers/iterationsteps/multigridstep.cc
new file mode 100644
index 0000000000000000000000000000000000000000..214ba2efaf8f5c19179b6e6d4cd4c9d559192029
--- /dev/null
+++ b/dune-solvers/iterationsteps/multigridstep.cc
@@ -0,0 +1,229 @@
+#include <typeinfo>
+
+#include <dune/ag-common/transferoperators/multigridtransfer.hh>
+#include <dune/ag-common/solvers/loopsolver.hh>
+#include "blockgsstep.hh"
+#include "dune/ag-common/genericvectortools.hh"
+
+#ifdef HAVE_IPOPT
+#include "quadraticipopt.hh"
+#endif
+
+template <class OperatorType, class VectorType, class BitVectorType>
+VectorType MultigridStep<OperatorType, VectorType, BitVectorType>::
+getSol()
+{
+ return this->x_[numLevels_-1];
+}
+
+template<class OperatorType, class VectorType, class BitVectorType>
+inline
+const OperatorType* MultigridStep<OperatorType, VectorType, BitVectorType>::
+getMatrix()
+{
+ return this->mat_[this->mat_.size()-1];
+}
+
+template<class OperatorType, class VectorType, class BitVectorType>
+inline
+void MultigridStep<OperatorType, VectorType, BitVectorType>::
+setMGType(int mu, int nu1, int nu2)
+{
+ mu_ = mu;
+ nu1_ = nu1;
+ nu2_ = nu2;
+}
+
+template<class OperatorType, class VectorType, class BitVectorType>
+void MultigridStep<OperatorType, VectorType, BitVectorType>::preprocess()
+{
+ if(preprocessCalled)
+ std::cout << "Preprocess has already been called without calling iterate afterwards!" << std::endl;
+ preprocessCalled = true;
+
+ int i;
+ this->level_ = this->numLevels_-1;
+
+ // //////////////////////////////////////////////////////////
+ // Check that transfer operators have been supplied
+ // //////////////////////////////////////////////////////////
+ if (this->mgTransfer_.size()!=this->numLevels_-1)
+ DUNE_THROW(SolverError, "You requested " << this->numLevels_ << " levels "
+ << "but you supplied " << this->mgTransfer_.size() << " transfer operators!");
+
+ for (i=0; i<this->mgTransfer_.size(); i++)
+ if (this->mgTransfer_[i] == NULL)
+ DUNE_THROW(SolverError, "You have not supplied a multigrid restriction operator "
+ "for level " << i);
+
+ // The matrices for the linear correction problems
+ for (i=this->numLevels_-2; i>=0; i--) {
+
+
+ // Delete coarse grid matrix if it exists already
+ if (this->mat_[i]) {
+ delete(this->mat_[i]);
+ this->mat_[i] = NULL;
+ }
+// int nDofs = this->mgTransfer_[i]->getMatrix().M();
+// this->mat_[i] = new OperatorType(nDofs, nDofs, OperatorType::row_wise);
+ this->mat_[i] = new OperatorType;
+
+ // Compute which entries are present in the (sparse) coarse grid stiffness
+ // matrices.
+ this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->mat_[i+1], *(const_cast<OperatorType*>(this->mat_[i])));
+
+ // Set solution vector sizes for the lower levels
+ GenericVector::resize(this->x_[i],*this->mat_[i]);
+ }
+
+ // /////////////////////////////////////////////////////
+ // Assemble the complete hierarchy of matrices
+ // /////////////////////////////////////////////////////
+ for (i=this->mgTransfer_.size()-1; i>=0; i--)
+ this->mgTransfer_[i]->galerkinRestrict(*this->mat_[i+1], *(const_cast<OperatorType*>(this->mat_[i])));
+
+ // /////////////////////////////////////////////////////
+ // Setup dirichlet bitfields
+ // /////////////////////////////////////////////////////
+ for (size_t i=0; i<ignoreNodesHierarchy_.size()-1; i++)
+ delete(ignoreNodesHierarchy_[i]);
+
+ if (this->ignoreNodes_!=0)
+ {
+ ignoreNodesHierarchy_[this->numLevels_-1] = const_cast<BitVectorType*>(this->ignoreNodes_);
+ for (i=this->numLevels_-2; i>=0; --i)
+ {
+ ignoreNodesHierarchy_[i] = new BitVectorType();
+ this->mgTransfer_[i]->restrictToFathers(*ignoreNodesHierarchy_[i+1], *ignoreNodesHierarchy_[i]);
+ }
+ } else
+ DUNE_THROW(SolverError, "We need a set of nodes to ignore");
+
+ for (i=this->mgTransfer_.size()-1; i>=0; i--)
+ this->mgTransfer_[i]->galerkinRestrict(*this->mat_[i+1], *(const_cast<OperatorType*>(this->mat_[i])));
+
+ // /////////////////////////////////////////////
+ // Set up base solver
+ // /////////////////////////////////////////////
+
+ if (this->basesolver_ == NULL)
+ DUNE_THROW(SolverError, "You have not provided a base solver!");
+
+ if (typeid(*this->basesolver_) == typeid(LoopSolver<VectorType>)) {
+
+ LoopSolver<VectorType>* loopBaseSolver = dynamic_cast<LoopSolver<VectorType>* > (this->basesolver_);
+
+ typedef LinearIterationStep<OperatorType, VectorType> SmootherType;
+ assert(dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_));
+
+ dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->setProblem(*(this->mat_[0]), this->x_[0], this->rhs_[0]);
+ dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->ignoreNodes_ = ignoreNodesHierarchy_[0];
+
+ }
+#ifdef HAVE_IPOPT
+ else if (typeid(*this->basesolver_) == typeid(QuadraticIPOptSolver<OperatorType,VectorType>)) {
+
+ QuadraticIPOptSolver<OperatorType,VectorType>* ipoptBaseSolver = dynamic_cast<QuadraticIPOptSolver<OperatorType,VectorType>*> (this->basesolver_);
+
+ ipoptBaseSolver->setProblem(*(this->mat_[0]), this->x_[0], this->rhs_[0]);
+ ipoptBaseSolver->ignoreNodes_ = ignoreNodesHierarchy_[0];
+
+ }
+#endif
+else {
+ DUNE_THROW(SolverError, "You can't use " << typeid(*this->basesolver_).name()
+ << " as a base solver for contact problems!");
+ }
+}
+
+
+
+template<class OperatorType, class VectorType, class BitVectorType>
+void MultigridStep<OperatorType, VectorType, BitVectorType>::nestedIteration()
+{
+
+ for (level_ = 0; level_<numLevels_; level_++) {
+
+ iterate();
+
+ mgTransfer_[level_]->prolong(x_[level_], x_[level_+1]);
+
+ }
+
+ iterate();
+}
+
+
+template<class OperatorType, class VectorType, class BitVectorType>
+void MultigridStep<OperatorType, VectorType, BitVectorType>::postprocess()
+{
+
+}
+
+
+template<class OperatorType, class VectorType, class BitVectorType>
+void MultigridStep<OperatorType, VectorType, BitVectorType>::iterate()
+{
+ preprocessCalled = false;
+
+ int& level = this->level_;
+
+ // Define references just for ease of notation
+ std::vector<const OperatorType*>& mat = this->mat_;
+ std::vector<VectorType>& x = this->x_;
+ std::vector<VectorType>& rhs = this->rhs_;
+
+ // Solve directly if we're looking at the coarse problem
+ if (level == 0) {
+ basesolver_->solve();
+ return;
+ }
+
+ // Presmoothing
+ presmoother_->setProblem(*(this->mat_[level]), x[level], rhs[level]);
+ presmoother_->ignoreNodes_ = ignoreNodesHierarchy_[level];
+
+ for (int i=0; i<this->nu1_; i++)
+ presmoother_->iterate();
+
+ // /////////////////////////
+ // Restriction
+
+ // compute residual
+ // fineResidual = rhs[level] - mat[level] * x[level];
+ VectorType fineResidual = rhs[level];
+ mat[level]->mmv(x[level], fineResidual);
+
+ // restrict residuum
+ this->mgTransfer_[level-1]->restrict(fineResidual, rhs[level-1]);
+
+
+ // Set Dirichlet values.
+ GenericVector::truncate(rhs[level-1], *ignoreNodesHierarchy_[level-1]);
+
+ // Choose all zeros as the initial correction
+ x[level-1] = 0;
+
+ // ///////////////////////////////////////
+ // Recursively solve the coarser system
+ level--;
+ iterate();
+ level++;
+
+ // ////////////////////////////////////////
+ // Prolong
+
+ // add correction to the presmoothed solution
+ VectorType tmp;
+ this->mgTransfer_[level-1]->prolong(x[level-1], tmp);
+ x[level] += tmp;
+
+ // Postsmoothing
+ postsmoother_->setProblem(*(mat[level]), x[level], rhs[level]);
+ postsmoother_->ignoreNodes_ = ignoreNodesHierarchy_[level];
+
+ for (int i=0; i<this->nu2_; i++)
+ postsmoother_->iterate();
+
+}
diff --git a/dune-solvers/iterationsteps/multigridstep.hh b/dune-solvers/iterationsteps/multigridstep.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8913f97f6ec5932b1efcda4a28bd87c191d8a17f
--- /dev/null
+++ b/dune-solvers/iterationsteps/multigridstep.hh
@@ -0,0 +1,214 @@
+#ifndef MULTIGRID_STEP_HH
+#define MULTIGRID_STEP_HH
+
+#include <vector>
+#include <dune/common/bitsetvector.hh>
+#include "iterativesolver.hh"
+#include "lineariterationstep.hh"
+
+#include <dune/ag-common/transferoperators/multigridtransfer.hh>
+
+ template<
+ class OperatorType,
+ class VectorType,
+ class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> >
+ class MultigridStep : public LinearIterationStep<OperatorType, VectorType, BitVectorType>
+ {
+
+ static const int blocksize = VectorType::block_type::dimension;
+
+ public:
+
+ MultigridStep() :
+ numLevels_(0),
+ presmoother_(0),
+ postsmoother_(0),
+ basesolver_(0),
+ ignoreNodesHierarchy_(0),
+ preprocessCalled(false)
+ {}
+
+ MultigridStep(const OperatorType& mat,
+ VectorType& x,
+ VectorType& rhs, int numLevels,
+ int mu, int nu1, int nu2,
+ LinearIterationStep<OperatorType, VectorType>* preSmoother,
+ LinearIterationStep<OperatorType, VectorType>* postSmoother,
+ Solver* baseSolver,
+ const BitVectorType* ignoreNodes) :
+ LinearIterationStep<OperatorType, VectorType>(mat, x, rhs),
+ mat_(numLevels),
+ ignoreNodesHierarchy_(numLevels),
+ x_(numLevels),
+ rhs_(numLevels),
+ preprocessCalled(false)
+ {
+ numLevels_ = numLevels;
+ level_ = numLevels-1;
+
+ mat_[level_] = &mat;
+ x_[level_] = x;
+ rhs_[level_] = rhs;
+
+ mu_ = mu;
+ nu1_ = nu1;
+ nu2_ = nu2;
+
+ presmoother_ = preSmoother;
+ postsmoother_ = postSmoother;
+ basesolver_ = baseSolver;
+
+ this->ignoreNodes_ = ignoreNodes;
+ }
+
+ MultigridStep(const OperatorType& mat,
+ VectorType& x,
+ VectorType& rhs, int numLevels) :
+ LinearIterationStep<OperatorType, VectorType>(mat, x, rhs),
+ presmoother_(0),
+ postsmoother_(0),
+ basesolver_(0),
+ mat_(numLevels),
+ ignoreNodesHierarchy_(numLevels),
+ x_(numLevels),
+ rhs_(numLevels),
+ preprocessCalled(false)
+ {
+ numLevels_ = numLevels;
+ level_ = numLevels-1;
+
+ mat_[level_] = &mat;
+ x_[level_] = x;
+ rhs_[level_] = rhs;
+ }
+
+
+ virtual ~MultigridStep() {
+ // Delete all matrices but the top one, which has been supplied from the outside
+ for (size_t i=0; i<mat_.size()-1; i++)
+ {
+ delete(mat_[i]);
+ delete(ignoreNodesHierarchy_[i]);
+ }
+ }
+
+ virtual void setProblem(const OperatorType& mat,
+ VectorType& x,
+ VectorType& rhs,
+ int numLevels)
+ {
+ setNumberOfLevels(numLevels);
+ setProblem(mat, x, rhs);
+ }
+
+ virtual void setNumberOfLevels(int numLevels)
+ {
+ numLevels_ = numLevels;
+
+ for (int i=0; i<int(mat_.size())-1; i++)
+ if (mat_[i])
+ delete(mat_[i]);
+
+ for (int i=0; i<int(ignoreNodesHierarchy_.size())-1; i++)
+ if (ignoreNodesHierarchy_[i])
+ delete(ignoreNodesHierarchy_[i]);
+
+ mat_.resize(numLevels);
+ ignoreNodesHierarchy_.resize(numLevels);
+
+ for (int i=0; i<int(mat_.size())-1; i++)
+ {
+ mat_[i] = NULL;
+ ignoreNodesHierarchy_[i] = NULL;
+ }
+
+ x_.resize(numLevels);
+ rhs_.resize(numLevels);
+ }
+
+ virtual void setProblem(const OperatorType& mat,
+ VectorType& x,
+ VectorType& rhs)
+ {
+ level_ = numLevels_-1;
+
+ mat_[level_] = &mat;
+ x_[level_] = x;
+ rhs_[level_] = rhs;
+ }
+
+ template <class DerivedTransferHierarchy>
+ void setTransferOperators(const DerivedTransferHierarchy& transfer)
+ {
+ mgTransfer_.resize(transfer.size());
+ for(int j=0; j<transfer.size(); ++j)
+ mgTransfer_[j] = transfer[j];
+ }
+
+ virtual void iterate();
+
+ virtual void nestedIteration();
+
+ virtual void preprocess();
+
+ virtual void postprocess();
+
+ virtual VectorType getSol();
+
+ virtual const OperatorType* getMatrix();
+
+ virtual int level() const {return level_;}
+
+ /** \brief Sets the number of pre- and postsmoothing steps
+ and of coarse corrections.
+ \param mu Number of coarse corrections
+ \param nu1 Number of presmoothing steps
+ \param nu2 Number of postsmoothing steps
+ */
+ virtual void setMGType(int mu, int nu1, int nu2);
+
+ LinearIterationStep<OperatorType, VectorType>* presmoother_;
+
+ LinearIterationStep<OperatorType, VectorType>* postsmoother_;
+
+ Solver* basesolver_;
+
+ protected:
+ //! Number of presmoothing steps
+ int nu1_;
+
+ //! Number of postsmoothing steps
+ int nu2_;
+
+ //! Number of coarse corrections
+ int mu_;
+ public:
+ //! Variable used to store the current level
+ int level_;
+
+ //! The total number of levels of the underlying grid
+ int numLevels_;
+
+ public:
+ //! The linear operators on each level
+ std::vector<const OperatorType*> mat_;
+
+ protected:
+ //! Flags specifying the dirichlet nodes on each level
+ std::vector<BitVectorType*> ignoreNodesHierarchy_;
+
+ public:
+ std::vector<VectorType> x_;
+
+ std::vector<VectorType> rhs_;
+
+ //protected:
+ std::vector<MultigridTransfer<VectorType, BitVectorType, OperatorType>* > mgTransfer_;
+
+ protected:
+ bool preprocessCalled;
+ };
+
+#include "multigridstep.cc"
+
+#endif