diff --git a/dune/solvers/test/Makefile.am b/dune/solvers/test/Makefile.am
index 81abd3f17c35a2df474f3451d86d553561ad37cf..a4cdd073bb13ceef6c8958394d0a3bb8e596663b 100644
--- a/dune/solvers/test/Makefile.am
+++ b/dune/solvers/test/Makefile.am
@@ -2,6 +2,7 @@
# list of tests to run
TESTS = genericvectortoolstest \
lowrankoperatortest \
+ mmgtest \
multigridtest \
nulloperatortest \
obstacletnnmgtest \
@@ -34,6 +35,11 @@ lowrankoperatortest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
lowrankoperatortest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
lowrankoperatortest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS)
+mmgtest_SOURCES = mmgtest.cc
+mmgtest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS) $(IPOPT_CPPFLAGS)
+mmgtest_LDADD = $(COMMON_LDADD) $(GRID_LDADD) $(IPOPT_LIBS)
+mmgtest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS) $(IPOPT_LIBS)
+
multigridtest_SOURCES = multigridtest.cc
multigridtest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
multigridtest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
diff --git a/dune/solvers/test/mmgtest.cc b/dune/solvers/test/mmgtest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ce8d5e4294e4847746ae8df45ea100b6e4d2f3f4
--- /dev/null
+++ b/dune/solvers/test/mmgtest.cc
@@ -0,0 +1,217 @@
+#include <config.h>
+
+#include <iostream>
+#include <sstream>
+
+// dune-common includes
+#include <dune/common/bitsetvector.hh>
+
+// dune-istl includes
+#include <dune/istl/bcrsmatrix.hh>
+
+// dune-grid includes
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+// dune-solver includes
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/iterationsteps/mmgstep.hh>
+#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
+
+
+#include "common.hh"
+
+template <class DomainType, class RangeType>
+class ConstantFunction :
+ public Dune::VirtualFunction<DomainType, RangeType>
+{
+ public:
+ ConstantFunction(double c):
+ c_(c)
+ {}
+
+ void evaluate(const DomainType& x, RangeType& y) const
+ {
+ y = c_;
+ }
+
+ private:
+ double c_;
+};
+
+
+template<class GridType, class MatrixType, class VectorType, class BitVector>
+void solveObstacleProblemByMMGSolver(const GridType& grid, const MatrixType& mat, VectorType& x, const VectorType& rhs, const BitVector& ignore, int maxIterations=100, double tolerance=1.0e-10)
+{
+ typedef VectorType Vector;
+ typedef MatrixType Matrix;
+ typedef EnergyNorm<Matrix, Vector> Norm;
+ typedef QuadraticIPOptSolver<MatrixType,VectorType> Solver;
+
+ // create double obstacle constraints
+#if 0
+ const int blockSize = VectorType::block_type::dimension;
+
+ BoxConstraintVector boxConstraints(rhs.size());
+ for (int i = 0; i < boxConstraints.size(); ++i)
+ {
+ for (int j = 0; j < blockSize; ++j)
+ {
+ boxConstraints[i].lower(j) = -1;
+ boxConstraints[i].upper(j) = +1;
+ }
+ }
+#endif
+
+ /////////////////////////////////////////////////////////////
+ // Create a monotone multigrid solver
+ /////////////////////////////////////////////////////////////
+
+ // First create a Gauss-seidel base solver
+ ProjectedBlockGSStep<MatrixType, VectorType>* baseSolverStep = new ProjectedBlockGSStep<MatrixType, VectorType>;
+
+ EnergyNorm<MatrixType, VectorType>* baseEnergyNorm = new EnergyNorm<MatrixType, VectorType>(*baseSolverStep);
+
+ ::LoopSolver<VectorType>* baseSolver = new ::LoopSolver<VectorType>(baseSolverStep,
+ 100, // base iterations
+ 1e-8, // tolerance
+ baseEnergyNorm,
+ Solver::QUIET);
+
+ // Make pre and postsmoothers
+ ProjectedBlockGSStep<MatrixType, VectorType>* presmoother = new ProjectedBlockGSStep<MatrixType, VectorType>;
+ ProjectedBlockGSStep<MatrixType, VectorType>* postsmoother = new ProjectedBlockGSStep<MatrixType, VectorType>;
+
+ MonotoneMGStep<MatrixType, VectorType>* mmgStep = new MonotoneMGStep<MatrixType, VectorType>();
+
+ mmgStep->setMGType(3, 3, 1);
+ mmgStep->ignoreNodes_ = &ignore;
+ mmgStep->basesolver_ = baseSolver;
+ mmgStep->setSmoother(presmoother, postsmoother);
+ mmgStep->obstacleRestrictor_= new MandelObstacleRestrictor<VectorType>();
+ mmgStep->hasObstacle_ = NULL;
+ mmgStep->verbosity_ = Solver::QUIET;
+
+ EnergyNorm<MatrixType,VectorType>* h1SemiNorm = new EnergyNorm<MatrixType,VectorType>(mat);
+
+ LoopSolver<VectorType> solver (mmgStep,
+ 100, // iterations
+ 1e-6, // tolerance
+ h1SemiNorm,
+ Solver::REDUCED);
+ // solve problem
+ solver.solve();
+}
+
+
+
+
+
+template <class GridType>
+bool checkWithGrid(const GridType& grid, const std::string fileName="")
+{
+ bool passed = true;
+
+ static const int dim = GridType::dimension;
+
+ typedef typename Dune::FieldMatrix<double, 1, 1> MatrixBlock;
+ typedef typename Dune::FieldVector<double, 1> VectorBlock;
+ typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
+ typedef typename Dune::BlockVector<VectorBlock> Vector;
+ typedef typename Dune::BitSetVector<1> BitVector;
+
+ typedef typename GridType::LeafGridView GridView;
+ typedef typename Dune::FieldVector<double, dim> DomainType;
+ typedef typename Dune::FieldVector<double, 1> RangeType;
+
+
+ const GridView gridView = grid.leafView();
+
+ Matrix A;
+ constructPQ1Pattern(gridView, A);
+ A=0.0;
+ assemblePQ1Stiffness(gridView, A);
+
+ Vector rhs(A.N());
+ rhs = 0;
+ ConstantFunction<DomainType, RangeType> f(50);
+ assemblePQ1RHS(gridView, rhs, f);
+
+ Vector u(A.N());
+ u = 0;
+
+ BitVector ignore(A.N());
+ ignore.unsetAll();
+ markBoundaryDOFs(gridView, ignore);
+
+
+ solveObstacleProblemByMMGSolver(grid, A, u, rhs, ignore);
+
+ if (fileName!="")
+ {
+ typename Dune::VTKWriter<GridView> vtkWriter(gridView);
+ vtkWriter.addVertexData(u, "solution");
+ vtkWriter.write(fileName);
+ }
+
+ return passed;
+}
+
+
+template <int dim>
+bool checkWithYaspGrid(int refine, const std::string fileName="")
+{
+ bool passed = true;
+
+ typedef Dune::YaspGrid<dim> GridType;
+
+ typename Dune::FieldVector<typename GridType::ctype,dim> L(1.0);
+ typename Dune::FieldVector<int,dim> s(1);
+ typename Dune::FieldVector<bool,dim> periodic(false);
+
+ GridType grid(L, s, periodic, 0);
+
+ for (int i = 0; i < refine; ++i)
+ grid.globalRefine(1);
+
+ std::cout << "Testing with YaspGrid<" << dim << ">" << std::endl;
+ std::cout << "Number of levels: " << (grid.maxLevel()+1) << std::endl;
+ std::cout << "Number of leaf nodes: " << grid.leafView().size(dim) << std::endl;
+
+ passed = passed and checkWithGrid(grid, fileName);
+
+
+ return passed;
+}
+
+
+
+
+int main(int argc, char** argv) try
+{
+ bool passed(true);
+
+
+ int refine1d = 16;
+ int refine2d = 8;
+ int refine3d = 5;
+
+ if (argc>1)
+ std::istringstream(argv[1]) >> refine1d;
+ if (argc>2)
+ std::istringstream(argv[2]) >> refine2d;
+ if (argc>3)
+ std::istringstream(argv[3]) >> refine3d;
+
+ passed = passed and checkWithYaspGrid<1>(refine1d, "obstacletnnmgtest_yasp_1d_solution");
+ passed = passed and checkWithYaspGrid<2>(refine2d, "obstacletnnmgtest_yasp_2d_solution");
+ passed = passed and checkWithYaspGrid<3>(refine3d, "obstacletnnmgtest_yasp_3d_solution");
+
+ return passed ? 0 : 1;
+}
+catch (Dune::Exception e) {
+
+ std::cout << e << std::endl;
+
+}
+