dune/tectonic/tectonic.hh

-#ifndef DUNE_TECTONIC_TECTONIC_HH
-#define DUNE_TECTONIC_TECTONIC_HH
-
-// add your classes here
-#endif

dune/tectonic/tests/CMakeLists.txt

+add_subdirectory("contacttest")
+add_subdirectory("hdf5test")
+#add_subdirectory("tnnmgtest")
+
+dune_add_test(SOURCES globalfrictioncontainertest.cc)
+dune_add_test(SOURCES gridgluefrictiontest.cc)
+dune_add_test(SOURCES nodalweightstest.cc)
+dune_add_test(SOURCES supportpatchfactorytest.cc)
+dune_add_test(SOURCES solverfactorytest.cc)

dune/tectonic/tests/common.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_SRC_TESTS_COMMON_HH
+#define DUNE_TECTONIC_SRC_TESTS_COMMON_HH
+
+#include <dune/grid/common/gridfactory.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#include <dune/fufem/referenceelementhelper.hh>
+
+// utility functions
+bool isClose(double a, double b) {
+  return std::abs(a - b) < 1e-14;
+}
+
+template <class Vector>
+bool xyCollinear(Vector const &a, Vector const &b, Vector const &c) {
+  return isClose((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
+}
+
+template <class Vector>
+bool xyBoxed(Vector const &v1, Vector const &v2, Vector const &x) {
+  auto const minmax0 = std::minmax(v1[0], v2[0]);
+  auto const minmax1 = std::minmax(v1[1], v2[1]);
+
+  if (minmax0.first - 1e-14 > x[0] or
+      x[0] > minmax0.second + 1e-14)
+    return false;
+  if (minmax1.first - 1e-14 > x[1] or
+      x[1] > minmax1.second + 1e-14)
+    return false;
+
+  return true;
+}
+
+template <class Vector>
+bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x) {
+  return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
+}
+
+// build cube grid given by lowerLeft and upperRight point in global coordinates
+template <class GlobalCoords, class GridType>
+void buildGrid(const GlobalCoords& lowerLeft, const GlobalCoords& upperRight, const int n, std::shared_ptr<GridType>& grid, bool simplexGrid = true) {
+    std::array<unsigned int, GridType::dimension> elements;
+    std::fill(elements.begin(), elements.end(), n);
+
+    Dune::GridFactory<GridType> factory;
+    if (!simplexGrid) {
+        Dune::StructuredGridFactory<GridType>::createCubeGrid(factory, lowerLeft, upperRight, elements);
+        grid = std::move(factory.createGrid());
+    } else {
+        Dune::StructuredGridFactory<GridType>::createSimplexGrid(factory, lowerLeft, upperRight, elements);
+        grid = std::move(factory.createGrid());
+    }
+}
+
+template <class Intersection>
+bool containsInsideSubentity(const Intersection& nIt, int subEntity, int codim) {
+    return ReferenceElementHelper<double, Intersection::GridView::dim>::subEntityContainsSubEntity(nIt.inside().type(), nIt.indexInInside(), 1, subEntity, codim);
+}
+
+#endif

dune/tectonic/tests/contacttest/CMakeLists.txt

+add_custom_target(tectonic_tests_contacttest SOURCES
+  staticcontacttest.cfg
+  staticcontacttest-2D.cfg
+  staticcontacttest-3D.cfg
+)
+
+set(CONTACTTEST_SOURCE_FILES
+  ../../assemblers.cc
+  ../../data-structures/body/body.cc
+  ../../data-structures/network/levelcontactnetwork.cc
+  ../../data-structures/network/contactnetwork.cc
+  ../../data-structures/enumparser.cc
+  ../../factories/twoblocksfactory.cc
+  ../../factories/threeblocksfactory.cc
+  ../../factories/stackedblocksfactory.cc
+  #../../io/vtk.cc
+  ../../problem-data/grid/cuboidgeometry.cc
+  ../../problem-data/grid/mygrids.cc
+  ../../problem-data/grid/simplexmanager.cc
+  ../../spatial-solving/solverfactory.cc
+  staticcontacttest.cc
+)
+
+foreach(_dim 2 3)
+  set(_contacttest_target staticcontacttest-${_dim}D)
+
+  dune_add_test(NAME ${_contacttest_target} SOURCES ${CONTACTTEST_SOURCE_FILES})
+
+  add_dune_ug_flags(${_contacttest_target})
+  add_dune_hdf5_flags(${_contacttest_target})
+
+  set_property(TARGET ${_contacttest_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+endforeach()

dune/tectonic/tests/contacttest/staticcontacttest-2D.cfg

+# -*- mode:conf -*-
+[boundary.friction]
+smallestDiameter = 0.3  # 2e-3 [m]
+
+[u0.solver]
+tolerance         = 1e-8
+
+[solver.tnnmg.preconditioner.basesolver]
+tolerance          = 1e-10
+
+[solver.tnnmg.preconditioner.patchsolver]
+tolerance          = 1e-10

dune/tectonic/tests/contacttest/staticcontacttest-3D.cfg

+# -*- mode:conf -*-
+[body0]
+smallestDiameter = 0.005  # 2e-3 [m]
+
+[body1]
+smallestDiameter = 0.005  # 2e-3 [m]
+
+[u0.solver]
+tolerance         = 1e-8
+
+[solver.tnnmg.preconditioner.basesolver]
+tolerance          = 1e-10
+
+[solver.tnnmg.preconditioner.patchsolver]
+tolerance          = 1e-10

dune/tectonic/tests/contacttest/staticcontacttest.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#define MY_DIM 2
+
+#include <atomic>
+#include <cmath>
+#include <csignal>
+#include <exception>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <memory>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+
+#include <dune/fufem/formatstring.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+
+#include <dune/tectonic/assemblers.hh>
+#include <dune/tectonic/gridselector.hh>
+#include <dune/tectonic/explicitgrid.hh>
+#include <dune/tectonic/explicitvectors.hh>
+
+#include <dune/tectonic/data-structures/enumparser.hh>
+#include <dune/tectonic/data-structures/enums.hh>
+#include <dune/tectonic/data-structures/network/contactnetwork.hh>
+#include <dune/tectonic/data-structures/matrices.hh>
+
+#include <dune/tectonic/factories/twoblocksfactory.hh>
+#include <dune/tectonic/factories/threeblocksfactory.hh>
+#include <dune/tectonic/factories/stackedblocksfactory.hh>
+
+#include <dune/tectonic/io/vtk.hh>
+
+#include <dune/tectonic/spatial-solving/tnnmg/functional.hh>
+#include <dune/tectonic/spatial-solving/tnnmg/zerononlinearity.hh>
+#include <dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh>
+#include <dune/tectonic/spatial-solving/solverfactory.hh>
+#include <dune/tectonic/spatial-solving/contact/nbodycontacttransfer.hh>
+
+#include <dune/tectonic/utils/debugutils.hh>
+
+// for getcwd
+#include <unistd.h>
+
+size_t const dims = MY_DIM;
+
+Dune::ParameterTree getParameters(int argc, char *argv[]) {
+  Dune::ParameterTree parset;
+  Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/dune/tectonic/tests/contacttest/staticcontacttest.cfg", parset);
+  Dune::ParameterTreeParser::readINITree(
+      Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/dune/tectonic/tests/contacttest/staticcontacttest-%dD.cfg", dims), parset);
+  Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+  return parset;
+}
+
+int main(int argc, char *argv[]) {
+  try {
+    Dune::MPIHelper::instance(argc, argv);
+
+    char buffer[256];
+    char *val = getcwd(buffer, sizeof(buffer));
+    if (val) {
+        std::cout << buffer << std::endl;
+        std::cout << argv[0] << std::endl;
+    }
+
+    std::ofstream out("staticcontacttest.log");
+    std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+    std::cout.rdbuf(out.rdbuf()); //redirect std::cout to log.txt
+
+    auto const parset = getParameters(argc, argv);
+
+    // ----------------------
+    // set up contact network
+    // ----------------------
+    using BlocksFactory = StackedBlocksFactory<Grid, Vector>;
+    BlocksFactory blocksFactory(parset);
+
+    using ContactNetwork = typename BlocksFactory::ContactNetwork;
+    blocksFactory.build();
+
+    auto& contactNetwork = blocksFactory.contactNetwork();
+
+    const size_t bodyCount = contactNetwork.nBodies();
+
+   /* for (size_t i=0; i<contactNetwork.nLevels(); i++) {
+        // printDofLocation(contactNetwork.body(i)->gridView());
+
+        //Vector def(contactNetwork.deformedGrids()[i]->size(dims));
+        //def = 1;
+        //deformedGridComplex.setDeformation(def, i);
+
+        const auto& level = *contactNetwork.level(i);
+
+        for (size_t j=0; j<level.nBodies(); j++) {
+            writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
+        }
+    }*/
+
+    for (size_t i=0; i<bodyCount; i++) {
+        writeToVTK(contactNetwork.body(i)->gridView(), "", "initial_body_" + std::to_string(i) + "_leaf");
+    }
+
+    // ----------------------------
+    // assemble contactNetwork
+    // ----------------------------
+    contactNetwork.assemble();
+
+    //printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
+
+
+    // ----------------------------
+    // compute minimal stress solution
+    // ----------------------------
+    using Matrix = typename ContactNetwork::Matrix;
+    const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+
+    // Initial displacement: Start from a situation of minimal stress,
+    // which is automatically attained in the case [v = 0 = a].
+    // Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
+    BitVector dirichletNodes;
+    contactNetwork.totalNodes("dirichlet", dirichletNodes);
+
+    // minimal stress displacement
+    std::vector<Vector> u(bodyCount);
+    for (size_t i=0; i<bodyCount; i++) {
+        u[i].resize(contactNetwork.body(i)->nVertices());
+        u[i] = 0.0;
+    }
+
+    // set rhs
+    std::vector<Vector> ell0(bodyCount);
+    for (size_t i=0; i<bodyCount; i++) {
+        ell0[i].resize(u[i].size());
+        ell0[i] = 0.0;
+
+        contactNetwork.body(i)->externalForce()(0.0, ell0[i]);
+    }
+
+    // set elasticity operator
+    const auto& matrices = contactNetwork.matrices().elasticity;
+    std::vector<const Matrix*> matrices_ptr(matrices.size());
+    for (size_t i=0; i<matrices_ptr.size(); i++) {
+        matrices_ptr[i] = matrices[i].get();
+    }
+
+    // assemble full global contact problem
+    Matrix bilinearForm;
+
+    nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+
+    Vector totalRhs;
+    nBodyAssembler.assembleRightHandSide(ell0, totalRhs);
+
+    Vector totalX;
+    nBodyAssembler.nodalToTransformed(u, totalX);
+
+    // get lower and upper obstacles
+    const auto& totalObstacles = nBodyAssembler.totalObstacles_;
+    Vector lower(totalObstacles.size());
+    Vector upper(totalObstacles.size());
+
+    for (size_t j=0; j<totalObstacles.size(); ++j) {
+        const auto& totalObstaclesj = totalObstacles[j];
+        auto& lowerj = lower[j];
+        auto& upperj = upper[j];
+        for (size_t d=0; d<dims; ++d) {
+            lowerj[d] = totalObstaclesj[d][0];
+            upperj[d] = totalObstaclesj[d][1];
+        }
+    }
+
+    // print problem
+    print(bilinearForm, "bilinearForm");
+    print(totalRhs, "totalRhs");
+    print(dirichletNodes, "ignore");
+
+    std::vector<const Dune::BitSetVector<1>*> frictionNodes;
+    contactNetwork.frictionNodes(frictionNodes);
+
+    for (size_t i=0; i<frictionNodes.size(); i++) {
+        print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
+    }
+
+    print(totalObstacles, "totalObstacles");
+    print(lower, "lower");
+    print(upper, "upper");
+
+    // set up functional
+    using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
+    Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper); //TODO
+
+       /* std::vector<BitVector> bodyDirichletNodes;
+        nBodyAssembler.postprocess(dirichletNodes, bodyDirichletNodes);
+        for (size_t i=0; i<bodyDirichletNodes.size(); i++) {
+          print(bodyDirichletNodes[i], "bodyDirichletNodes_" + std::to_string(i) + ": ");
+        }*/
+
+       /* print(bilinearForm, "matrix: ");
+        print(totalX, "totalX: ");
+        print(totalRhs, "totalRhs: ");*/
+
+    // make linear solver for linear correction in TNNMGStep
+    const auto& solverParset = parset.sub("u0.solver");
+
+    using Norm =  EnergyNorm<Matrix, Vector>;
+    using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+    // set multigrid solver
+    auto smoother = TruncatedBlockGSStep<Matrix, Vector>();
+
+    using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+    using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+    TransferOperators transfer(contactNetwork.nLevels()-1);
+    for (size_t i=0; i<transfer.size(); i++) {
+        transfer[i] = std::make_shared<TransferOperator>();
+        transfer[i]->setup(contactNetwork, i, i+1);
+    }
+
+    // Remove any recompute filed so that initially the full transferoperator is assembled
+    for (size_t i=0; i<transfer.size(); i++)
+        std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+    auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+    linearMultigridStep->setMGType(1, 3, 3);
+    linearMultigridStep->setSmoother(smoother);
+    linearMultigridStep->setTransferOperators(transfer);
+
+    Norm norm(*linearMultigridStep);
+
+    auto linearSolver = std::make_shared<LinearSolver>(linearMultigridStep, parset.get<int>("solver.tnnmg.main.multi"), parset.get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
+
+    // set up TNMMG solver
+    using Factory = SolverFactory<Functional, BitVector>;
+    Factory factory(parset.sub("solver.tnnmg"), J, dirichletNodes);
+
+    factory.build(linearSolver);
+    auto tnnmgStep = factory.step();
+    factory.setProblem(totalX);
+
+    //const EnergyNorm<Matrix, Vector> norm(bilinearForm);
+
+    std::cout << "solving linear problem for u..." << std::endl;
+    LoopSolver<Vector> solver(
+            *tnnmgStep.get(), solverParset.get<size_t>("maximumIterations"),
+            solverParset.get<double>("tolerance"), norm,
+            solverParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
+
+    solver.preprocess();
+    solver.solve();
+
+    nBodyAssembler.postprocess(tnnmgStep->getSol(), u);
+
+    print(tnnmgStep->getSol(), "totalX:");
+    print(u, "u:");
+
+    for (size_t i=0; i<bodyCount; i++) {
+      auto& body = contactNetwork.body(i);
+      body->setDeformation(u[i]);
+      writeToVTK(body->gridView(), "", "body_" + std::to_string(i) + "_leaf");
+    }
+
+    /*
+    using BoundaryFunctions = typename ContactNetwork::BoundaryFunctions;
+    using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
+    using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
+                               StateUpdater<ScalarVector, Vector>>;
+
+    BoundaryFunctions velocityDirichletFunctions;
+    contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
+
+    BoundaryNodes dirichletNodes;
+    contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
+
+    for (size_t i=0; i<dirichletNodes.size(); i++) {
+        for (size_t j=0; j<dirichletNodes[i].size(); j++) {
+        print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
+        }
+    }*/
+
+    std::cout.rdbuf(coutbuf); //reset to standard output again
+
+  } catch (Dune::Exception &e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+  } catch (std::exception &e) {
+    std::cerr << "Standard exception: " << e.what() << std::endl;
+  }
+}

dune/tectonic/tests/contacttest/staticcontacttest.cfg

+# -*- mode:conf -*-
+gravity         = 9.81     # [m/s^2]
+
+[body]
+length          = 0.4      # [m]
+height          = 0.04     # [m]
+depth           = 0.04     # [m]
+bulkModulus     = 2190    # [Pa] #2190
+poissonRatio    = 0.11     # [1]  #0.11
+[body.elastic]
+density         = 750      # [kg/m^3] #750
+shearViscosity  = 0     # [Pas]
+bulkViscosity   = 0     # [Pas]
+[body.viscoelastic]
+density         = 1000     # [kg/m^3]
+shearViscosity  = 0     # [Pas]
+bulkViscosity   = 0     # [Pas]
+
+[boundary.friction]
+C               = 6       # [Pa]
+mu0             = 0.48      # [ ]
+V0              = 1e-3     # [m/s]
+L               = 1e-6  # [m]
+initialAlpha    = 0        # [ ]
+stateModel      = AgeingLaw
+frictionModel   = Truncated #Regularised
+[boundary.friction.weakening]
+a               = 0.054    # [ ]
+b               = 0.074    # [ ]
+[boundary.friction.strengthening]
+a               = 0.054    # [ ]
+b               = 0.074    # [ ]
+
+[boundary.neumann]
+sigmaN          = 0.0      # [Pa]
+
+[boundary.dirichlet]
+finalVelocity   = 5e-5     # [m/s]
+
+[problem]
+bodyCount       = 2
+
+[u0.solver]
+maximumIterations = 100
+verbosity         = full
+
+[solver.tnnmg.preconditioner]
+mode         = additive
+patchDepth   = 1
+maximumIterations = 2
+verbosity         = quiet
+[solver.tnnmg.preconditioner.patchsolver]
+maximumIterations = 100
+verbosity         = quiet
+[solver.tnnmg.preconditioner.basesolver]
+maximumIterations = 10000
+verbosity         = quiet
+
+[solver.tnnmg.main]
+pre   = 1
+multi = 5 # number of multigrid steps
+post  = 0
+
+
+

dune/tectonic/tests/couplingtest.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_GRIDGLUE_TEST_COUPLINGTEST_HH
+#define DUNE_GRIDGLUE_TEST_COUPLINGTEST_HH
+
+#include <iostream>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/version.hh>
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <dune/grid-glue/gridglue.hh>
+
+template <class IntersectionIt>
+bool testIntersection(const IntersectionIt & rIIt)
+{
+  bool success = true;
+
+  typedef typename IntersectionIt::value_type Intersection;
+  // Dimension of the intersection
+  const int dim = Intersection::mydim;
+
+  // Dimension of world coordinates
+  const int coorddim = Intersection::coorddim;
+
+  // Create a set of test points
+  const Dune::QuadratureRule<double, dim>& quad = Dune::QuadratureRules<double, dim>::rule(rIIt->type(), 3);
+
+  for (unsigned int l=0; l<quad.size(); l++) {
+    const auto inside = rIIt->inside();
+    const auto outside = rIIt->outside();
+
+    Dune::FieldVector<double, dim> quadPos = quad[l].position();
+
+    Dune::FieldVector<double, Intersection::InsideGridView::dimensionworld> localGrid0Pos =
+      inside.geometry().global(rIIt->geometryInInside().global(quadPos));
+
+    // currently the intersection maps to the GV::dimworld, this will hopefully change soon
+    Dune::FieldVector<double, Intersection::InsideGridView::dimensionworld> globalGrid0Pos =
+      rIIt->geometry().global(quadPos);
+
+    Dune::FieldVector<double, Intersection::OutsideGridView::dimensionworld> localGrid1Pos =
+      outside.geometry().global(rIIt->geometryInOutside().global(quadPos));
+
+    // currently the intersection maps to the GV::dimworld, this will hopefully change soon
+    Dune::FieldVector<double, Intersection::OutsideGridView::dimensionworld> globalGrid1Pos =
+      rIIt->geometryOutside().global(quadPos);
+
+    // Test whether local grid0 position is consistent with global grid0 position
+    if ( (localGrid0Pos-globalGrid0Pos).two_norm() >= 1e-6 )
+    {
+      std::cout << __FILE__ << ":" << __LINE__ << ": error: assert( (localGrid0Pos-globalGrid0Pos).two_norm() < 1e-6 ) failed\n";
+      std::cerr << "localGrid0Pos  = " << localGrid0Pos << "\n";
+      std::cerr << "globalGrid0Pos = " << globalGrid0Pos << "\n";
+      success = false;
+    }
+
+    // Test whether local grid1 position is consistent with global grid1 position
+    if ( (localGrid1Pos-globalGrid1Pos).two_norm() >= 1e-6 )
+    {
+      std::cout << __FILE__ << ":" << __LINE__ << ": error: assert( (localGrid1Pos-globalGrid1Pos).two_norm() < 1e-6 ) failed\n";
+      std::cerr << "localGrid1Pos  = " << localGrid1Pos << "\n";
+      std::cerr << "globalGrid1Pos = " << globalGrid1Pos << "\n";
+      success = false;
+    }
+
+    // Here we assume that the two interfaces match geometrically:
+    if ( (globalGrid0Pos-globalGrid1Pos).two_norm() >= 1e-4 )
+    {
+      std::cout << __FILE__ << ":" << __LINE__ << ": error: assert( (globalGrid0Pos-globalGrid1Pos).two_norm() < 1e-4 ) failed\n";
+      std::cerr << "localGrid0Pos  = " << localGrid0Pos << "\n";
+      std::cerr << "globalGrid0Pos = " << globalGrid0Pos << "\n";
+      std::cerr << "localGrid1Pos  = " << localGrid1Pos << "\n";
+      std::cerr << "globalGrid1Pos = " << globalGrid1Pos << "\n";
+      success = false;
+    }
+
+    // Test the normal vector methods.  At least test whether they don't crash
+    if (coorddim - dim == 1) // only test for codim 1
+    {
+      rIIt->outerNormal(quadPos);
+      rIIt->unitOuterNormal(quadPos);
+      rIIt->integrationOuterNormal(quadPos);
+      rIIt->centerUnitOuterNormal();
+    }
+  }
+
+  return success;
+}
+
+
+template <class GlueType>
+void testCoupling(const GlueType& glue)
+{
+  bool success = true;
+#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 6)
+  using View0Mapper = Dune::MultipleCodimMultipleGeomTypeMapper< typename GlueType::template GridView<0> >;
+  using View1Mapper = Dune::MultipleCodimMultipleGeomTypeMapper< typename GlueType::template GridView<1> >;
+  View0Mapper view0mapper(glue.template gridView<0>(), Dune::mcmgElementLayout());
+  View1Mapper view1mapper(glue.template gridView<1>(), Dune::mcmgElementLayout());
+#else
+  using View0Mapper = Dune::MultipleCodimMultipleGeomTypeMapper< typename GlueType::template GridView<0>, Dune::MCMGElementLayout >;
+  using View1Mapper = Dune::MultipleCodimMultipleGeomTypeMapper< typename GlueType::template GridView<1>, Dune::MCMGElementLayout >;
+  View0Mapper view0mapper(glue.template gridView<0>());
+  View1Mapper view1mapper(glue.template gridView<1>());
+#endif
+
+  std::vector<unsigned int> countInside0(view0mapper.size());
+  std::vector<unsigned int> countOutside1(view1mapper.size());
+  std::vector<unsigned int> countInside1(view1mapper.size(), 0);
+  std::vector<unsigned int> countOutside0(view0mapper.size(), 0);
+
+  // ///////////////////////////////////////
+  //   IndexSet
+  // ///////////////////////////////////////
+
+  {
+    size_t count = 0;
+    for (auto rIIt = glue.template ibegin<0>(); rIIt != glue.template iend<0>(); ++rIIt) count ++;
+    typename GlueType::IndexSet is = glue.indexSet();
+    if(is.size() != glue.size())
+      DUNE_THROW(Dune::Exception,
+        "Inconsistent size information: indexSet.size() " << is.size() << " != GridGlue.size() " << glue.size());
+    if(is.size() != count)
+      DUNE_THROW(Dune::Exception,
+        "Inconsistent size information: indexSet.size() " << is.size() << " != iterator count " << count);
+    std::vector<bool> visited(count, false);
+    for (auto rIIt = glue.template ibegin<0>(); rIIt != glue.template iend<0>(); ++rIIt) {
+      size_t idx = is.index(*rIIt);
+      if(idx >= count)
+        DUNE_THROW(Dune::Exception,
+          "Inconsistent IndexSet: index " << idx << " out of range, size is " << count);
+      if(visited[idx] != false)
+        DUNE_THROW(Dune::Exception,
+          "Inconsistent IndexSet: visited index " << idx << " twice");
+      visited[idx] = true;
+    }
+    // make sure that we have a consecutive zero starting index set
+    for (size_t i = 0; i<count; i++)
+    {
+      if (visited[i] != true)
+        DUNE_THROW(Dune::Exception,
+          "Non-consective IndexSet: " << i << " missing.");
+    }
+  }
+
+
+  // ///////////////////////////////////////
+  //   MergedGrid centric Grid0->Grid1
+  // ///////////////////////////////////////
+
+  {
+    for (auto rIIt = glue.template ibegin<0>(); rIIt != glue.template iend<0>(); ++rIIt)
+    {
+      if (rIIt->self() && rIIt->neighbor())
+      {
+        const auto index0 = view0mapper.index(rIIt->inside());
+        const auto index1 = view1mapper.index(rIIt->outside());
+
+        countInside0[index0]++;
+        countOutside1[index1]++;
+        success = success && testIntersection(rIIt);
+      }
+    }
+  }
+
+  // ///////////////////////////////////////
+  //   MergedGrid centric Grid1->Grid0
+  // ///////////////////////////////////////
+
+  {
+    for (auto rIIt = glue.template ibegin<1>(); rIIt != glue.template iend<1>(); ++rIIt)
+    {
+      if (rIIt->self() && rIIt->neighbor())
+      {
+        const auto index1 = view1mapper.index(rIIt->inside());
+        const auto index0 = view0mapper.index(rIIt->outside());
+
+        countInside1[index1]++;
+        countOutside0[index0]++;
+        success = success && testIntersection(rIIt);
+      }
+    }
+  }
+
+  if (! success)
+    DUNE_THROW(Dune::Exception, "Test failed, see above for details.");
+}
+
+#endif // DUNE_GRIDGLUE_TEST_COUPLINGTEST_HH

dune/tectonic/tests/globalfrictioncontainertest.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <iostream>
+#include <vector>
+
+#include <exception>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/stdstreams.hh>
+
+#include "../data-structures/globalfrictioncontainer.hh"
+
+struct Body {
+    size_t id;
+    size_t localAlpha;
+
+    void updateAlpha(size_t lAlpha) {
+        localAlpha = lAlpha;
+    }
+};
+
+int main(int argc, char *argv[]) {
+    try {
+        Dune::MPIHelper::instance(argc, argv);
+
+        std::cout << "--------------------------" << std::endl;
+        std::cout << "-- GlobalFriction Test: --" << std::endl;
+        std::cout << "--------------------------" << std::endl << std::endl;
+
+        // have to be <10
+        size_t n = 5;
+        size_t m = 3;
+        size_t l = 2;
+
+        std::vector<std::vector<int>> ref;
+        ref.resize(n);
+        for (size_t i=0; i<n; i++) {
+            ref[i].resize(m);
+            for (size_t j=0; j<m; j++) {
+                ref[i][j] = 10*i + j;
+            }
+        }
+
+        // resize() test
+        bool resizeTest = true;
+        {
+            GlobalFrictionContainer<Body,1> globalFriction1;
+            globalFriction1.resize({n});
+            resizeTest = resizeTest & (globalFriction1.size() == n);
+
+            GlobalFrictionContainer<Body,2> globalFriction2;
+            globalFriction2.resize({n});
+            resizeTest = resizeTest & (globalFriction2.size() == n);
+
+            for (size_t i=0; i<globalFriction2.size(); i++) {
+                resizeTest = resizeTest & (globalFriction2[i].size() == 0);
+            }
+        }
+
+        // operator[] test
+        bool operatorTest = true;
+        {
+            GlobalFrictionContainer<Body,2> globalFriction;
+            globalFriction.resize({n, m});
+            operatorTest = operatorTest & (globalFriction.size() == n);
+
+            for (size_t i=0; i<n; i++) {
+                for (size_t j=0; j<m; j++) {
+                    operatorTest = operatorTest & (globalFriction[i][j] == nullptr);
+                }
+            }
+        }
+
+        // updateAlpha() test
+        bool updateTest = true;
+        {
+            GlobalFrictionContainer<Body,3> globalFriction;
+            globalFriction.resize({n,m,l});
+
+            for (size_t i=0; i<n; i++) {
+                for (size_t j=0; j<m; j++) {
+                    for (size_t k=0; k<l; k++) {
+                        globalFriction[i][j][k] = std::make_shared<Body>();
+                        globalFriction[i][j][k]->id = 10*i + j;
+                    }
+                }
+            }
+
+            globalFriction.updateAlpha(ref);
+
+            for (size_t i=0; i<n; i++) {
+                for (size_t j=0; j<m; j++) {
+                    for (size_t k=0; k<l; k++) {
+                        updateTest = updateTest & (globalFriction[i][j][k]->id == globalFriction[i][j][k]->localAlpha);
+                    }
+                }
+            }
+
+        }
+
+        std::cout << "resize(): " << resizeTest << std::endl;
+        std::cout << "operator[]: " << operatorTest << std::endl;
+        std::cout << "updateAlpha(): " << updateTest << std::endl;
+        std::cout << "-------------------------" << std::endl << std::endl;
+
+        bool passed = resizeTest & operatorTest & updateTest;
+
+        std::cout << "Overall the test " << (passed ? "was successful!" : "failed!") << std::endl;
+
+        return passed ? 0 : 1;
+
+    } catch (Dune::Exception &e) {
+        Dune::derr << "Dune reported error: " << e << std::endl;
+    } catch (std::exception &e) {
+        std::cerr << "Standard exception: " << e.what() << std::endl;
+    }
+}

dune/tectonic/tests/gridgluefrictiontest.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <exception>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/stdstreams.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/bitsetvector.hh>
+
+//#include <dune/fufem/assemblers/assembler.hh>
+#include <dune/fufem/functiontools/basisinterpolator.hh>
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
+#include <dune/fufem/functions/constantfunction.hh>
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+
+#include <dune/grid/uggrid.hh>
+
+#include <dune/grid-glue/adapter/gridgluevtkwriter.hh>
+#include <dune/grid-glue/merging/merger.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/assemblers/dualmortarcoupling.hh>
+#include <dune/contact/assemblers/nbodyassembler.hh>
+#include <dune/contact/common/dualbasisadapter.hh>
+
+#include "../utils/debugutils.hh"
+
+//#include <dune/tectonic/transformedglobalratestatefriction.hh>
+
+#include "common.hh"
+
+const int dim = 2;
+const int n = 5;
+const bool simplexGrid = true;
+
+const std::string path = "";
+const std::string outputFile = "gridgluefrictiontest.log";
+
+#if HAVE_UG
+    using GridType = typename Dune::UGGrid<dim>;
+#else
+#error No UG!
+#endif
+    using LevelView = typename GridType::LevelGridView;
+    using LeafView = typename GridType::LeafGridView;
+    using LevelBoundaryPatch = BoundaryPatch<LevelView>;
+    using LeafBoundaryPatch = BoundaryPatch<LeafView>;
+
+    using c_type = double;
+    using GlobalCoords = Dune::FieldVector<c_type, dim>;
+    using Vector = Dune::FieldVector<c_type, 1>;
+    using BlockVector = Dune::BlockVector<Vector>;
+
+    using CouplingPair = Dune::Contact::CouplingPair<GridType, GridType, c_type>;
+    using CouplingType = Dune::Contact::DualMortarCoupling<c_type, GridType>;
+    using NBodyAssembler = Dune::Contact::NBodyAssembler<GridType, BlockVector>;
+
+
+
+// boundary function along y=1.0 line
+template <class VectorType>
+class BoundaryFunction0 : public Dune::VirtualFunction<VectorType, c_type> {
+public:
+    void evaluate(const VectorType& x, c_type& res) const override {
+        if (isClose(x[1], 1.0)) {
+            res = 1.0 ;// +std::sin(x[0]);
+        } else {
+            res = 1.0;
+        }
+    }
+};
+
+// boundary function along y=1.0 line
+template <class VectorType>
+class BoundaryFunction1 : public Dune::VirtualFunction<VectorType, c_type> {
+public:
+    void evaluate(const VectorType& x, c_type& res) const override{
+        if (isClose(x[1], 1.0)) {
+            res = 2.0 ;//+ std::cos(x[0]);
+        } else {
+            res = 2.0;
+        }
+    }
+};
+/*
+template <class GridView>
+void dualWeights(const BoundaryPatch<GridView>& boundaryPatch, BlockVector& weights, bool initializeVector = true) {
+
+    typedef typename BoundaryPatch<GridView>::iterator BoundaryIterator;
+
+    /*typedef typename LocalBoundaryFunctionalAssemblerType::LocalVector LocalVector;
+    typedef typename TrialBasis::LinearCombination LinearCombination;
+
+    int rows = tBasis_.size();
+
+
+        const auto inside = it->inside();
+
+        // get shape functions
+        DualP1LocalFiniteElement<class D, class R, dim>
+  class
+        const typename TrialBasis::LocalFiniteElement& lFE = tBasis_.getLocalFiniteElement(inside);
+
+        LocalVector localB(tFE.localBasis().size());
+
+    }*/
+/*
+    // cache for the dual functions on the boundary
+    using DualCache = Dune::Contact::DualBasisAdapter<GridView, field_type>;
+    std::unique_ptr<DualCache> dualCache;
+    dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView, field_type> >();
+
+    if (initializeVector) {
+        weights.resize(rows);
+        std::fill(weights.begin(), weights.end(), 0.0);
+    }
+
+    // loop over all boundary intersections
+    BoundaryIterator it = boundaryPatch.begin();
+    BoundaryIterator end = boundaryPatch.end();
+
+    for (; it != end; ++it) {
+        const auto& inside = it->inside();
+
+        if (!boundaryPatch.contains(inside, it->indexInInside()))
+            continue;
+
+        // types of the elements supporting the boundary segments in question
+        Dune::GeometryType nonmortarEType = inside.type();
+
+        const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(nonmortarEType);
+
+        int noOfMortarVec = targetRefElement.size(dim);
+
+        Dune::GeometryType nmFaceType = domainRefElement.type(rIs.indexInInside(),1);
+
+
+        // Select a quadrature rule
+        // 2 in 2d and for integration over triangles in 3d.  If one (or both) of the two faces involved
+        // are quadrilaterals, then the quad order has to be risen to 3 (4).
+        int quadOrder = 2 + (!nmFaceType.isSimplex());
+        const auto& quadRule = Dune::QuadratureRules<ctype, dim-1>::rule(it->type(), quadOrder);
+
+        dualCache->bind(inside, it->indexInInside());
+
+
+
+        const auto& rGeomInInside = it->geometryInInside();
+
+        int nNonmortarFaceNodes = domainRefElement.size(it->indexInInside(),1,dim);
+        std::vector<int> nonmortarFaceNodes;
+        for (int i=0; i<nNonmortarFaceNodes; i++) {
+          int faceIdxi = domainRefElement.subEntity(it->indexInInside(), 1, i, dim);
+          nonmortarFaceNodes.push_back(faceIdxi);
+        }
+
+        // store values of shape functions
+        std::vector<Dune::FieldVector<field_type,1>> values(tFE.localBasis().size());
+
+        // get geometry and store it
+        const auto& rGeom = it->geometry();
+
+        localVector = 0.0;
+
+        // get quadrature rule
+        QuadratureRuleKey tFEquad(tFE);
+        QuadratureRuleKey quadKey = tFEquad.product(functionQuadKey_);
+        const auto& quad = QuadratureRuleCache<double, dim>::rule(quadKey);
+
+        // loop over quadrature points
+        for (const auto& pt : quadRule) {
+
+            // get quadrature point
+            const Dune::FieldVector<ctype,dim>& quadPos = quad[pt].position();
+
+            // get integration factor
+            const ctype integrationElement = rGeom.integrationElement(quadPos);
+
+            // evaluate basis functions
+            dualCache->evaluateFunction(quadPos, values);
+
+            // compute values of function
+            T f_pos;
+            const GridFunction* gf = dynamic_cast<const GridFunction*>(&f_);
+            if (gf and gf->isDefinedOn(element))
+                gf->evaluateLocal(element, quadPos, f_pos);
+            else
+                f_.evaluate(geometry.global(quadPos), f_pos);
+
+            // and vector entries
+            for(size_t i=0; i<values.size(); ++i)
+            {
+                localVector[i].axpy(values[i]*quad[pt].weight()*integrationElement, f_pos);
+            }
+        }
+
+        for (size_t i=0; i<tFE.localBasis().size(); ++i) {
+            int idx = tBasis_.index(inside, i);
+            b[idx] += localB[i];
+        }
+}*/
+
+int main(int argc, char *argv[]) { try {
+    Dune::MPIHelper::instance(argc, argv);
+
+    std::ofstream out(path + outputFile);
+    std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+    std::cout.rdbuf(out.rdbuf()); //redirect std::cout to outputFile
+
+    std::cout << "------------------------------" << std::endl;
+    std::cout << "-- Grid-Glue Friction Test: --" << std::endl;
+    std::cout << "------------------------------" << std::endl << std::endl;
+
+    // building grids
+    std::vector<std::shared_ptr<GridType>> grids(2);
+
+    GlobalCoords lowerLeft0({0, 0});
+    GlobalCoords upperRight0({2, 1});
+    buildGrid(lowerLeft0, upperRight0, n, grids[0]);
+
+    GlobalCoords lowerLeft1({0, 1});
+    GlobalCoords upperRight1({2, 2});
+    buildGrid(lowerLeft1, upperRight1, n, grids[1]);
+
+    // writing grids
+    for (size_t i=0; i<grids.size(); i++) {
+        const auto& levelView = grids[i]->levelGridView(0);
+        writeToVTK(levelView, path, "body_" + std::to_string(i) + "_level0");
+    }
+
+    // compute coupling boundaries
+    LevelView gridView0 = grids[0]->levelGridView(0);
+    LevelView gridView1 = grids[1]->levelGridView(0);
+    LevelBoundaryPatch upper(gridView0);
+    LevelBoundaryPatch lower(gridView1);
+
+    lower.insertFacesByProperty([&](typename LevelView::Intersection const &in) {
+        return xyBetween(lowerLeft1, {upperRight1[0], lowerLeft1[1]}, in.geometry().center());
+    });
+
+    upper.insertFacesByProperty([&](typename LevelView::Intersection const &in) {
+        return xyBetween({lowerLeft0[0], upperRight0[1]}, upperRight0, in.geometry().center());
+    });
+
+    // set contact coupling
+    Dune::Contact::NormalProjection<LeafBoundaryPatch> contactProjection;
+
+    CouplingPair coupling;
+    coupling.set(0, 1, upper, lower, 0.1, CouplingPair::CouplingType::STICK_SLIP, contactProjection, nullptr);
+
+   /* double coveredArea_ = 0.8;
+    CouplingType contactCoupling;
+    contactCoupling.setGrids(*grids[0], *grids[1]);
+    contactCoupling.setupContactPatch(*coupling.patch0(),*coupling.patch1());
+    contactCoupling.gridGlueBackend_ = coupling.backend();
+    contactCoupling.setCoveredArea(coveredArea_);
+    contactCoupling.setup();*/
+
+    // set nBodyAssembler
+    using NBodyAssembler = Dune::Contact::NBodyAssembler<GridType, BlockVector>;
+    NBodyAssembler nBodyAssembler(grids.size(), 1);
+
+    std::vector<const GridType*> grids_ptr(grids.size());
+    for (size_t i=0; i<grids_ptr.size(); i++) {
+        grids_ptr[i] = grids[i].get();
+    }
+
+    nBodyAssembler.setGrids(grids_ptr);
+    nBodyAssembler.setCoupling(coupling, 0);
+
+    nBodyAssembler.assembleTransferOperator();
+    nBodyAssembler.assembleObstacle();
+
+    // define basis
+    using Basis = P1NodalBasis<LevelView, double>;
+    Basis basis0(grids[0]->levelGridView(0));
+    Basis basis1(grids[1]->levelGridView(0));
+
+    std::cout << "--------------" << std::endl;
+    std::cout << "--- Basis0 ---" << std::endl;
+    std::cout << "--------------" << std::endl;
+    printBasisDofLocation(basis0);
+
+    std::cout << "--------------" << std::endl;
+    std::cout << "--- Basis1 ---" << std::endl;
+    std::cout << "--------------" << std::endl;
+    printBasisDofLocation(basis1);
+
+    // set grid functions on coupling boundary
+    std::vector<BlockVector> f(2);
+
+    BasisInterpolator<Basis, BlockVector, BoundaryFunction0<GlobalCoords>> interpolator0;
+    interpolator0.interpolate(basis0, f[0], BoundaryFunction0<GlobalCoords>());
+
+    BasisInterpolator<Basis, BlockVector, BoundaryFunction1<GlobalCoords>> interpolator1;
+    interpolator1.interpolate(basis1, f[1], BoundaryFunction1<GlobalCoords>());
+
+    BlockVector transformedF;
+    nBodyAssembler.nodalToTransformed(f, transformedF);
+
+    std::vector<Dune::BitSetVector<1>> boundaryVertices(2);
+    const auto& mortarCoupling = nBodyAssembler.getContactCouplings()[0];
+    const auto& nonmortarBoundary = mortarCoupling->nonmortarBoundary();
+    const auto& mortarBoundary = mortarCoupling->mortarBoundary();
+
+    nonmortarBoundary.getVertices(boundaryVertices[0]);
+    mortarBoundary.getVertices(boundaryVertices[1]);
+
+    print(f[0], "f0: ");
+    print(boundaryVertices[0], "nonmortarBoundary: ");
+    print(f[1], "f1: ");
+    print(boundaryVertices[1], "mortarBoundary: ");
+    print(transformedF, "transformedF: ");
+    writeToVTK(basis0, f[0], path, "body_0_level0");
+    writeToVTK(basis1, f[1], path, "body_1_level0");
+
+    print(mortarCoupling->mortarLagrangeMatrix(), "M: ");
+    print(mortarCoupling->nonmortarLagrangeMatrix(), "D: ");
+
+    std::vector<BlockVector> postprocessed(2);
+    nBodyAssembler.postprocess(transformedF, postprocessed);
+
+    print(postprocessed, "postprocessed: ");
+
+    const auto& contactCouplings = nBodyAssembler.getContactCouplings();
+    const auto& contactCoupling = contactCouplings[0];
+    const auto& glue = *contactCoupling->getGlue();
+
+    size_t isCount = 0;
+    for (const auto& rIs : intersections(glue)) {
+        std::cout << "intersection id: " << isCount
+                  << " insideElement: " << gridView0.indexSet().index(rIs.inside())
+                  << " outsideElement: " << gridView1.indexSet().index(rIs.outside()) << std::endl;
+        isCount++;
+    }
+
+    for (size_t i=0; i<isCount; i++) {
+        const auto& is = glue.getIntersection(i);
+        std::cout << "intersection id: " << i
+                  << " insideElement: " << gridView0.indexSet().index(is.inside())
+                  << " outsideElement: " << gridView1.indexSet().index(is.outside()) << std::endl;
+    }
+
+   /* using DualBasis = ;
+    DualBasis dualBasis;
+    BoundaryFunctionalAssembler<DualBasis> bAssembler(dualBasis, nonmortarBoundary);
+
+    ConstantFunction<GlobalCoords, Vector> oneFunction(1);
+
+    BlockVector b;
+    L2FunctionalAssembler localAssembler(oneFunction);
+    bAssembler.assemble(localAssembler, b);
+
+    print(b, "b: ");*/
+
+    /*
+    std::vector<BlockVector> g(2);
+    g[0].resize(f[0].size());
+    g[1].resize(f[1].size());
+
+    g[1][6] = 2;
+    BlockVector transformedG;
+    nBodyAssembler.nodalToTransformed(g, transformedG);
+
+    print(g[1], "g1: ");
+    print(transformedG, "transformedG: ");
+
+    // merged gridGlue coupling boundary
+    auto& gridGlue = *contactCoupling.getGlue();
+
+    // make basis grid functions
+    auto&& gridFunction0 = Functions::makeFunction(basis0, f[0]); */
+
+   /* for (const auto& bIs : intersections(upper)) {
+        const auto& inside = bIs.inside();
+        const auto& bGeometry = bIs.geometry();
+
+        for (size_t i=0; i<bGeometry.corners(); i++) {
+
+
+            typename BasisGridFunction<Basis, BlockVector>::RangeType y;
+            gridFunction1.evaluateLocal(outside, outGeometry.local(rGeometry.corner(i)), y);
+            print(rGeometry.corner(i), "corner " + std::to_string(i));
+            std::cout << "f1(corner) = " << y << std::endl;
+        }
+    }*/
+
+    /*
+    auto&& gridFunction1 = Functions::makeFunction(basis1, f[1]);
+
+    for (const auto& rIs : intersections(gridGlue)) {
+        const auto& inside = rIs.inside();
+        const auto& outside = rIs.outside();
+
+        const auto& rGeometry = rIs.geometry();
+        const auto& outGeometry = outside.geometry();
+
+        for (size_t i=0; i<rGeometry.corners(); i++) {
+            typename BasisGridFunction<Basis, BlockVector>::RangeType y;
+            gridFunction1.evaluateLocal(outside, outGeometry.local(rGeometry.corner(i)), y);
+            print(rGeometry.corner(i), "corner " + std::to_string(i));
+            std::cout << "f1(corner) = " << y << std::endl;
+
+            std::cout << std::endl;
+        }
+        std::cout << "---------"<< std::endl;
+    }*/
+
+    for (const auto& elem : elements(gridView0)) {
+        std::cout << "seed element corners:" << std::endl;
+
+        const auto& eRefElement = Dune::ReferenceElements<double, dim>::general(elem.type());
+        size_t vertices = eRefElement.size(dim);
+
+        for (size_t j=0; j<vertices; j++) {
+            print(elem.geometry().corner(j), "corner: ");
+        }
+
+
+        for (auto&& is : intersections(gridView0, elem)) {
+            const auto& isGeo = is.geometry();
+
+            if (!is.neighbor())
+                continue;
+
+            std::cout << "neighbor corners:" << std::endl;
+
+            const auto& outside = is.outside();
+
+            const auto& refElement = Dune::ReferenceElements<double, dim>::general(outside.type());
+            size_t nVertices = refElement.size(dim);
+
+            const auto& isRefElement = Dune::ReferenceElements<double, dim-1>::general(isGeo.type());
+
+            for (size_t j=0; j<nVertices; j++) {
+                print(outside.geometry().corner(j), "corner: ");
+                const auto& local = elem.geometry().local(outside.geometry().corner(j));
+
+                print(local, "local:");
+            }
+
+        }
+
+        break;
+    }
+
+
+    bool passed = true;
+
+    std::cout << "Overall, the test " << (passed ? "was successful!" : "failed!") << std::endl;
+
+    std::cout.rdbuf(coutbuf); //reset to standard output again
+    return passed ? 0 : 1;
+
+} catch (Dune::Exception &e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+} catch (std::exception &e) {
+    std::cerr << "Standard exception: " << e.what() << std::endl;
+} // end try
+} // end main

dune/tectonic/tests/hdf5test/CMakeLists.txt

+add_custom_target(tectonic_tests_hdf5test SOURCES
+  hdf5test.cfg
+  hdf5test-2D.cfg
+)
+
+set(HDF5TEST_SOURCE_FILES
+  ../../assemblers.cc
+  ../../data-structures/body/body.cc
+  ../../data-structures/network/levelcontactnetwork.cc
+  ../../data-structures/network/contactnetwork.cc
+  ../../data-structures/enumparser.cc
+  ../../factories/strikeslipfactory.cc
+  #../../io/vtk.cc
+  #../../io/hdf5/frictionalboundary-writer.cc
+  #../../io/hdf5/iteration-writer.cc
+  #../../io/hdf5/time-writer.cc
+  ../../problem-data/grid/simplexmanager.cc
+  ../../problem-data/grid/trianglegeometry.cc
+  ../../problem-data/grid/trianglegrids.cc
+  ../../spatial-solving/solverfactory.cc
+  hdf5test.cc
+)
+
+foreach(_dim 2)
+  set(_hdf5test_target hdf5test-${_dim}D)
+
+  dune_add_test(NAME ${_hdf5test_target} SOURCES ${HDF5TEST_SOURCE_FILES})
+
+  add_dune_ug_flags(${_hdf5test_target})
+  add_dune_hdf5_flags(${_hdf5test_target})
+
+  set_property(TARGET ${_hdf5test_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+endforeach()

dune/tectonic/tests/hdf5test/hdf5test-2D.cfg

+# -*- mode:conf -*-
+[body0]
+smallestDiameter = 0.3  # 2e-3 [m]
+
+[body1]
+smallestDiameter = 0.3 # 2e-3 [m]
+
+[timeSteps]
+refinementTolerance = 1e-5 # 1e-5
+
+[u0.solver]
+tolerance         = 1e-8
+
+[a0.solver]
+tolerance         = 1e-8
+
+[v.solver]
+tolerance         = 1e-8
+
+[v.fpi]
+tolerance         = 1e-5
+
+[solver.tnnmg.preconditioner.basesolver]
+tolerance          = 1e-10
+
+[solver.tnnmg.preconditioner.patchsolver]
+tolerance          = 1e-10

dune/tectonic/tests/hdf5test/hdf5test.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#define MY_DIM 2
+
+#include <atomic>
+#include <cmath>
+#include <csignal>
+#include <exception>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <memory>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/fufem/formatstring.hh>
+#include <dune/fufem/hdf5/file.hh>
+
+#include <dune/tectonic/assemblers.hh>
+#include <dune/tectonic/gridselector.hh>
+#include <dune/tectonic/explicitgrid.hh>
+#include <dune/tectonic/explicitvectors.hh>
+
+#include <dune/tectonic/data-structures/enumparser.hh>
+#include <dune/tectonic/data-structures/enums.hh>
+#include <dune/tectonic/data-structures/network/contactnetwork.hh>
+#include <dune/tectonic/data-structures/program_state.hh>
+
+#include <dune/tectonic/factories/strikeslipfactory.hh>
+
+#include <dune/tectonic/io/vtk.hh>
+#include <dune/tectonic/io/hdf5-writer.hh>
+
+#include <dune/tectonic/utils/geocoordinate.hh>
+#include <dune/tectonic/utils/debugutils.hh>
+
+// for getcwd
+#include <unistd.h>
+
+size_t const dims = MY_DIM;
+
+Dune::ParameterTree getParameters(int argc, char *argv[]) {
+  Dune::ParameterTree parset;
+  Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/dune/tectonic/tests/hdf5test/hdf5test.cfg", parset);
+  Dune::ParameterTreeParser::readINITree(
+      Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/dune/tectonic/tests/hdf5test/hdf5test-%dD.cfg", dims), parset);
+  Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+  return parset;
+}
+
+int main(int argc, char *argv[]) {
+  try {
+    Dune::MPIHelper::instance(argc, argv);
+
+    char buffer[256];
+    char *val = getcwd(buffer, sizeof(buffer));
+    if (val) {
+        std::cout << buffer << std::endl;
+        std::cout << argv[0] << std::endl;
+    }
+
+    std::ofstream out("hdf5test.log");
+    std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+    std::cout.rdbuf(out.rdbuf()); //redirect std::cout to log.txt
+
+    auto const parset = getParameters(argc, argv);
+
+    // ----------------------
+    // set up contact network
+    // ----------------------
+    using BlocksFactory = StrikeSlipFactory<Grid, Vector>;
+    BlocksFactory blocksFactory(parset);
+
+    using ContactNetwork = typename BlocksFactory::ContactNetwork;
+    blocksFactory.build();
+
+    auto& contactNetwork = blocksFactory.contactNetwork();
+
+    const size_t bodyCount = contactNetwork.nBodies();
+
+   /* for (size_t i=0; i<contactNetwork.nLevels(); i++) {
+        // printDofLocation(contactNetwork.body(i)->gridView());
+
+        //Vector def(contactNetwork.deformedGrids()[i]->size(dims));
+        //def = 1;
+        //deformedGridComplex.setDeformation(def, i);
+
+        const auto& level = *contactNetwork.level(i);
+
+        for (size_t j=0; j<level.nBodies(); j++) {
+            writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
+        }
+    }*/
+
+    for (size_t i=0; i<bodyCount; i++) {
+        writeToVTK(contactNetwork.body(i)->gridView(), "", "initial_body_" + std::to_string(i) + "_leaf");
+    }
+
+    // ----------------------------
+    // assemble contactNetwork
+    // ----------------------------
+    contactNetwork.assemble();
+
+    //printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
+
+    std::vector<size_t> nVertices(bodyCount);
+    for (size_t i=0; i<bodyCount; i++) {
+        nVertices[i] = contactNetwork.body(i)->nVertices();
+    }
+
+    using MyProgramState = ProgramState<Vector, ScalarVector>;
+    MyProgramState programState(nVertices);
+    programState.setupInitialConditions(parset, contactNetwork);
+
+    auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+    for (size_t i=0; i<bodyCount; i++) {
+      contactNetwork.body(i)->setDeformation(programState.u[i]);
+    }
+    nBodyAssembler.assembleTransferOperator();
+    nBodyAssembler.assembleObstacle();
+
+    // ------------------------
+    // assemble global friction
+    // ------------------------
+    contactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
+
+    auto& globalFriction = contactNetwork.globalFriction();
+    globalFriction.updateAlpha(programState.alpha);
+
+    // -----------------
+    // init input/output
+    // -----------------
+    using Assembler = MyAssembler<DefLeafGridView, dims>;
+    using MyVertexBasis = typename Assembler::VertexBasis;
+    using MyCellBasis = typename Assembler::CellBasis;
+    std::vector<Vector> vertexCoordinates(bodyCount);
+    std::vector<const MyVertexBasis* > vertexBases(bodyCount);
+    std::vector<const MyCellBasis* > cellBases(bodyCount);
+
+    auto& wPatches = blocksFactory.weakPatches();
+    std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches(bodyCount);
+
+    for (size_t i=0; i<bodyCount; i++) {
+      const auto& body = contactNetwork.body(i);
+      vertexBases[i] = &(body->assembler()->vertexBasis);
+      cellBases[i] = &(body->assembler()->cellBasis);
+
+      weakPatches[i].resize(1);
+      weakPatches[i][0] = wPatches[i].get();
+
+      auto& vertexCoords = vertexCoordinates[i];
+      vertexCoords.resize(nVertices[i]);
+
+      auto hostLeafView = body->grid()->hostGrid().leafGridView();
+      Dune::MultipleCodimMultipleGeomTypeMapper<
+          LeafGridView, Dune::MCMGVertexLayout> const vertexMapper(hostLeafView, Dune::mcmgVertexLayout());
+      for (auto &&v : vertices(hostLeafView))
+        vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
+    }
+
+    std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches;
+    contactNetwork.frictionPatches(frictionPatches);
+
+    auto const writeData = parset.get<bool>("io.data.write");
+    auto dataFile = writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
+
+    auto dataWriter =
+        writeData ? std::make_unique<
+                        HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
+                        *dataFile, vertexCoordinates, vertexBases,
+                        frictionPatches, weakPatches)
+                  : nullptr;
+
+    IterationRegister iterationCount;
+
+    auto const report = [&](bool initial = false) {
+      if (writeData) {
+        dataWriter->reportSolution(programState, contactNetwork, globalFriction);
+        if (!initial)
+          dataWriter->reportIterations(programState, iterationCount);
+        dataFile->flush();
+      }
+    };
+    report(true);
+
+    std::cout.rdbuf(coutbuf); //reset to standard output again
+
+  } catch (Dune::Exception &e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+  } catch (std::exception &e) {
+    std::cerr << "Standard exception: " << e.what() << std::endl;
+  }
+}

dune/tectonic/tests/hdf5test/hdf5test.cfg

+# -*- mode:conf -*-
+gravity         = 9.81     # [m/s^2]
+
+[body0]
+length          = 0.5      # [m]
+height          = 0.5     # [m]
+depth           = 0.12     # [m]
+bulkModulus     = 1.5e5    # [Pa] #2190
+poissonRatio    = 0.11     # [1]  #0.11
+[body0.elastic]
+density         = 1300      # [kg/m^3] #750
+shearViscosity  = 0     # [Pas]
+bulkViscosity   = 0     # [Pas]
+[body0.viscoelastic]
+density         = 1300     # [kg/m^3]
+shearViscosity  = 0     # [Pas]
+bulkViscosity   = 0     # [Pas]
+
+[body1]
+length          = 0.5     # [m]
+height          = 0.5     # [m]
+depth           = 0.12     # [m]
+bulkModulus     = 1.5e5    # [Pa]
+poissonRatio    = 0.11     # [1]
+[body1.elastic]
+density         = 1300      # [kg/m^3]
+shearViscosity  = 0     # [Pas]
+bulkViscosity   = 0     # [Pas]
+[body1.viscoelastic]
+density         = 1300     # [kg/m^3]
+shearViscosity  = 0     # [Pas]
+bulkViscosity   = 0     # [Pas]
+
+[boundary.friction]
+C               = 6       # [Pa]
+mu0             = 0.48      # [ ]
+V0              = 1e-3     # [m/s]
+L               = 1e-6  # [m]
+initialAlpha    = 0        # [ ]
+stateModel      = AgeingLaw
+frictionModel   = Truncated #Regularised
+[boundary.friction.weakening]
+a               = 0.054    # [ ]
+b               = 0.074    # [ ]
+[boundary.friction.strengthening]
+a               = 0.054    # [ ]
+b               = 0.074    # [ ]
+
+[boundary.neumann]
+sigmaN          = 200.0      # [Pa]
+
+[boundary.dirichlet]
+finalVelocity   = 1e-5     # [m/s]
+
+[io]
+data.write      = true
+
+[problem]
+finalTime       = 1000     # [s] #1000
+bodyCount       = 2
+
+[initialTime]
+timeStep = 0
+relativeTime = 0.0
+relativeTau = 2e-2 # 1e-6
+
+[timeSteps]
+scheme = newmark
+timeSteps = 1000
+
+[u0.solver]
+maximumIterations = 100
+verbosity         = full
+
+[a0.solver]
+maximumIterations = 100
+verbosity         = full
+
+[v.solver]
+maximumIterations = 100
+verbosity         = quiet
+
+[v.fpi]
+maximumIterations = 10000
+lambda            = 0.5
+
+[solver.tnnmg.preconditioner]
+mode         = additive
+patchDepth   = 1
+maximumIterations = 2
+verbosity         = quiet
+[solver.tnnmg.preconditioner.patchsolver]
+maximumIterations = 100
+verbosity         = quiet
+[solver.tnnmg.preconditioner.basesolver]
+maximumIterations = 10000
+verbosity         = quiet
+
+[solver.tnnmg.main]
+pre   = 1
+multi = 5 # number of multigrid steps
+post  = 0
+
+
+

dune/tectonic/tests/nodalweights.cc

+#include "utils/almostequal.hh"
+
+#include "nodalweights.hh"
+
+template <class Basis0, class Basis1>
+NodalWeights<Basis0, Basis1>::NodalWeights(const Basis0& basis0, const Basis1& basis1)
+    : basis0_(basis0),
+      basis1_(basis1)
+    {}
+
+template <class Basis0, class Basis1>
+template <class Basis, class Element, class GlobalCoords>
+auto NodalWeights<Basis0, Basis1>::basisDof(const Basis& basis, const Element& elem, const GlobalCoords& vertex) const {
+    int dof = -1;
+
+    const typename Basis::LocalFiniteElement& lFE = basis.getLocalFiniteElement(elem);
+    const auto& geometry = elem.geometry();
+    const auto& localVertex = geometry.local(vertex);
+
+    const size_t localBasisSize = lFE.localBasis().size();
+    std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
+    lFE.localBasis().evaluateFunction(localVertex, localBasisRep);
+
+    for(size_t i=0; i<localBasisSize; i++) {
+        if (almost_equal(localBasisRep[i][0], 1.0, 2)) {
+            dof = basis.index(elem, i);
+            break;
+        }
+    }
+
+    return dof;
+}
+
+template <class Basis0, class Basis1>
+template <class GridGlue, class ScalarVector>
+void NodalWeights<Basis0, Basis1>::assemble(const GridGlue& glue, ScalarVector& weights0, ScalarVector& weights1, bool initializeVector) const {
+    using ctype = typename ScalarVector::field_type;
+
+    if (initializeVector==true) {
+        weights0.resize(basis0_.size());
+        weights1.resize(basis1_.size());
+    }
+
+    // loop over all intersections
+    for (const auto& rIs : intersections(glue)) {
+        const auto& inside = rIs.inside();
+        const auto& outside = rIs.outside();
+
+        /*if (!nmBoundary.contains(inside, rIs.indexInInside())) {
+            std::cout << "it happened" << std::endl;
+            continue;
+        }*/
+
+        const auto& geometry = rIs.geometry();
+        const ctype val = 1.0/(geometry.mydimension + 1)*geometry.volume();
+
+        std::cout << geometry.mydimension  << " " << geometry.volume() << " " << val << std::endl;
+
+        for (int i=0; i<geometry.corners(); i++) {
+            const auto& vertex = geometry.corner(i);
+
+            const int inIdx = basisDof(basis0_, inside, vertex);
+            if (inIdx >= 0) {
+                weights0[inIdx] += val;
+            }
+
+            const int outIdx = basisDof(basis1_, outside, vertex);
+            if (outIdx >= 0) {
+                weights1[outIdx] += val;
+            }
+
+            std::cout << inIdx << " " << outIdx << std::endl;
+        }
+    }
+}

dune/tectonic/tests/nodalweights.hh

+#ifndef SRC_NODALWEIGHTS_HH
+#define SRC_NODALWEIGHTS_HH
+
+
+/**
+ * let merged contact boundary Gamma be given by GridGlue object;
+ * computes
+ *                  int_Gamma lambda_i dx,
+ * where lambda_i is nodal basis of merged contact boundary whose
+ * dofs are given by nonmortar and mortar dofs;
+ *
+ *
+ * NOTE: works only for P1 nodal bases
+ **/
+
+template <class Basis0, class Basis1>
+class NodalWeights {
+private:
+    enum {dim = Basis0::GridView::Grid::dimension};
+
+    template <class Basis, class Element, class GlobalCoords>
+    auto basisDof(const Basis& basis, const Element& elem, const GlobalCoords& vertex) const;
+
+public:
+    NodalWeights(const Basis0& basis0, const Basis1& basis1);
+
+    template <class GridGlue, class ScalarVector>
+    void assemble(const GridGlue& glue, ScalarVector& weights0, ScalarVector& weights1, bool initializeVector = true) const;
+
+private:
+    const Basis0& basis0_;
+    const Basis1& basis1_;
+};
+
+#include "nodalweights.cc"
+
+#endif

dune/tectonic/tests/nodalweightstest.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <iostream>
+#include <exception>
+#include <fstream>
+
+#include <dune/grid/uggrid.hh>
+
+#include <dune/fufem/assemblers/assembler.hh>
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/functions/constantfunction.hh>
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/assemblers/dualmortarcoupling.hh>
+
+#include "../utils/debugutils.hh"
+#include "common.hh"
+
+#include "../nodalweights.hh"
+
+const int dim = 2;
+const int n = 5;
+const bool simplexGrid = true;
+
+const std::string path = "";
+const std::string outputFile = "nodalweightstest.log";
+
+#if HAVE_UG
+    using GridType = typename Dune::UGGrid<dim>;
+#else
+#error No UG!
+#endif
+    using LevelView = typename GridType::LevelGridView;
+    using LeafView = typename GridType::LeafGridView;
+    using LevelBoundaryPatch = BoundaryPatch<LevelView>;
+    using LeafBoundaryPatch = BoundaryPatch<LeafView>;
+
+    using c_type = double;
+    using GlobalCoords = Dune::FieldVector<c_type, dim>;
+    using BlockVector = Dune::BlockVector<GlobalCoords>;
+
+    using ScalarVector = Dune::FieldVector<c_type, 1>;
+    using ScalarBlockVector = Dune::BlockVector<ScalarVector>;
+
+    using CouplingPair = Dune::Contact::CouplingPair<GridType, GridType, c_type>;
+    using CouplingType = Dune::Contact::DualMortarCoupling<c_type, GridType>;
+
+int main(int argc, char *argv[]) { try {
+    Dune::MPIHelper::instance(argc, argv);
+
+    std::ofstream out(path + outputFile);
+    std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+    std::cout.rdbuf(out.rdbuf()); //redirect std::cout to outputFile
+
+    std::cout << "------------------------" << std::endl;
+    std::cout << "-- NodalWeights Test: --" << std::endl;
+    std::cout << "------------------------" << std::endl << std::endl;
+
+    // building grids
+    std::vector<std::shared_ptr<GridType>> grids(2);
+
+    GlobalCoords lowerLeft0({0, 0});
+    GlobalCoords upperRight0({2, 1});
+    buildGrid(lowerLeft0, upperRight0, n, grids[0]);
+
+    GlobalCoords lowerLeft1({0, 1});
+    GlobalCoords upperRight1({2, 2});
+    buildGrid(lowerLeft1, upperRight1, n, grids[1]);
+
+    // writing grids
+    for (size_t i=0; i<grids.size(); i++) {
+        const auto& levelView = grids[i]->levelGridView(0);
+        writeToVTK(levelView, path, "body_" + std::to_string(i) + "_level0");
+    }
+
+    // compute coupling boundaries
+    LevelView gridView0 = grids[0]->levelGridView(0);
+    LevelView gridView1 = grids[1]->levelGridView(0);
+    LevelBoundaryPatch upper(gridView0);
+    LevelBoundaryPatch lower(gridView1);
+
+    lower.insertFacesByProperty([&](typename LevelView::Intersection const &in) {
+        return xyBetween(lowerLeft1, {upperRight1[0], lowerLeft1[1]}, in.geometry().center());
+    });
+
+    upper.insertFacesByProperty([&](typename LevelView::Intersection const &in) {
+        return xyBetween({lowerLeft0[0], upperRight0[1]}, upperRight0, in.geometry().center());
+    });
+
+    // set contact coupling
+    Dune::Contact::NormalProjection<LeafBoundaryPatch> contactProjection;
+
+    CouplingPair coupling;
+    coupling.set(0, 1, upper, lower, 0.1, CouplingPair::CouplingType::STICK_SLIP, contactProjection, nullptr);
+
+    double coveredArea_ = 0.8;
+    CouplingType contactCoupling;
+    contactCoupling.setGrids(*grids[0], *grids[1]);
+    contactCoupling.setupContactPatch(*coupling.patch0(),*coupling.patch1());
+    contactCoupling.gridGlueBackend_ = coupling.backend();
+    contactCoupling.setCoveredArea(coveredArea_);
+    contactCoupling.setup();
+
+    using Basis = P1NodalBasis<LevelView, c_type>;
+    Basis basis0(grids[0]->levelGridView(0));
+    Basis basis1(grids[1]->levelGridView(0));
+
+    printBasisDofLocation(basis0);
+    printBasisDofLocation(basis1);
+
+    ScalarBlockVector nWeights0, nWeights1;
+    NodalWeights<Basis, Basis> nodalWeights(basis0, basis1);
+    nodalWeights.assemble(*contactCoupling.getGlue(), nWeights0, nWeights1, true);
+
+    ScalarBlockVector weights0;
+    {
+        Assembler<Basis, Basis> assembler(basis0, basis0);
+        NeumannBoundaryAssembler<GridType, typename ScalarBlockVector::block_type> boundaryAssembler(
+                     std::make_shared<ConstantFunction<
+                     GlobalCoords, typename ScalarBlockVector::block_type>>(1));
+             assembler.assembleBoundaryFunctional(boundaryAssembler, weights0, upper);
+    }
+
+    ScalarBlockVector weights1;
+    {
+        Assembler<Basis, Basis> assembler(basis1, basis1);
+        NeumannBoundaryAssembler<GridType, typename ScalarBlockVector::block_type> boundaryAssembler(
+                     std::make_shared<ConstantFunction<
+                     GlobalCoords, typename ScalarBlockVector::block_type>>(1));
+             assembler.assembleBoundaryFunctional(boundaryAssembler, weights1, lower);
+    }
+
+    print(weights0, "assembled weights0: ");
+    print(nWeights0, "nWeights0: ");
+
+    print(weights1, "assembled weights1: ");
+    print(nWeights1, "nWeights1: ");
+
+    bool sizePassed = true;
+    if (weights0.size() != nWeights0.size() | weights1.size() != nWeights1.size()) {
+        sizePassed = false;
+    }
+
+    bool entriesPassed0 = true;
+    for (size_t i=0; i<weights0.size(); i++) {
+        entriesPassed0 = entriesPassed0 & isClose(weights0[i], nWeights0[i]);
+    }
+
+    bool entriesPassed1 = true;
+    for (size_t i=0; i<weights1.size(); i++) {
+        entriesPassed1 = entriesPassed1 & isClose(weights1[i], nWeights1[i]);
+    }
+
+    bool passed = sizePassed & entriesPassed0 & entriesPassed1;
+
+    std::cout << "Overall, the test " << (passed ? "was successful!" : "failed!") << std::endl;
+
+    std::cout.rdbuf(coutbuf); //reset to standard output again
+    return passed ? 0 : 1;
+
+} catch (Dune::Exception &e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+} catch (std::exception &e) {
+    std::cerr << "Standard exception: " << e.what() << std::endl;
+} // end try
+} // end main

dune/tectonic/tests/nonoverlappingcouplingtest.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include "config.h"
+
+#include <array>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/stdstreams.hh>
+#include <iostream>
+
+#include <dune/common/fvector.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/geometrygrid.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/grid-glue/extractors/codim1extractor.hh>
+
+#include <dune/grid-glue/merging/contactmerge.hh>
+#include <dune/grid-glue/gridglue.hh>
+
+#include <dune/grid-glue/test/couplingtest.hh>
+#include <dune/grid-glue/test/communicationtest.hh>
+
+using namespace Dune;
+using namespace Dune::GridGlue;
+
+template<typename GridView>
+typename Dune::GridGlue::Codim1Extractor<GridView>::Predicate
+makeVerticalFacePredicate(double sliceCoord)
+{
+  using Element = typename GridView::Traits::template Codim<0>::Entity;
+  auto predicate = [sliceCoord](const Element& element, unsigned int face) -> bool {
+    const int dim = GridView::dimension;
+    const auto& refElement = Dune::ReferenceElements<double, dim>::general(element.type());
+
+    int numVertices = refElement.size(face, 1, dim);
+
+    for (int i=0; i<numVertices; i++)
+      if ( std::abs(element.geometry().corner(refElement.subEntity(face,1,i,dim))[0] - sliceCoord) > 1e-6 )
+        return false;
+
+    return true;
+  };
+  return predicate;
+}
+
+/** \brief trafo used for yaspgrids */
+template<int dim, typename ctype>
+class ShiftTrafo
+  : public AnalyticalCoordFunction< ctype, dim, dim, ShiftTrafo<dim,ctype> >
+{
+  double shift;
+public:
+  ShiftTrafo(double x) : shift(x) {};
+
+  //! evaluate method for global mapping
+  void evaluate ( const Dune::FieldVector<ctype, dim> &x, Dune::FieldVector<ctype, dim> &y ) const
+  {
+    y = x;
+    y[0] += shift;
+  }
+};
+
+template <int dim>
+void testMatchingCubeGrids()
+{
+
+  // ///////////////////////////////////////
+  //   Make two cube grids
+  // ///////////////////////////////////////
+
+  using GridType = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<double, dim> >;
+
+  std::array<int, dim> elements;
+  elements.fill(8);
+  FieldVector<double,dim> lower(0);
+  FieldVector<double,dim> upper(1);
+
+  GridType cubeGrid0(lower, upper, elements);
+
+  lower[0] += 1;
+  upper[0] += 1;
+
+  GridType cubeGrid1(lower, upper, elements);
+
+
+  // ////////////////////////////////////////
+  //   Set up coupling at their interface
+  // ////////////////////////////////////////
+
+  typedef typename GridType::LevelGridView DomGridView;
+  typedef typename GridType::LevelGridView TarGridView;
+
+  typedef Codim1Extractor<DomGridView> DomExtractor;
+  typedef Codim1Extractor<TarGridView> TarExtractor;
+
+  const typename DomExtractor::Predicate domdesc = makeVerticalFacePredicate<DomGridView>(1);
+  const typename TarExtractor::Predicate tardesc = makeVerticalFacePredicate<TarGridView>(1);
+
+  auto domEx = std::make_shared<DomExtractor>(cubeGrid0.levelGridView(0), domdesc);
+  auto tarEx = std::make_shared<TarExtractor>(cubeGrid1.levelGridView(0), tardesc);
+
+  typedef Dune::GridGlue::GridGlue<DomExtractor,TarExtractor> GlueType;
+
+  // Testing with ContactMerge
+  typedef ContactMerge<dim,double> ContactMergeImpl;
+
+  auto contactMerger = std::make_shared<ContactMergeImpl>(0.01);
+  GlueType contactGlue(domEx, tarEx, contactMerger);
+
+  contactGlue.build();
+
+  std::cout << "Gluing successful, " << contactGlue.size() << " remote intersections found!" << std::endl;
+  assert(contactGlue.size() > 0);
+
+  // ///////////////////////////////////////////
+  //   Test the coupling
+  // ///////////////////////////////////////////
+
+  testCoupling(contactGlue);
+  testCommunication(contactGlue);
+}
+
+
+template <int dim>
+void testNonMatchingCubeGrids()
+{
+
+  // ///////////////////////////////////////
+  //   Make two cube grids
+  // ///////////////////////////////////////
+
+  using GridType = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<double, dim> >;
+
+  std::array<int, dim> elements;
+  elements.fill(8);
+  FieldVector<double,dim> lower(0);
+  FieldVector<double,dim> upper(1);
+
+  GridType cubeGrid0(lower, upper, elements);
+
+  elements.fill(10);
+  lower[0] += 1;
+  upper[0] += 1;
+
+  GridType cubeGrid1(lower, upper, elements);
+
+
+  // ////////////////////////////////////////
+  //   Set up coupling at their interface
+  // ////////////////////////////////////////
+
+  typedef typename GridType::LevelGridView DomGridView;
+  typedef typename GridType::LevelGridView TarGridView;
+
+  typedef Codim1Extractor<DomGridView> DomExtractor;
+  typedef Codim1Extractor<TarGridView> TarExtractor;
+
+  const typename DomExtractor::Predicate domdesc = makeVerticalFacePredicate<DomGridView>(1);
+  const typename TarExtractor::Predicate tardesc = makeVerticalFacePredicate<TarGridView>(1);
+
+  auto domEx = std::make_shared<DomExtractor>(cubeGrid0.levelGridView(0), domdesc);
+  auto tarEx = std::make_shared<TarExtractor>(cubeGrid1.levelGridView(0), tardesc);
+
+  typedef Dune::GridGlue::GridGlue<DomExtractor,TarExtractor> GlueType;
+
+  // Testing with ContactMerge
+  typedef ContactMerge<dim,double> ContactMergeImpl;
+
+  auto contactMerger = std::make_shared<ContactMergeImpl>(0.01);
+  GlueType contactGlue(domEx, tarEx, contactMerger);
+
+  contactGlue.build();
+
+  std::cout << "Gluing successful, " << contactGlue.size() << " remote intersections found!" << std::endl;
+  assert(contactGlue.size() > 0);
+
+  // ///////////////////////////////////////////
+  //   Test the coupling
+  // ///////////////////////////////////////////
+
+  testCoupling(contactGlue);
+  testCommunication(contactGlue);
+}
+
+
+template<int dim, bool par>
+class MeshGenerator
+{
+  bool tar;
+public:
+  MeshGenerator(bool b) : tar(b) {};
+
+  using GridType = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<double, dim> >;
+
+  std::shared_ptr<GridType> generate()
+  {
+    std::array<int, dim> elements;
+    elements.fill(8);
+    FieldVector<double,dim> lower(0);
+    FieldVector<double,dim> upper(1);
+
+    if (tar)
+    {
+      elements.fill(10);
+      lower[0] += 1;
+      upper[0] += 1;
+    }
+
+    return std::make_shared<GridType>(lower, upper, elements);
+  }
+};
+
+
+template<int dim>
+class MeshGenerator<dim, true>
+{
+  bool tar;
+public:
+  MeshGenerator(bool b) : tar(b) {};
+
+  typedef YaspGrid<dim> HostGridType;
+  typedef GeometryGrid<HostGridType, ShiftTrafo<dim,double> > GridType;
+
+  std::shared_ptr<GridType> generate()
+  {
+    std::array<int,dim> elements;
+    std::fill(elements.begin(), elements.end(), 8);
+    std::bitset<dim> periodic(0);
+    FieldVector<double,dim> size(1);
+    int overlap = 1;
+    double shift = 0.0;
+
+    if (tar)
+    {
+      std::fill(elements.begin(), elements.end(), 10);
+      shift = 1.0;
+    }
+
+    HostGridType * hostgridp = new HostGridType(
+      size, elements, periodic, overlap
+#if HAVE_MPI
+      , MPI_COMM_WORLD
+#endif
+      );
+    ShiftTrafo<dim,double> * trafop = new ShiftTrafo<dim,double>(shift);
+    return std::make_shared<GridType>(*hostgridp, *trafop);
+  }
+};
+
+
+template <int dim, class DomGen, class TarGen>
+void testParallelCubeGrids()
+{
+  // ///////////////////////////////////////
+  //   Make two cube grids
+  // ///////////////////////////////////////
+
+  typedef typename DomGen::GridType GridType0;
+  typedef typename TarGen::GridType GridType1;
+
+  DomGen domGen(0);
+  TarGen tarGen(1);
+
+  double slice = 1.0;
+
+  std::shared_ptr<GridType0> cubeGrid0 = domGen.generate();
+  std::shared_ptr<GridType1> cubeGrid1 = tarGen.generate();
+
+  // ////////////////////////////////////////
+  //   Set up Traits
+  // ////////////////////////////////////////
+
+  typedef typename GridType0::LevelGridView DomGridView;
+  typedef typename GridType1::LevelGridView TarGridView;
+
+  typedef Codim1Extractor<DomGridView> DomExtractor;
+  typedef Codim1Extractor<TarGridView> TarExtractor;
+
+  const typename DomExtractor::Predicate domdesc = makeVerticalFacePredicate<DomGridView>(slice);
+  const typename TarExtractor::Predicate tardesc = makeVerticalFacePredicate<TarGridView>(slice);
+
+  auto domEx = std::make_shared<DomExtractor>(cubeGrid0->levelGridView(0), domdesc);
+  auto tarEx = std::make_shared<TarExtractor>(cubeGrid1->levelGridView(0), tardesc);
+
+  // ////////////////////////////////////////
+  //   Set up coupling at their interface
+  // ////////////////////////////////////////
+
+  typedef Dune::GridGlue::GridGlue<DomExtractor,TarExtractor> GlueType;
+
+  // Testing with ContactMerge
+  typedef ContactMerge<dim,double> ContactMergeImpl;
+
+  auto contactMerger = std::make_shared<ContactMergeImpl>(0.01);
+  GlueType contactGlue(domEx, tarEx, contactMerger);
+
+  contactGlue.build();
+
+  std::cout << "Gluing successful, " << contactGlue.size() << " remote intersections found!" << std::endl;
+  assert(contactGlue.size() > 0);
+
+  // ///////////////////////////////////////////
+  //   Test the coupling
+  // ///////////////////////////////////////////
+
+  testCoupling(contactGlue);
+  testCommunication(contactGlue);
+}
+
+#if HAVE_MPI
+void eh( MPI_Comm *comm, int *err, ... )
+{
+  int len = 1024;
+  char error_txt[len];
+
+  MPI_Error_string(*err, error_txt, &len);
+  assert(len <= 1024);
+  DUNE_THROW(Dune::Exception, "MPI ERROR -- " << error_txt);
+}
+#endif // HAVE_MPI
+
+int main(int argc, char *argv[])
+{
+  Dune::MPIHelper::instance(argc, argv);
+  Dune::dinfo.attach(std::cout);
+
+#if HAVE_MPI
+  MPI_Errhandler errhandler;
+  MPI_Comm_create_errhandler(eh, &errhandler);
+  MPI_Comm_set_errhandler(MPI_COMM_WORLD, errhandler);
+#endif
+
+  // 2d Tests
+  typedef MeshGenerator<2,false>  Seq;
+  typedef MeshGenerator<2,true>   Par;
+
+  // Test two unit squares
+  std::cout << "==== 2D hybrid =============================================\n";
+  testMatchingCubeGrids<2>();
+  std::cout << "============================================================\n";
+  testNonMatchingCubeGrids<2>();
+  std::cout << "============================================================\n";
+  testParallelCubeGrids<2,Seq,Seq>();
+  std::cout << "============================================================\n";
+  testParallelCubeGrids<2,Par,Seq>();
+  std::cout << "============================================================\n";
+  testParallelCubeGrids<2,Seq,Par>();
+  std::cout << "============================================================\n";
+  testParallelCubeGrids<2,Par,Par>();
+  std::cout << "============================================================\n";
+
+  // 3d Tests
+#if ! HAVE_MPI
+  typedef MeshGenerator<3,false>  Seq3d;
+  typedef MeshGenerator<3,true>   Par3d;
+
+  // Test two unit cubes
+  std::cout << "==== 3D hybrid =============================================\n";
+  testMatchingCubeGrids<3>();
+  std::cout << "============================================================\n";
+  testNonMatchingCubeGrids<3>();
+  std::cout << "============================================================\n";
+  testParallelCubeGrids<3,Seq3d,Seq3d>();
+  std::cout << "============================================================\n";
+  testParallelCubeGrids<3,Par3d,Seq3d>();
+  std::cout << "============================================================\n";
+  testParallelCubeGrids<3,Seq3d,Par3d>();
+  std::cout << "============================================================\n";
+  testParallelCubeGrids<3,Par3d,Par3d>();
+  std::cout << "============================================================\n";
+#endif // HAVE_MPI
+
+  return 0;
+}

dune/tectonic/tests/solverfactorytest.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#define MY_DIM 2
+
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <exception>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/stdstreams.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/fufem/formatstring.hh>
+
+#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
+#include <dune/tnnmg/functionals/bcqfconstrainedlinearization.hh>
+#include <dune/tnnmg/localsolvers/scalarobstaclesolver.hh>
+#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+
+#include <dune/tectonic/geocoordinate.hh>
+
+#include "assemblers.hh"
+#include "gridselector.hh"
+#include "explicitgrid.hh"
+#include "explicitvectors.hh"
+
+#include "data-structures/enumparser.hh"
+#include "data-structures/enums.hh"
+#include "data-structures/contactnetwork.hh"
+#include "data-structures/matrices.hh"
+#include "data-structures/program_state.hh"
+
+#include "io/vtk.hh"
+
+#include "spatial-solving/tnnmg/linesearchsolver.hh"
+#include "spatial-solving/preconditioners/nbodycontacttransfer.hh"
+
+#include "factories/stackedblocksfactory.hh"
+
+#include "time-stepping/rate.hh"
+#include "time-stepping/state.hh"
+#include "time-stepping/updaters.hh"
+
+#include "utils/tobool.hh"
+#include "utils/debugutils.hh"
+
+
+const int dim = MY_DIM;
+
+Dune::ParameterTree parset;
+size_t bodyCount;
+
+std::vector<BodyState<Vector, ScalarVector>* > bodyStates;
+std::vector<Vector> u;
+std::vector<Vector> v;
+std::vector<Vector> a;
+std::vector<ScalarVector> alpha;
+std::vector<ScalarVector> weightedNormalStress;
+double relativeTime;
+double relativeTau;
+size_t timeStep = 0;
+size_t timeSteps = 0;
+
+const std::string path = "";
+const std::string outputFile = "solverfactorytest.log";
+
+std::vector<std::vector<double>> allReductionFactors;
+
+template<class IterationStepType, class NormType, class ReductionFactorContainer>
+Dune::Solvers::Criterion reductionFactorCriterion(
+      const IterationStepType& iterationStep,
+      const NormType& norm,
+      ReductionFactorContainer& reductionFactors)
+{
+  double normOfOldCorrection = 1;
+  auto lastIterate = std::make_shared<typename IterationStepType::Vector>(*iterationStep.getIterate());
+
+  return Dune::Solvers::Criterion(
+      [&, lastIterate, normOfOldCorrection] () mutable {
+        double normOfCorrection = norm.diff(*lastIterate, *iterationStep.getIterate());
+        double convRate = (normOfOldCorrection > 0) ? normOfCorrection / normOfOldCorrection : 0.0;
+
+        if (convRate>1.0)
+            std::cout << "Solver convergence rate of " << convRate << std::endl;
+
+        normOfOldCorrection = normOfCorrection;
+        *lastIterate = *iterationStep.getIterate();
+
+        reductionFactors.push_back(convRate);
+        return std::make_tuple(convRate < 0, Dune::formatString(" % '.5f", convRate));
+      },
+      " reductionFactor");
+}
+
+
+template<class IterationStepType, class Functional, class ReductionFactorContainer>
+Dune::Solvers::Criterion energyCriterion(
+      const IterationStepType& iterationStep,
+      const Functional& f,
+      ReductionFactorContainer& reductionFactors)
+{
+  double normOfOldCorrection = 1;
+  auto lastIterate = std::make_shared<typename IterationStepType::Vector>(*iterationStep.getIterate());
+
+  return Dune::Solvers::Criterion(
+      [&, lastIterate, normOfOldCorrection] () mutable {
+        double normOfCorrection = std::abs(f(*lastIterate) - f(*iterationStep.getIterate())); //norm.diff(*lastIterate, *iterationStep.getIterate());
+
+        double convRate = (normOfOldCorrection != 0.0) ? 1.0 - (normOfCorrection / normOfOldCorrection) : 0.0;
+
+        if (convRate>1.0)
+            std::cout << "Solver convergence rate of " << convRate << std::endl;
+
+        normOfOldCorrection = normOfCorrection;
+        *lastIterate = *iterationStep.getIterate();
+
+        reductionFactors.push_back(convRate);
+        return std::make_tuple(convRate < 0, Dune::formatString(" % '.5f", convRate));
+      },
+      " reductionFactor");
+}
+
+template <class ContactNetwork>
+void solveProblem(const ContactNetwork& contactNetwork,
+                  const Matrix& mat, const Vector& rhs, Vector& x,
+                  const BitVector& ignore, const Vector& lower, const Vector& upper, bool initial = false) {
+
+    using Solver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+    using Norm =  EnergyNorm<Matrix, Vector>;
+
+    //using LocalSolver = LocalBisectionSolver;
+    //using Linearization = Linearization<Functional, BitVector>;
+
+    /*print(ignore, "ignore:");
+    for (size_t i=0; i<x.size(); i++) {
+        std::cout << x[i] << std::endl;
+    }*/
+
+    // set up reference solver
+    Vector refX = x;
+    using ContactFunctional = Dune::TNNMG::BoxConstrainedQuadraticFunctional<Matrix&, Vector&, Vector&, Vector&, double>;
+    auto I = ContactFunctional(mat, rhs, lower, upper);
+
+    std::cout << "Energy start iterate: " << I(x) << std::endl;
+
+    using LocalSolver = Dune::TNNMG::ScalarObstacleSolver;
+    auto localSolver = Dune::TNNMG::gaussSeidelLocalSolver(LocalSolver());
+
+    using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<ContactFunctional, decltype(localSolver), BitVector>;
+    auto nonlinearSmoother = std::make_shared<NonlinearSmoother>(I, refX, localSolver);
+
+    using Linearization = Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<ContactFunctional, BitVector>;
+    using DefectProjection = Dune::TNNMG::ObstacleDefectProjection;
+    using Step = Dune::TNNMG::TNNMGStep<ContactFunctional, BitVector, Linearization, DefectProjection, LocalSolver>;
+
+    // set multigrid solver
+    auto smoother = TruncatedBlockGSStep<Matrix, Vector>();
+
+    using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+    using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+    TransferOperators transfer(contactNetwork.nLevels()-1);
+    for (size_t i=0; i<transfer.size(); i++) {
+        transfer[i] = std::make_shared<TransferOperator>();
+        transfer[i]->setup(contactNetwork, i, i+1);
+    }
+
+    // Remove any recompute filed so that initially the full transferoperator is assembled
+    for (size_t i=0; i<transfer.size(); i++)
+        std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+    auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+    linearMultigridStep->setMGType(1, 3, 3);
+    linearMultigridStep->setSmoother(&smoother);
+    linearMultigridStep->setTransferOperators(transfer);
+
+    int mu = parset.get<int>("solver.tnnmg.main.multi"); // #multigrid steps in Newton step
+    auto step = Step(I, refX, nonlinearSmoother, linearMultigridStep, mu, DefectProjection(), LocalSolver());
+
+    // compute reference solution with generic functional and solver
+    auto norm = Norm(mat);
+
+    auto refSolver = Solver(step, parset.get<size_t>("u0.solver.maximumIterations"),
+                            parset.get<double>("u0.solver.tolerance"), norm, Solver::FULL);
+
+    step.setIgnore(ignore);
+    step.setPreSmoothingSteps(parset.get<int>("solver.tnnmg.main.pre"));
+
+    refSolver.addCriterion(
+            [&](){
+            return Dune::formatString("   % 12.5e", I(refX));
+            },
+            "   energy      ");
+
+    double initialEnergy = I(refX);
+    refSolver.addCriterion(
+            [&](){
+            static double oldEnergy=initialEnergy;
+            double currentEnergy = I(refX);
+            double decrease = currentEnergy - oldEnergy;
+            oldEnergy = currentEnergy;
+            return Dune::formatString("   % 12.5e", decrease);
+            },
+            "   decrease    ");
+
+    refSolver.addCriterion(
+            [&](){
+            return Dune::formatString("   % 12.5e", step.lastDampingFactor());
+            },
+            "   damping     ");
+
+
+    refSolver.addCriterion(
+            [&](){
+            return Dune::formatString("   % 12d", step.linearization().truncated().count());
+            },
+            "   truncated   ");
+
+    if (timeStep>0 and initial) {
+        allReductionFactors[timeStep].resize(0);
+        refSolver.addCriterion(reductionFactorCriterion(step, norm, allReductionFactors[timeStep]));
+    }
+
+    refSolver.preprocess();
+    refSolver.solve();
+
+    //print(refX, "refX: ");
+
+    // set up solver factory solver
+    using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
+
+    const auto& preconditionerParset = parset.sub("solver.tnnmg.preconditioner");
+
+    Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
+    Preconditioner preconditioner(preconditionerParset, contactNetwork, activeLevels);
+    preconditioner.setPatchDepth(preconditionerParset.get<size_t>("patchDepth"));
+    preconditioner.build();
+
+    auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector> >();
+    cgStep->setPreconditioner(preconditioner);
+
+    // set up functional
+    auto& globalFriction = contactNetwork.globalFriction();
+
+  /*  std::vector<const Dune::BitSetVector<1>*> nodes;
+     contactNetwork.frictionNodes(nodes);
+     print(*nodes[0], "frictionNodes: ");
+     print(*nodes[1], "frictionNodes: ");
+
+     print(ignore, "ignore: ");*/
+
+    //using MyFunctional = Functional<Matrix&, Vector&, std::decay_t<decltype(globalFriction)>&, Vector&, Vector&, typename Matrix::field_type>;
+    //MyFunctional J(mat, rhs, globalFriction, lower, upper);
+    using MyFunctional = Functional<Matrix&, Vector&, ZeroNonlinearity, Vector&, Vector&, typename Matrix::field_type>;
+    MyFunctional J(mat, rhs, ZeroNonlinearity(), lower, upper);
+
+    //std::cout << "ref energy: " << J(refX) << std::endl;
+
+    // set up TNMMG solver
+    // dummy solver, uses direct solver for linear correction no matter what is set here
+    //Norm mgNorm(*linearMultigridStep);
+    //auto mgSolver = std::make_shared<Solver>(linearMultigridStep, parset.get<size_t>("solver.tnnmg.linear.maximumIterations"), parset.get<double>("solver.tnnmg.linear.tolerance"), mgNorm, Solver::QUIET);
+    Norm mgNorm(*cgStep);
+    auto mgSolver = std::make_shared<Solver>(cgStep, parset.get<int>("solver.tnnmg.main.multi"), parset.get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), mgNorm, Solver::QUIET);
+
+    using Factory = SolverFactory<MyFunctional, BitVector>;
+    Factory factory(parset.sub("solver.tnnmg"), J, ignore);
+    factory.build(mgSolver);
+
+   /* std::vector<BitVector> bodyDirichletNodes;
+    nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
+    for (size_t i=0; i<bodyDirichletNodes.size(); i++) {
+      print(bodyDirichletNodes[i], "bodyDirichletNodes_" + std::to_string(i) + ": ");
+    }*/
+
+   /* print(bilinearForm, "matrix: ");
+    print(totalX, "totalX: ");
+    print(totalRhs, "totalRhs: ");*/
+
+    auto tnnmgStep = factory.step();
+    factory.setProblem(x);
+
+    LoopSolver<Vector> solver(
+        tnnmgStep.get(), parset.get<size_t>("u0.solver.maximumIterations"),
+                parset.get<double>("u0.solver.tolerance"), &norm,
+        Solver::FULL); //, true, &refX); // absolute error
+
+    solver.addCriterion(
+            [&](){
+            return Dune::formatString("   % 12.5e", J(x));
+            },
+            "   energy      ");
+
+    initialEnergy = J(x);
+    solver.addCriterion(
+            [&](){
+            static double oldEnergy=initialEnergy;
+            double currentEnergy = J(x);
+            double decrease = currentEnergy - oldEnergy;
+            oldEnergy = currentEnergy;
+            return Dune::formatString("   % 12.5e", decrease);
+            },
+            "   decrease    ");
+
+    solver.addCriterion(
+            [&](){
+            return Dune::formatString("   % 12.5e", tnnmgStep->lastDampingFactor());
+            },
+            "   damping     ");
+
+    solver.addCriterion(
+            [&](){
+            return Dune::formatString("   % 12d", tnnmgStep->linearization().truncated().count());
+            },
+            "   truncated   ");
+
+
+    /*std::vector<double> factors;
+    solver.addCriterion(reductionFactorCriterion(*tnnmgStep, norm, factors));
+
+    solver.addCriterion(energyCriterion(*tnnmgStep, J, factors));*/
+
+    solver.preprocess();
+    solver.solve();
+
+    auto diff = x;
+    diff -= refX;
+    std::cout << "Solver error in energy norm: " << norm(diff) << std::endl;
+
+    std::cout << "Energy end iterate: " << J(x) << std::endl;
+}
+
+
+template <class ContactNetwork>
+void setupInitialConditions(const ContactNetwork& contactNetwork) {
+  using Matrix = typename ContactNetwork::Matrix;
+  const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+
+  std::cout << std::endl << "-- setupInitialConditions --" << std::endl;
+  std::cout << "----------------------------" << std::endl;
+
+  // Solving a linear problem with a multigrid solver
+  auto const solveLinearProblem = [&](
+      const BitVector& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
+      const std::vector<Vector>& _rhs, std::vector<Vector>& _x) {
+
+    std::vector<const Matrix*> matrices_ptr(_matrices.size());
+    for (size_t i=0; i<matrices_ptr.size(); i++) {
+          matrices_ptr[i] = _matrices[i].get();
+    }
+
+    // assemble full global contact problem
+    Matrix bilinearForm;
+
+    nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+
+    Vector totalRhs, oldTotalRhs;
+    nBodyAssembler.assembleRightHandSide(_rhs, totalRhs);
+    oldTotalRhs = totalRhs;
+
+    Vector totalX, oldTotalX;
+    nBodyAssembler.nodalToTransformed(_x, totalX);
+    oldTotalX = totalX;
+
+    // get lower and upper obstacles
+    const auto& totalObstacles = nBodyAssembler.totalObstacles_;
+    Vector lower(totalObstacles.size());
+    Vector upper(totalObstacles.size());
+
+    for (size_t j=0; j<totalObstacles.size(); ++j) {
+        const auto& totalObstaclesj = totalObstacles[j];
+        auto& lowerj = lower[j];
+        auto& upperj = upper[j];
+      for (size_t d=0; d<dim; ++d) {
+          lowerj[d] = totalObstaclesj[d][0];
+          upperj[d] = totalObstaclesj[d][1];
+      }
+    }
+
+    // print problem
+   /* print(bilinearForm, "bilinearForm");
+    print(totalRhs, "totalRhs");
+    print(_dirichletNodes, "ignore");
+    print(totalObstacles, "totalObstacles");
+    print(lower, "lower");
+    print(upper, "upper");*/
+
+    solveProblem(contactNetwork, bilinearForm, totalRhs, totalX, _dirichletNodes, lower, upper, true);
+
+    nBodyAssembler.postprocess(totalX, _x);
+  };
+
+  timeStep = parset.get<size_t>("initialTime.timeStep");
+  relativeTime = parset.get<double>("initialTime.relativeTime");
+  relativeTau = parset.get<double>("initialTime.relativeTau");
+
+  bodyStates.resize(bodyCount);
+  u.resize(bodyCount);
+  v.resize(bodyCount);
+  a.resize(bodyCount);
+  alpha.resize(bodyCount);
+  weightedNormalStress.resize(bodyCount);
+
+  for (size_t i=0; i<bodyCount; i++) {
+    size_t leafVertexCount =  contactNetwork.body(i)->nVertices();
+
+    u[i].resize(leafVertexCount),
+    v[i].resize(leafVertexCount),
+    a[i].resize(leafVertexCount),
+    alpha[i].resize(leafVertexCount),
+    weightedNormalStress[i].resize(leafVertexCount),
+
+    bodyStates[i] = new BodyState<Vector, ScalarVector>(&u[i], &v[i], &a[i], &alpha[i], &weightedNormalStress[i]);
+  }
+
+  std::vector<Vector> ell0(bodyCount);
+  for (size_t i=0; i<bodyCount; i++) {
+    // Initial velocity
+    u[i] = 0.0;
+    v[i] = 0.0;
+
+    ell0[i].resize(u[i].size());
+    ell0[i] = 0.0;
+
+    contactNetwork.body(i)->externalForce()(relativeTime, ell0[i]);
+  }
+
+  // Initial displacement: Start from a situation of minimal stress,
+  // which is automatically attained in the case [v = 0 = a].
+  // Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
+  BitVector dirichletNodes;
+  contactNetwork.totalNodes("dirichlet", dirichletNodes);
+  /*for (size_t i=0; i<dirichletNodes.size(); i++) {
+      bool val = false;
+      for (size_t d=0; d<dims; d++) {
+          val = val || dirichletNodes[i][d];
+      }
+
+      dirichletNodes[i] = val;
+      for (size_t d=0; d<dims; d++) {
+          dirichletNodes[i][d] = val;
+      }
+  }*/
+
+  std::cout << "solving linear problem for u..." << std::endl;
+
+  solveLinearProblem(dirichletNodes, contactNetwork.matrices().elasticity, ell0, u);
+
+  //print(u, "initial u:");
+
+  // Initial acceleration: Computed in agreement with Ma = ell0 - Au
+  // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
+  std::vector<Vector> accelerationRHS = ell0;
+  for (size_t i=0; i<bodyCount; i++) {
+    // Initial state
+    alpha[i] = parset.get<double>("boundary.friction.initialAlpha");
+
+    // Initial normal stress
+
+    const auto& body = contactNetwork.body(i);
+    /*std::vector<std::shared_ptr<typename ContactNetwork::LeafBody::BoundaryCondition>> frictionBoundaryConditions;
+    body->boundaryConditions("friction", frictionBoundaryConditions);
+    for (size_t j=0; j<frictionBoundaryConditions.size(); j++) {
+        ScalarVector frictionBoundaryStress(weightedNormalStress[i].size());
+
+        body->assembler()->assembleWeightedNormalStress(
+          *frictionBoundaryConditions[j]->boundaryPatch(), frictionBoundaryStress, body->data()->getYoungModulus(),
+          body->data()->getPoissonRatio(), u[i]);
+
+        weightedNormalStress[i] += frictionBoundaryStress;
+    }*/
+
+    Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, u[i]);
+  }
+
+  for (size_t i=0; i<contactNetwork.nCouplings(); i++) {
+    const auto& coupling = contactNetwork.coupling(i);
+    const auto& contactCoupling = contactNetwork.nBodyAssembler().getContactCouplings()[i];
+
+    const auto nonmortarIdx = coupling->gridIdx_[0];
+    const auto& body = contactNetwork.body(nonmortarIdx);
+
+    ScalarVector frictionBoundaryStress(weightedNormalStress[nonmortarIdx].size());
+
+    body->assembler()->assembleWeightedNormalStress(
+      contactCoupling->nonmortarBoundary(), frictionBoundaryStress, body->data()->getYoungModulus(),
+      body->data()->getPoissonRatio(), u[nonmortarIdx]);
+
+    weightedNormalStress[nonmortarIdx] += frictionBoundaryStress;
+  }
+
+  //print(weightedNormalStress, "weightedNormalStress: ");
+
+  std::cout << "solving linear problem for a..." << std::endl;
+
+  BitVector noNodes(dirichletNodes.size(), false);
+  solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, a);
+
+  //print(a, "initial a:");
+}
+
+template <class ContactNetwork>
+void relativeVelocities(const ContactNetwork& contactNetwork, const Vector& v, std::vector<Vector>& v_rel) {
+
+    const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+    const size_t nBodies = nBodyAssembler.nGrids();
+   // const auto& contactCouplings = nBodyAssembler.getContactCouplings();
+
+    std::vector<const Dune::BitSetVector<1>*> bodywiseNonmortarBoundaries;
+    contactNetwork.frictionNodes(bodywiseNonmortarBoundaries);
+
+    size_t globalIdx = 0;
+    v_rel.resize(nBodies);
+    for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
+        const auto& nonmortarBoundary = *bodywiseNonmortarBoundaries[bodyIdx];
+        auto& v_rel_ = v_rel[bodyIdx];
+
+        v_rel_.resize(nonmortarBoundary.size());
+
+        for (size_t i=0; i<v_rel_.size(); i++) {
+            if (toBool(nonmortarBoundary[i])) {
+                v_rel_[i] = v[globalIdx];
+            }
+            globalIdx++;
+        }
+    }
+
+}
+
+template <class Updaters, class ContactNetwork>
+void run(Updaters& updaters, ContactNetwork& contactNetwork,
+    const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
+    std::vector<Vector>& velocityIterates) {
+
+  const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+  const auto nBodies = nBodyAssembler.nGrids();
+
+  std::vector<const Matrix*> matrices_ptr(nBodies);
+  for (int i=0; i<nBodies; i++) {
+      matrices_ptr[i] = &velocityMatrices[i];
+  }
+
+  // assemble full global contact problem
+  Matrix bilinearForm;
+  nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+
+  Vector totalRhs;
+  nBodyAssembler.assembleRightHandSide(velocityRHSs, totalRhs);
+
+  // get lower and upper obstacles
+  const auto& totalObstacles = nBodyAssembler.totalObstacles_;
+  Vector lower(totalObstacles.size());
+  Vector upper(totalObstacles.size());
+
+  for (size_t j=0; j<totalObstacles.size(); ++j) {
+      const auto& totalObstaclesj = totalObstacles[j];
+      auto& lowerj = lower[j];
+      auto& upperj = upper[j];
+    for (size_t d=0; d<dim; ++d) {
+        lowerj[d] = totalObstaclesj[d][0];
+        upperj[d] = totalObstaclesj[d][1];
+    }
+  }
+
+  // compute velocity obstacles
+  /*Vector vLower, vUpper;
+  std::vector<Vector> u0, v0;
+  updaters.rate_->extractOldVelocity(v0);
+  updaters.rate_->extractOldDisplacement(u0);
+
+  Vector totalu0, totalv0;
+  nBodyAssembler.nodalToTransformed(u0, totalu0);
+  nBodyAssembler.nodalToTransformed(v0, totalv0);
+
+  updaters.rate_->velocityObstacles(totalu0, lower, totalv0, vLower);
+  updaters.rate_->velocityObstacles(totalu0, upper, totalv0, vUpper);*/
+
+  const auto& errorNorms = contactNetwork.stateEnergyNorms();
+
+  auto& globalFriction = contactNetwork.globalFriction();
+
+  BitVector totalDirichletNodes;
+  contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
+
+  size_t fixedPointMaxIterations = parset.get<size_t>("v.fpi.maximumIterations");
+  double fixedPointTolerance = parset.get<double>("v.fpi.tolerance");
+  double lambda = parset.get<double>("v.fpi.lambda");
+
+  size_t fixedPointIteration;
+  size_t multigridIterations = 0;
+
+  std::vector<ScalarVector> alpha(nBodies);
+  updaters.state_->extractAlpha(alpha);
+
+  Vector totalVelocityIterate;
+  nBodyAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
+
+  // project in onto admissible set
+  const size_t blocksize = Vector::block_type::dimension;
+  for (size_t i=0; i<totalVelocityIterate.size(); i++) {
+      for (size_t j=0; j<blocksize; j++) {
+          if (totalVelocityIterate[i][j] < lower[i][j]) {
+              totalVelocityIterate[i][j] = lower[i][j];
+          }
+
+          if (totalVelocityIterate[i][j] > upper[i][j]) {
+              totalVelocityIterate[i][j] = upper[i][j];
+          }
+      }
+  }
+
+  for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations;
+       ++fixedPointIteration) {
+
+    // contribution from nonlinearity
+    //print(alpha, "alpha: ");
+
+    globalFriction.updateAlpha(alpha);
+
+    solveProblem(contactNetwork, bilinearForm, totalRhs, totalVelocityIterate, totalDirichletNodes, lower, upper, false);
+
+    nBodyAssembler.postprocess(totalVelocityIterate, velocityIterates);
+
+    std::vector<Vector> v_m; //TODO : wrong, isnt used atm;
+    updaters.rate_->extractOldVelocity(v_m);
+
+    for (size_t i=0; i<v_m.size(); i++) {
+      v_m[i] *= 1.0 - lambda;
+      Dune::MatrixVector::addProduct(v_m[i], lambda, velocityIterates[i]);
+    }
+
+    // extract relative velocities in mortar basis
+    std::vector<Vector> v_rel;
+    relativeVelocities(contactNetwork, totalVelocityIterate, v_rel);
+
+    //print(v_rel, "v_rel");
+
+    std::cout << "- State problem set" << std::endl;
+
+    // solve a state problem
+    updaters.state_->solve(v_rel);
+
+    std::cout << "-- Solved" << std::endl;
+
+    std::vector<ScalarVector> newAlpha(nBodies);
+    updaters.state_->extractAlpha(newAlpha);
+
+    bool breakCriterion = true;
+    for (int i=0; i<nBodies; i++) {
+        if (alpha[i].size()==0 || newAlpha[i].size()==0)
+            continue;
+
+        //print(alpha[i], "alpha i:");
+        //print(newAlpha[i], "new alpha i:");
+        if (errorNorms[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance) {
+            breakCriterion = false;
+            std::cout << "fixedPoint error: " << errorNorms[i]->diff(alpha[i], newAlpha[i]) << std::endl;
+            break;
+        }
+    }
+
+    if (lambda < 1e-12 or breakCriterion) {
+      std::cout << "-FixedPointIteration finished! " << (lambda < 1e-12 ? "lambda" : "breakCriterion") << std::endl;
+      fixedPointIteration++;
+      break;
+    }
+    alpha = newAlpha;
+  }
+
+  std::cout << "-FixedPointIteration finished! " << std::endl;
+
+  if (fixedPointIteration == fixedPointMaxIterations)
+    DUNE_THROW(Dune::Exception, "FPI failed to converge");
+
+  updaters.rate_->postProcess(velocityIterates);
+}
+
+
+template <class Updaters, class ContactNetwork>
+void step(Updaters& updaters, ContactNetwork& contactNetwork)  {
+  updaters.state_->nextTimeStep();
+  updaters.rate_->nextTimeStep();
+
+  auto const newRelativeTime = relativeTime + relativeTau;
+  typename ContactNetwork::ExternalForces externalForces;
+  contactNetwork.externalForces(externalForces);
+  std::vector<Vector> ell(externalForces.size());
+  for (size_t i=0; i<externalForces.size(); i++) {
+    (*externalForces[i])(newRelativeTime, ell[i]);
+  }
+
+  std::vector<Matrix> velocityMatrix;
+  std::vector<Vector> velocityRHS;
+  std::vector<Vector> velocityIterate;
+
+  double finalTime = parset.get<double>("problem.finalTime");
+  auto const tau = relativeTau * finalTime;
+  updaters.state_->setup(tau);
+  updaters.rate_->setup(ell, tau, newRelativeTime, velocityRHS, velocityIterate, velocityMatrix);
+
+  run(updaters, contactNetwork, velocityMatrix, velocityRHS, velocityIterate);
+}
+
+template <class Updaters, class ContactNetwork>
+void advanceStep(Updaters& current, ContactNetwork& contactNetwork) {
+    step(current, contactNetwork);
+    relativeTime += relativeTau;
+}
+
+void getParameters(int argc, char *argv[]) {
+  Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem.cfg", parset);
+  Dune::ParameterTreeParser::readINITree(
+      Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem-%dD.cfg", dim), parset);
+  Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+}
+
+int main(int argc, char *argv[]) { try {
+    Dune::MPIHelper::instance(argc, argv);
+
+    std::ofstream out(path + outputFile);
+    std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+    std::cout.rdbuf(out.rdbuf()); //redirect std::cout to outputFile
+
+    std::cout << "-------------------------" << std::endl;
+    std::cout << "-- SolverFactory Test: --" << std::endl;
+    std::cout << "-------------------------" << std::endl << std::endl;
+
+    getParameters(argc, argv);
+
+    // ----------------------
+    // set up contact network
+    // ----------------------
+    StackedBlocksFactory<Grid, Vector> stackedBlocksFactory(parset);
+    using ContactNetwork = typename StackedBlocksFactory<Grid, Vector>::ContactNetwork;
+    stackedBlocksFactory.build();
+
+    ContactNetwork& contactNetwork = stackedBlocksFactory.contactNetwork();
+    bodyCount = contactNetwork.nBodies();
+
+    for (size_t i=0; i<contactNetwork.nLevels(); i++) {
+        const auto& level = *contactNetwork.level(i);
+
+        for (size_t j=0; j<level.nBodies(); j++) {
+            writeToVTK(level.body(j)->gridView(), "debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
+        }
+    }
+
+    for (size_t i=0; i<bodyCount; i++) {
+        writeToVTK(contactNetwork.body(i)->gridView(), "debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
+    }
+
+    // ----------------------------
+    // assemble contactNetwork
+    // ----------------------------
+    contactNetwork.assemble();
+
+    //printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
+
+    // -----------------
+    // init input/output
+    // -----------------
+    timeSteps = parset.get<size_t>("timeSteps.timeSteps");
+    allReductionFactors.resize(timeSteps+1);
+
+    setupInitialConditions(contactNetwork);
+
+    auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+    for (size_t i=0; i<bodyCount; i++) {
+      contactNetwork.body(i)->setDeformation(u[i]);
+    }
+    nBodyAssembler.assembleTransferOperator();
+    nBodyAssembler.assembleObstacle();
+
+    // ------------------------
+    // assemble global friction
+    // ------------------------
+    contactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), weightedNormalStress);
+
+    auto& globalFriction = contactNetwork.globalFriction();
+    globalFriction.updateAlpha(alpha);
+
+    using Assembler = MyAssembler<DefLeafGridView, dim>;
+    using field_type = Matrix::field_type;
+    using MyVertexBasis = typename Assembler::VertexBasis;
+    using MyCellBasis = typename Assembler::CellBasis;
+    std::vector<Vector> vertexCoordinates(bodyCount);
+    std::vector<const MyVertexBasis* > vertexBases(bodyCount);
+    std::vector<const MyCellBasis* > cellBases(bodyCount);
+
+    for (size_t i=0; i<bodyCount; i++) {
+      const auto& body = contactNetwork.body(i);
+      vertexBases[i] = &(body->assembler()->vertexBasis);
+      cellBases[i] = &(body->assembler()->cellBasis);
+
+      auto& vertexCoords = vertexCoordinates[i];
+      vertexCoords.resize(body->nVertices());
+
+      Dune::MultipleCodimMultipleGeomTypeMapper<
+          DefLeafGridView, Dune::MCMGVertexLayout> const vertexMapper(body->gridView(), Dune::mcmgVertexLayout());
+      for (auto &&v : vertices(body->gridView()))
+        vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
+    }
+
+    const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/storage/mi/podlesny/software/dune/dune-tectonic/body");
+
+    auto const report = [&](bool initial = false) {
+      if (parset.get<bool>("io.printProgress"))
+        std::cout << "timeStep = " << std::setw(6) << timeStep
+                  << ", time = " << std::setw(12) << relativeTime
+                  << ", tau = " << std::setw(12) << relativeTau
+                  << std::endl;
+
+      if (parset.get<bool>("io.vtk.write")) {
+        std::vector<ScalarVector> stress(bodyCount);
+
+        for (size_t i=0; i<bodyCount; i++) {
+          const auto& body = contactNetwork.body(i);
+          body->assembler()->assembleVonMisesStress(body->data()->getYoungModulus(),
+                                           body->data()->getPoissonRatio(),
+                                           u[i], stress[i]);
+
+        }
+
+        vtkWriter.write(timeStep, u, v, alpha, stress);
+      }
+    };
+    report(true);
+
+    // -------------------
+    // Set up TNNMG solver
+    // -------------------
+
+    /*BitVector totalDirichletNodes;
+    contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
+
+    print(totalDirichletNodes, "totalDirichletNodes:");*/
+
+    //using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, field_type>;
+    //using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
+    //using NonlinearFactory = SolverFactory<Functional, BitVector>;
+
+    using BoundaryFunctions = typename ContactNetwork::BoundaryFunctions;
+    using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
+    using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
+                               StateUpdater<ScalarVector, Vector>>;
+
+    BoundaryFunctions velocityDirichletFunctions;
+    contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
+
+    BoundaryNodes dirichletNodes;
+    contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
+
+    /*for (size_t i=0; i<dirichletNodes.size(); i++) {
+        for (size_t j=0; j<dirichletNodes[i].size(); j++) {
+        print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
+        }
+    }*/
+
+
+    /*for (size_t i=0; i<frictionNodes.size(); i++) {
+        print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
+    }*/
+
+    Updaters current(
+        initRateUpdater(
+            parset.get<Config::scheme>("timeSteps.scheme"),
+            velocityDirichletFunctions,
+            dirichletNodes,
+            contactNetwork.matrices(),
+            u,
+            v,
+            a),
+        initStateUpdater<ScalarVector, Vector>(
+            parset.get<Config::stateModel>("boundary.friction.stateModel"),
+            alpha,
+            nBodyAssembler.getContactCouplings(),
+            contactNetwork.couplings())
+            );
+
+
+
+    //const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();
+
+    for (timeStep=1; timeStep<=timeSteps; timeStep++) {
+
+      advanceStep(current, contactNetwork);
+
+      relativeTime += relativeTau;
+
+      current.rate_->extractDisplacement(u);
+      current.rate_->extractVelocity(v);
+      current.rate_->extractAcceleration(a);
+      current.state_->extractAlpha(alpha);
+
+      contactNetwork.setDeformation(u);
+
+      report();
+    }
+
+    // output reduction factors
+    size_t count = 0;
+    for (size_t i=0; i<allReductionFactors.size(); i++) {
+        count = std::max(count, allReductionFactors[i].size());
+    }
+    std::vector<double> avgReductionFactors(count);
+    for (size_t i=0; i<count; i++) {
+        avgReductionFactors[i] = 1;
+        size_t c = 0;
+
+        for (size_t j=0; j<allReductionFactors.size(); j++) {
+            if (!(i<allReductionFactors[j].size()))
+                continue;
+
+            avgReductionFactors[i] *= allReductionFactors[j][i];
+            c++;
+        }
+
+        avgReductionFactors[i] = std::pow(avgReductionFactors[i], 1.0/((double)c));
+    }
+
+    print(avgReductionFactors, "average reduction factors: ");
+
+    bool passed = true;
+
+    std::cout << "Overall, the test " << (passed ? "was successful!" : "failed!") << std::endl;
+
+    std::cout.rdbuf(coutbuf); //reset to standard output again
+    return passed ? 0 : 1;
+
+} catch (Dune::Exception &e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+} catch (std::exception &e) {
+    std::cerr << "Standard exception: " << e.what() << std::endl;
+} // end try
+} // end main