diff --git a/dune/solvers/test/Makefile.am b/dune/solvers/test/Makefile.am
index 2e23ad000fb1edc8d382c400d10d92156690f02f..3b8b4eaf102fdb6a56e073b87875f545e1477788 100644
--- a/dune/solvers/test/Makefile.am
+++ b/dune/solvers/test/Makefile.am
@@ -2,12 +2,13 @@
# list of tests to run
TESTS = lowrankoperatortest \
nulloperatortest \
- sumoperatortest
+ sumoperatortest \
+ obstacletnnmgtest
# programs just to build when "make check" is used
check_PROGRAMS = $(TESTS)
-noinst_HEADERS =
+noinst_HEADERS = common.hh
# collect some common flags
COMMON_CPPFLAGS = $(AM_CPPFLAGS) $(DUNEMPICPPFLAGS) -I$(top_srcdir)
@@ -35,4 +36,9 @@ sumoperatortest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
sumoperatortest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
sumoperatortest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS)
+obstacletnnmgtest_SOURCES = obstacletnnmgtest.cc
+obstacletnnmgtest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
+obstacletnnmgtest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)
+obstacletnnmgtest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS)
+
include $(top_srcdir)/am/global-rules
diff --git a/dune/solvers/test/common.hh b/dune/solvers/test/common.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7b68d5fbc1e9f341d9ac746e0572a01e1a4dddc5
--- /dev/null
+++ b/dune/solvers/test/common.hh
@@ -0,0 +1,345 @@
+#ifndef DUNE_SOLVERS_TESTS_COMMON_HH
+#define DUNE_SOLVERS_TESTS_COMMON_HH
+
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/istl/matrixindexset.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+
+template<class GridView, class Matrix>
+void constructPQ1Pattern(const GridView& gridView, Matrix& matrix)
+{
+ static const int dim = GridView::Grid::dimension;
+
+ typedef typename GridView::IndexSet IndexSet;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
+ typedef typename FiniteElementCache::FiniteElementType FiniteElement;
+
+ const IndexSet& indexSet = gridView.indexSet();
+ FiniteElementCache cache;
+
+ int size = indexSet.size(dim);
+
+ Dune::MatrixIndexSet indices(size, size);
+
+ ElementIterator it = gridView.template begin<0>();
+ ElementIterator end = gridView.template end<0>();
+ for (; it != end; ++it)
+ {
+ Element& e = *it;
+ const FiniteElement& fe = cache.get(e.type());
+
+ int localSize = fe.localBasis().size();
+ for (int i = 0; i < localSize; ++i)
+ {
+ int iGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(i).subEntity(), dim);
+ for (int j = 0; j < localSize; ++j)
+ {
+ int jGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(j).subEntity(), dim);
+ indices.add(iGlobal, jGlobal);
+ indices.add(jGlobal, iGlobal);
+ }
+ }
+ }
+ indices.exportIdx(matrix);
+}
+
+
+template<class GridView, class Matrix>
+void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix)
+{
+ static const int dim = GridView::Grid::dimension;
+ static const int dimworld = GridView::Grid::dimensionworld;
+
+ typedef typename GridView::IndexSet IndexSet;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
+ typedef typename FiniteElementCache::FiniteElementType FiniteElement;
+
+ typedef typename Dune::QuadratureRule<double,dim> QuadratureRule;
+ typedef typename Dune::template FieldVector<double,dim> LocalCoordinate;
+ typedef typename Dune::template FieldVector<double,dimworld> GlobalCoordinate;
+ typedef typename Dune::template FieldMatrix<double,dimworld,dim> InverseJacobianTransposed;
+ typedef typename FiniteElement::Traits::LocalBasisType::Traits::JacobianType JacobianType;
+
+ const IndexSet& indexSet = gridView.indexSet();
+ FiniteElementCache cache;
+
+ ElementIterator it = gridView.template begin<0>();
+ ElementIterator end = gridView.template end<0>();
+ for (; it != end; ++it)
+ {
+ Element& e = *it;
+ const FiniteElement& fe = cache.get(e.type());
+
+ int localSize = fe.localBasis().size();
+
+ // get quadrature rule of appropiate order (P1/Q1)
+ int order = 0;
+ if (e.type().isSimplex())
+ order = 2*(1-1);
+ else
+ order = 2*(dim-1);
+ const QuadratureRule& quad = Dune::template QuadratureRules<double, dim>::rule(e.type(), order);
+
+ // store gradients of shape functions and base functions
+ std::vector<JacobianType> referenceGradients(localSize);
+ std::vector<GlobalCoordinate> gradients(localSize);
+
+ for (size_t pt=0; pt < quad.size(); ++pt)
+ {
+ // get quadrature point
+ const LocalCoordinate& quadPos = quad[pt].position();
+
+ // get transposed inverse of Jacobian of transformation
+ const InverseJacobianTransposed& invJacobian = e.geometry().jacobianInverseTransposed(quadPos);
+
+ // get integration factor
+ const double integrationElement = e.geometry().integrationElement(quadPos);
+
+ // evaluate gradients of shape functions
+ fe.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+ // transform gradients
+ for (size_t i=0; i<gradients.size(); ++i)
+ invJacobian.mv(referenceGradients[i][0], gradients[i]);
+
+ // compute matrix entries
+ double z = quad[pt].weight() * integrationElement;
+ for (int i = 0; i < localSize; ++i)
+ {
+ int iGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(i).subEntity(), dim);
+ for (int j = i+1; j < localSize; ++j)
+ {
+ int jGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(j).subEntity(), dim);
+
+ double zij = (gradients[i] * gradients[j]) * z;
+ matrix[iGlobal][jGlobal] += zij;
+ matrix[jGlobal][iGlobal] += zij;
+ }
+ double zii = (gradients[i] * gradients[i]) * z;
+ matrix[iGlobal][iGlobal] += zii;
+ }
+ }
+ }
+}
+
+
+template<class GridView, class Matrix>
+void assemblePQ1Mass(const GridView& gridView, Matrix& matrix)
+{
+ static const int dim = GridView::Grid::dimension;
+
+ typedef typename GridView::IndexSet IndexSet;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
+ typedef typename FiniteElementCache::FiniteElementType FiniteElement;
+
+ typedef typename Dune::QuadratureRule<double,dim> QuadratureRule;
+ typedef typename Dune::template FieldVector<double,dim> LocalCoordinate;
+ typedef typename FiniteElement::Traits::LocalBasisType::Traits::RangeType RangeType;
+
+ const IndexSet& indexSet = gridView.indexSet();
+ FiniteElementCache cache;
+
+ ElementIterator it = gridView.template begin<0>();
+ ElementIterator end = gridView.template end<0>();
+ for (; it != end; ++it)
+ {
+ Element& e = *it;
+ const FiniteElement& fe = cache.get(e.type());
+
+ int localSize = fe.localBasis().size();
+
+ // get quadrature rule of appropiate order (P1/Q1)
+ int order = 0;
+ if (e.type().isSimplex())
+ order = 2*1;
+ else
+ order = 2*dim;
+ const QuadratureRule& quad = Dune::template QuadratureRules<double, dim>::rule(e.type(), order);
+
+ // store values of shape functions
+ std::vector<RangeType> values(localSize);
+
+ for (size_t pt=0; pt < quad.size(); ++pt)
+ {
+ // get quadrature point
+ const LocalCoordinate& quadPos = quad[pt].position();
+
+ // get integration factor
+ const double integrationElement = e.geometry().integrationElement(quadPos);
+
+ // evaluate basis functions
+ fe.localBasis().evaluateFunction(quadPos, values);
+
+ // compute matrix entries
+ double z = quad[pt].weight() * integrationElement;
+ for (int i = 0; i < localSize; ++i)
+ {
+ int iGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(i).subEntity(), dim);
+ double zi = values[i]*z;
+ for (int j = i+1; j < localSize; ++j)
+ {
+ int jGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(j).subEntity(), dim);
+
+ double zij = values[j] * zi;
+ matrix[iGlobal][jGlobal] += zij;
+ matrix[jGlobal][iGlobal] += zij;
+ }
+ double zii = values[i] * zi;
+ matrix[iGlobal][iGlobal] += zii;
+ }
+ }
+ }
+}
+
+
+template<class GridView, class Vector, class Function>
+void assemblePQ1RHS(const GridView& gridView, Vector& r, const Function& f)
+{
+ static const int dim = GridView::Grid::dimension;
+
+ typedef typename GridView::IndexSet IndexSet;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
+ typedef typename FiniteElementCache::FiniteElementType FiniteElement;
+
+ typedef typename Dune::QuadratureRule<double,dim> QuadratureRule;
+ typedef typename Dune::template FieldVector<double,dim> LocalCoordinate;
+ typedef typename FiniteElement::Traits::LocalBasisType::Traits::RangeType RangeType;
+ typedef typename Function::RangeType FunctionRangeType;
+
+ const IndexSet& indexSet = gridView.indexSet();
+ FiniteElementCache cache;
+
+ int size = indexSet.size(dim);
+
+ ElementIterator it = gridView.template begin<0>();
+ ElementIterator end = gridView.template end<0>();
+ for (; it != end; ++it)
+ {
+ Element& e = *it;
+ const FiniteElement& fe = cache.get(e.type());
+
+ int localSize = fe.localBasis().size();
+
+ // get quadrature rule of appropiate order (P1/Q1)
+ int order = 0;
+ if (e.type().isSimplex())
+ order = 2*1;
+ else
+ order = 2*dim;
+ const QuadratureRule& quad = Dune::template QuadratureRules<double, dim>::rule(e.type(), order);
+
+ // store values of shape functions
+ std::vector<RangeType> values(localSize);
+
+ for (size_t pt=0; pt < quad.size(); ++pt)
+ {
+ // get quadrature point
+ const LocalCoordinate& quadPos = quad[pt].position();
+
+ // get integration factor
+ const double integrationElement = e.geometry().integrationElement(quadPos);
+
+ // evaluate basis functions
+ fe.localBasis().evaluateFunction(quadPos, values);
+
+ // evaluate function
+ FunctionRangeType fAtPos;
+ f.evaluate(e.geometry().global(quadPos), fAtPos);
+
+ // add vector entries
+ double z = quad[pt].weight() * integrationElement;
+ for (int i = 0; i < localSize; ++i)
+ {
+ int iGlobal = indexSet.subIndex(e, fe.localCoefficients().localKey(i).subEntity(), dim);
+
+ r[iGlobal].axpy(values[i]*z, fAtPos);
+ }
+ }
+ }
+}
+
+
+template<class Intersection>
+bool intersectionContainsVertex(const Intersection& i, int vertexInInsideElement)
+{
+ static const int dim = Intersection::dimension;
+
+ int faceIdx = i.indexInInside();
+
+ const Dune::GenericReferenceElement<double,dim>& refElement
+ = Dune::GenericReferenceElements<double, dim>::general(i.inside()->type());
+
+ for (int j = 0; j<refElement.size(faceIdx, 1, dim); ++j)
+ {
+ int intersectionVertexInElement = refElement.subEntity(faceIdx, 1, j, dim);
+
+ if (intersectionVertexInElement == vertexInInsideElement)
+ return true;
+ }
+ return false;
+}
+
+
+template<class GridView, class BitVector>
+void markBoundaryDOFs(const GridView& gridView, BitVector& isBoundary)
+{
+ static const int dim = GridView::Grid::dimension;
+
+ typedef typename GridView::IndexSet IndexSet;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename GridView::Intersection Intersection;
+ typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
+ typedef typename FiniteElementCache::FiniteElementType FiniteElement;
+
+ const IndexSet& indexSet = gridView.indexSet();
+ FiniteElementCache cache;
+
+ ElementIterator it = gridView.template begin<0>();
+ ElementIterator end = gridView.template end<0>();
+ for (; it != end; ++it)
+ {
+ Element& e = *it;
+ const FiniteElement& fe = cache.get(e.type());
+ int localSize = fe.localBasis().size();
+
+ IntersectionIterator nIt = gridView.ibegin(e);
+ IntersectionIterator nEnd = gridView.iend(e);
+
+ for (; nIt!=nEnd; ++nIt)
+ {
+ const Intersection& i = *nIt;
+
+ if (i.boundary())
+ {
+ for (int j = 0; j < localSize; ++j)
+ {
+ unsigned int vertexInInsideElement = fe.localCoefficients().localKey(j).subEntity();
+ int iGlobal = indexSet.subIndex(e, vertexInInsideElement, dim);
+ if (intersectionContainsVertex(i, vertexInInsideElement))
+ isBoundary[iGlobal] = true;
+ }
+ }
+ }
+ }
+}
+
+
+
+#endif
diff --git a/dune/solvers/test/obstacletnnmgtest.cc b/dune/solvers/test/obstacletnnmgtest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b3126ba02eda942394a9d493c1053b06c52536ea
--- /dev/null
+++ b/dune/solvers/test/obstacletnnmgtest.cc
@@ -0,0 +1,206 @@
+#include <config.h>
+
+// 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/iterationsteps/obstacletnnmgstep.hh>
+#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+
+#include "common.hh"
+
+#define DIMENSION 2
+
+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>
+void solveObstacleProblemByTNNMG(const GridType& grid, const MatrixType& mat, VectorType& x, const VectorType& rhs, int maxIterations=100, double tolerance=1.0e-10)
+{
+ // bit vector marking dofs to ignore
+ typedef typename ObstacleTNNMGStep<MatrixType, VectorType>::BitVector BitVector;
+ BitVector ignore(rhs.size());
+ ignore.unsetAll();
+
+ solveObstacleProblemByTNNMG(grid, mat, x, rhs, ignore, maxIterations, tolerance);
+}
+
+
+template<class GridType, class MatrixType, class VectorType, class BitVector>
+void solveObstacleProblemByTNNMG(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 LoopSolver<Vector> Solver;
+ typedef ObstacleTNNMGStep<Matrix, Vector> TNNMGStep;
+ typedef typename TNNMGStep::Transfer Transfer;
+ typedef typename TNNMGStep::BoxConstraintVector BoxConstraintVector;
+ typedef CompressedMultigridTransfer<Vector, BitVector, Matrix> TransferImplementation;
+
+ const int blockSize = VectorType::block_type::dimension;
+
+ // we need a vector of pointers to the transfer operator base class
+ std::vector<Transfer*> transfer(grid.maxLevel());
+ for (int i = 0; i < transfer.size(); ++i)
+ {
+ // create transfer operator from level i to i+1
+ TransferImplementation* t = new TransferImplementation;
+ t->setup(grid, i, i+1);
+ transfer[i] = t;
+ }
+
+ // create double obstacle constraints
+ 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;
+ }
+ }
+
+ // we want to control errors with respect to the energy norm induced by the matrix
+ Norm norm(mat);
+
+ // create iteration step
+ TNNMGStep tnnmgStep(mat, x, rhs, ignore, boxConstraints, transfer);
+ tnnmgStep.preprocess();
+ tnnmgStep.nestedIteration();
+
+ // cretae loop solver
+ Solver solver(&tnnmgStep, maxIterations, tolerance, &norm, Solver::FULL);
+
+ // solve problem
+ solver.solve();
+}
+
+
+
+
+
+template <class GridType>
+bool checkWithGrid(const GridType& grid)
+{
+ 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);
+
+
+ solveObstacleProblemByTNNMG(grid, A, u, rhs, ignore);
+
+ typename Dune::VTKWriter<GridView> vtkWriter(gridView);
+ vtkWriter.addVertexData(u, "solution");
+ vtkWriter.write("obstacletnnmgtest_solution");
+
+
+ return passed;
+}
+
+
+template <int dim>
+bool checkWithYaspGrid(int refine)
+{
+ 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);
+
+
+ return passed;
+}
+
+
+
+
+int main(int argc, char** argv) try
+{
+ static const int dim = DIMENSION;
+
+ bool passed(true);
+
+
+ passed = passed and checkWithYaspGrid<1>(15);
+ passed = passed and checkWithYaspGrid<2>(10);
+ passed = passed and checkWithYaspGrid<3>(5);
+
+ return passed ? 0 : 1;
+}
+catch (Dune::Exception e) {
+
+ std::cout << e << std::endl;
+
+}
+