Compare revisions

dc9fbb1c · dc9fbb1c · dc9fbb1c · dc9fbb1c · dc9fbb1c · dc9fbb1c
--- a/dune/tectonic/factories/levelcontactnetworkfactory.hh
+++ b/dune/tectonic/factories/levelcontactnetworkfactory.hh
+#ifndef SRC_CONTACTNETWORKFACTORY_HH
+#define SRC_CONTACTNETWORKFACTORY_HH
+
+#include <dune/common/parametertree.hh>
+
+#include "../data-structures/levelcontactnetwork.hh"
+
+template <class HostGridType>
+class ContactNetworkFactory {
+public:
+    using LevelContactNetwork = LevelContactNetwork<GridType, dims>;
+
+protected:
+    using Body = typename LevelContactNetwork::Body;
+    using FrictionCouplingPair = typename LevelContactNetwork::FrictionCouplingPair;
+
+    const Dune::ParameterTree& parset_;
+    const size_t bodyCount_;
+    const size_t couplingCount_;
+
+    std::vector<std::shared_ptr<Body>> bodies_;
+    std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
+
+    LevelContactNetwork levelContactNetwork_;
+
+private:
+    virtual void setBodies() = 0;
+    virtual void setCouplings() = 0;
+    virtual void setBoundaryConditions() = 0;
+
+public:
+    LevelContactNetworkFactory(const Dune::ParameterTree& parset, size_t bodyCount, size_t couplingCount) :
+        parset_(parset),
+        bodyCount_(bodyCount),
+        couplingCount_(couplingCount),
+        bodies_(bodyCount_),
+        couplings_(couplingCount_),
+        levelContactNetwork_(bodyCount_, couplingCount_) {}
+
+    void build() {
+        setBodies();
+        levelContactNetwork_.setBodies(bodies_);
+
+        setCouplings();
+        levelContactNetwork_.setCouplings(couplings_);
+
+        setBoundaryConditions();
+    }
+
+    LevelContactNetwork& levelContactNetwork() {
+        return levelContactNetwork_;
+    }
+};
+#endif
--- a/dune/tectonic/factories/stackedblocksfactory.cc
+++ b/dune/tectonic/factories/stackedblocksfactory.cc
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include "../problem-data/bc.hh"
+#include "../problem-data/grid/cuboidgeometry.hh"
+#include "../problem-data/myglobalfrictiondata.hh"
+
+#include "../utils/diameter.hh"
+
+#include "stackedblocksfactory.hh"
+
+
+template <class HostGridType, class VectorType>
+void StackedBlocksFactory<HostGridType, VectorType>::setBodies() {
+
+    // set up cuboid geometries
+    const auto& bodyParset = this->parset_.sub("body");
+
+    double const length = bodyParset.template get<double>("length");
+    double const height = bodyParset.template get<double>("height");
+    double const heightScaling = bodyParset.template get<double>("heightScaling");
+    double const weakLength = 1.00;
+
+    std::vector<GlobalCoords> origins(this->bodyCount_);
+
+    const auto& subParset = this->parset_.sub("boundary.friction.weakening");
+
+    std::vector<double> heights(this->bodyCount_);
+
+#if MY_DIM == 3
+        double const depth = 0.60;
+
+        auto weakBound = [&] (const GlobalCoords& origin) {
+            GlobalCoords h = {origin[0] + weakLength, origin[1], origin[2]};
+            return h;
+        };
+
+        origins[0] = {0,0,0};
+        GlobalCoords lowerWeakOrigin = {0.2, 0, 0};
+        GlobalCoords upperWeakOrigin = {0.2, height, 0};
+
+        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], length, height, depth);
+        cuboidGeometries_[0]->addWeakeningPatch(subParset, upperWeakOrigin, weakBound(upperWeakOrigin));
+
+        for (size_t i=1; i<this->bodyCount_-1; i++) {
+            origins[i] = cuboidGeometries_[i-1]->upperLeft();
+            GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
+            GlobalCoords upperWeakOrigin_i = upperWeakOrigin + origins[i];
+
+            cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height, depth);
+            cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));
+            cuboidGeometries_[i]->addWeakeningPatch(subParset, lowerWeakOrigin_i, weakBound(lowerWeakOrigin_i));
+        }
+
+        const size_t idx = this->bodyCount_-1;
+        origins[idx] = cuboidGeometries_[idx-1]->upperLeft();
+        lowerWeakOrigin += origins[idx];
+
+        cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], length, height, depth);
+        cuboidGeometries_[idx]->addWeakeningPatch(subParset, lowerWeakOrigin, weakBound(lowerWeakOrigin));
+#elif MY_DIM == 2
+        auto weakBound = [&] (const GlobalCoords& origin) {
+            GlobalCoords h = {origin[0] + weakLength*length, origin[1]};
+            return h;
+        };
+
+        origins[0] = {0,0};
+        GlobalCoords lowerWeakOrigin = {0.0, 0};
+        GlobalCoords upperWeakOrigin = {0.0, height};
+        heights[0] = height;
+
+        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], length, height);
+        cuboidGeometries_[0]->addWeakeningPatch(subParset, upperWeakOrigin, weakBound(upperWeakOrigin));
+
+        for (size_t i=1; i<this->bodyCount_-1; i++) {
+            origins[i] = cuboidGeometries_[i-1]->upperLeft();
+
+            double power = (i<this->bodyCount_/2.0) ? i : (this->bodyCount_-i-1);
+            heights[i] = height*std::pow(heightScaling, power);
+
+            GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
+            GlobalCoords upperWeakOrigin_i = {upperWeakOrigin[0], heights[i]};
+            upperWeakOrigin_i += origins[i];
+
+            cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, heights[i]);
+
+            cuboidGeometries_[i]->addWeakeningPatch(subParset, lowerWeakOrigin_i, weakBound(lowerWeakOrigin_i));
+            cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));  
+        }
+
+        const size_t idx = this->bodyCount_-1;
+        origins[idx] = cuboidGeometries_[idx-1]->upperLeft();
+        lowerWeakOrigin += origins[idx];
+        heights[idx] = height;
+
+        cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], length, height);
+        cuboidGeometries_[idx]->addWeakeningPatch(subParset, lowerWeakOrigin, weakBound(lowerWeakOrigin));
+#else
+#error CuboidGeometry only supports 2D and 3D!"
+#endif
+
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& cuboidGeometry = *cuboidGeometries_[i];
+
+        // set up reference grid
+        std::array<unsigned int, 2> initialElements = {
+            { (unsigned int) length/heights[i], 1 } };
+        DefaultCuboidGridConstructor<HostGridType> gridConstructor(cuboidGeometry, initialElements);
+
+        // define weak patch and refine grid
+        const auto& weakeningRegions = cuboidGeometry.weakeningRegions();
+        for (size_t j=0; j<weakeningRegions.size(); j++) {
+            gridConstructor.refine(weakeningRegions[j], this->parset_.template get<double>("boundary.friction.smallestDiameter"), CuboidGeometry::lengthScale());
+        }
+
+        this->bodies_[i] = std::make_shared<typename Base::LeafBody>(bodyData_, gridConstructor.grid());
+
+        // determine minDiameter and maxDiameter
+        double minDiameter = std::numeric_limits<double>::infinity();
+        double maxDiameter = 0.0;
+        for (auto &&e : elements(gridConstructor.grid()->leafGridView())) {
+          auto const geometry = e.geometry();
+          auto const diam = diameter(geometry);
+          minDiameter = std::min(minDiameter, diam);
+          maxDiameter = std::max(maxDiameter, diam);
+        }
+        std::cout << "Grid" << i << " min diameter: " << minDiameter << std::endl;
+        std::cout << "Grid" << i << " max diameter: " << maxDiameter << std::endl;
+    }
+}
+
+
+template <class HostGridType, class VectorType>
+void StackedBlocksFactory<HostGridType, VectorType>::setCouplings() {
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& cuboidGeometry = *cuboidGeometries_[i];
+
+        leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), cuboidGeometry);
+        levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), cuboidGeometry);
+    }
+
+    auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
+    std::shared_ptr<typename Base::FrictionCouplingPair::BackEndType> backend = nullptr;
+
+
+    for (size_t i=0; i<this->couplings_.size(); i++) {
+      auto& coupling = this->couplings_[i];
+      coupling = std::make_shared<typename Base::FrictionCouplingPair>();
+
+      nonmortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i+1]->lower);
+      mortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i]->upper);
+      weakPatches_[i] = std::make_shared<WeakeningRegion>(cuboidGeometries_[i+1]->weakeningRegions()[0]);
+
+      coupling->set(i+1, i, nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
+      coupling->setWeakeningPatch(weakPatches_[i]);
+      coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<GlobalCoords>>(this->parset_.sub("boundary.friction"), *weakPatches_[i], CuboidGeometry::lengthScale()));
+    }
+}
+
+template <class HostGridType, class VectorType>
+void StackedBlocksFactory<HostGridType, VectorType>::setLevelBodies() {
+    const size_t nBodies = this->bodyCount_;
+
+    size_t maxLevel = 0;
+    std::vector<size_t> maxLevels(nBodies);
+    for (size_t i=0; i<nBodies; i++) {
+        const size_t bodyMaxLevel = this->bodies_[i]->grid()->maxLevel();
+        maxLevels[i] = bodyMaxLevel;
+        maxLevel = std::max(maxLevel, bodyMaxLevel);
+    }
+
+    for (int l=maxLevel; l>=0; l--) {
+        this->contactNetwork_.addLevel(maxLevels, l);
+
+        for (size_t i=0; i<nBodies; i++) {
+            maxLevels[i] = (maxLevels[i]>0) ? maxLevels[i]-1 : 0;
+        }
+    }
+}
+
+
+template <class HostGridType, class VectorType>
+void StackedBlocksFactory<HostGridType, VectorType>::setBoundaryConditions() {
+    using LeafBoundaryCondition = BoundaryCondition<LeafGridView, dim>;
+
+    using Function = Dune::VirtualFunction<double, double>;
+    std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
+    std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>(this->parset_.template get<double>("boundary.dirichlet.finalVelocity"));
+
+    const double lengthScale = CuboidGeometry::lengthScale();
+
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& body = this->contactNetwork_.body(i);
+        const auto& leafVertexCount = body->nVertices();
+
+        // Neumann boundary
+        std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->gridView()), neumannFunction, "neumann");
+        body->addBoundaryCondition(neumannBoundary);
+
+        // upper Dirichlet Boundary
+        // identify Dirichlet nodes on leaf level
+        if (i==this->bodyCount_-1) {
+            std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount);
+            for (int j=0; j<leafVertexCount; j++) {
+                if (leafFaces_[i]->upper.containsVertex(j))
+                    (*velocityDirichletNodes)[j][0] = true;
+            }
+
+            std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+
+            velocityDirichletBoundary->setBoundaryPatch(body->gridView(), velocityDirichletNodes);
+            velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
+            body->addBoundaryCondition(velocityDirichletBoundary);
+        } else {
+
+        }
+    }
+
+    // lower Dirichlet Boundary
+    const auto& firstBody = this->contactNetwork_.body(0);
+    const auto& firstLeafVertexCount = firstBody->nVertices();
+    std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(firstLeafVertexCount);
+    for (int j=0; j<firstLeafVertexCount; j++) {
+        if (leafFaces_[0]->lower.containsVertex(j)) {
+            for (size_t d=0; d<dim; d++) {
+              (*zeroDirichletNodes)[j][d] = true;
+            }
+        }
+
+        #if MY_DIM == 3 //TODO: wrong, needs revision
+        if (leafFaces_[0]->front.containsVertex(j) || leafFaces_[0]->back.containsVertex(j))
+            (*zeroDirichletNodes)[j][2] = true;
+        #endif
+    }
+
+    std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+    zeroDirichletBoundary->setBoundaryPatch(firstBody->gridView(), zeroDirichletNodes);
+
+    std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
+    zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
+    firstBody->addBoundaryCondition(zeroDirichletBoundary);
+}
+
+#include "stackedblocksfactory_tmpl.cc"
--- a/dune/tectonic/factories/stackedblocksfactory.hh
+++ b/dune/tectonic/factories/stackedblocksfactory.hh
+#ifndef SRC_STACKEDBLOCKSFACTORY_HH
+#define SRC_STACKEDBLOCKSFACTORY_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include "contactnetworkfactory.hh"
+
+#include "../problem-data/mybody.hh"
+#include "../problem-data/grid/cuboidgridconstructor.hh"
+#include "../problem-data/grid/cuboidgeometry.hh"
+
+template <class HostGridType, class VectorType> class StackedBlocksFactory : public ContactNetworkFactory<HostGridType, VectorType>{
+private:
+    using Base = ContactNetworkFactory<HostGridType, VectorType>;
+
+public:
+    using ContactNetwork = typename Base::ContactNetwork;
+
+private:
+    using GlobalCoords = typename ContactNetwork::LocalVector;
+
+    using LeafGridView = typename ContactNetwork::GridView;
+    using LevelGridView = typename ContactNetwork::GridType::LevelGridView;
+
+    using LevelBoundaryPatch = BoundaryPatch<LevelGridView>;
+    using LeafBoundaryPatch = BoundaryPatch<LeafGridView>;
+
+    using LeafFaces = MyFaces<LeafGridView>;
+    using LevelFaces = MyFaces<LevelGridView>;
+
+    using CuboidGeometry= CuboidGeometry<typename GlobalCoords::field_type>;
+    using WeakeningRegion = typename CuboidGeometry::WeakeningRegion;
+
+    static const int dim = ContactNetwork::dim;
+
+    const std::shared_ptr<MyBodyData<dim>> bodyData_;     // material properties of bodies
+
+    std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries_;
+
+    std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
+    std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
+
+    std::vector<std::shared_ptr<WeakeningRegion>> weakPatches_;
+
+    std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
+    std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
+
+public:
+    StackedBlocksFactory(const Dune::ParameterTree& parset) :
+        Base(parset, parset.get<size_t>("problem.bodyCount"), parset.get<size_t>("problem.bodyCount")-1),
+        bodyData_(std::make_shared<MyBodyData<dim>>(this->parset_.sub("body"), this->parset_.template get<double>("general.gravity"), zenith_())),
+        cuboidGeometries_(this->bodyCount_),
+        leafFaces_(this->bodyCount_),
+        levelFaces_(this->bodyCount_),
+        weakPatches_(this->bodyCount_),
+        nonmortarPatch_(this->couplingCount_),
+        mortarPatch_(this->couplingCount_)
+    {}
+
+    void setBodies();
+    void setLevelBodies();
+    void setCouplings();
+    void setLevelCouplings() {}
+    void setBoundaryConditions();
+
+    auto& weakPatches() {
+        return weakPatches_;
+    }
+
+private:
+    static constexpr Dune::FieldVector<double, MY_DIM> zenith_() {
+        #if MY_DIM == 2
+        return {0, 1};
+        #elif MY_DIM == 3
+        return {0, 1, 0};
+        #endif
+    }
+};
+#endif
+
+
--- a/src/hdf5/surface-writer_tmpl.cc
+++ b/src/hdf5/surface-writer_tmpl.cc
-#include "../explicitvectors.hh"
-#include "../explicitgrid.hh"
-
-#include "../program_state.hh"
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif

-using MyProgramState = ProgramState<Vector, ScalarVector>;
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"

-template class SurfaceWriter<MyProgramState, GridView>;
+template class StackedBlocksFactory<Grid, Vector>;
--- a/dune/tectonic/factories/strikeslipfactory.cc
+++ b/dune/tectonic/factories/strikeslipfactory.cc
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include "../problem-data/grid/cuboidgeometry.hh"
+#include "../problem-data/bc.hh"
+#include "../problem-data/myglobalfrictiondata.hh"
+
+#include "../utils/diameter.hh"
+
+#include "strikeslipfactory.hh"
+
+template <class HostGridType, class VectorTEMPLATE>
+void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setBodies() {
+    // set up cuboid geometries
+
+    std::array<double, 2> lengths = {this->parset_.template get<double>("body0.length"), this->parset_.template get<double>("body1.length")};
+    std::array<double, 2> heights = {this->parset_.template get<double>("body0.height"), this->parset_.template get<double>("body1.height")};
+
+    GlobalCoords origin(0);
+
+    const auto& frictionParset = this->parset_.sub("boundary.friction");
+
+#if MY_DIM == 3 // TODO: not implemented
+        //std::array<double, 2> depths = {this->parset_.template get<double>("body0.depth"), this->parset_.template get<double>("body1.depth")};
+
+        //triangleGeometries_[0] = std::make_shared<TriangleGeometry>(origins[0], lengths[0], heights[0], depths[0]);
+        //triangleGeometries_[0]->addWeakeningPatch(frictionParset, {origins[0][0], origins[0][1]+ heights[0], 0}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0], 0});
+
+        //triangleGeometries_[1] = std::make_shared<TriangleGeometry>(origins[1], lengths[1], heights[1], depths[1]);
+        //triangleGeometries_[1]->addWeakeningPatch(frictionParset, origins[1], {origins[1][0] + lengths[1], origins[1][1], 0});
+#elif MY_DIM == 2
+        triangleGeometries_[0] = std::make_shared<TriangleGeometry>(origin, lengths[0], heights[0]);
+        triangleGeometries_[0]->addWeakeningPatch(frictionParset, triangleGeometries_[0]->A(), triangleGeometries_[0]->C());
+
+        triangleGeometries_[1] = std::make_shared<TriangleGeometry>(triangleGeometries_[0]->C(), -lengths[1], -heights[1]);
+        triangleGeometries_[1]->addWeakeningPatch(frictionParset, triangleGeometries_[1]->A(), triangleGeometries_[1]->C());
+#else
+#error CuboidGeometry only supports 2D and 3D!"
+#endif
+
+    // set up reference grids
+    gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(triangleGeometries_);
+    auto& grids = gridConstructor_->getGrids();
+
+    std::array<double, 2> smallestDiameter = {this->parset_.template get<double>("body0.smallestDiameter"), this->parset_.template get<double>("body1.smallestDiameter")};
+
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& triangleGeometry = *triangleGeometries_[i];
+
+        // define weak patch and refine grid
+        const auto& weakeningRegions = triangleGeometry.weakeningRegions();
+        for (size_t j=0; j<weakeningRegions.size(); j++) {
+            refine(*grids[i], weakeningRegions[j], smallestDiameter[i], TriangleGeometry::lengthScale());
+        }
+
+        // determine minDiameter and maxDiameter
+        double minDiameter = std::numeric_limits<double>::infinity();
+        double maxDiameter = 0.0;
+        for (auto &&e : elements(grids[i]->leafGridView())) {
+          auto const geometry = e.geometry();
+          auto const diam = diameter(geometry);
+          minDiameter = std::min(minDiameter, diam);
+          maxDiameter = std::max(maxDiameter, diam);
+        }
+        std::cout << "Grid" << i << " min diameter: " << minDiameter << std::endl;
+        std::cout << "Grid" << i << " max diameter: " << maxDiameter << std::endl;
+    }
+
+#if MY_DIM == 2
+    bodyData_[0] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body0"), 0.0, zenith_());
+    this->bodies_[0] = std::make_shared<typename Base::LeafBody>(bodyData_[0], grids[0]);
+
+    bodyData_[1] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body1"), 0.0, zenith_());
+    this->bodies_[1] = std::make_shared<typename Base::LeafBody>(bodyData_[1], grids[1]);
+#else
+    bodyData_[0] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body0"), this->parset_.template get<double>("general.gravity"), zenith_());
+    this->bodies_[0] = std::make_shared<typename Base::LeafBody>(bodyData_[0], grids[0]);
+
+    bodyData_[1] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body1"), this->parset_.template get<double>("general.gravity"), zenith_());
+    this->bodies_[1] = std::make_shared<typename Base::LeafBody>(bodyData_[1], grids[1]);
+#endif
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setLevelBodies() {
+    const size_t maxLevel = std::max({this->bodies_[0]->grid()->maxLevel(), this->bodies_[1]->grid()->maxLevel()});
+
+    for (size_t l=0; l<=maxLevel; l++) {
+        std::vector<size_t> bodyLevels(2, l);
+        this->contactNetwork_.addLevel(bodyLevels, l);
+    }
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setCouplings() {
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& triangleGeometry = *triangleGeometries_[i];
+        leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), triangleGeometry);
+        levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), triangleGeometry);
+    }
+
+    auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
+
+    nonmortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->b);
+    mortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->b);
+    weakPatches_[0] = std::make_shared<WeakeningRegion>(triangleGeometries_[1]->weakeningRegions()[0]);
+
+    auto& coupling = this->couplings_[0];
+    coupling = std::make_shared<typename Base::FrictionCouplingPair>();
+
+    coupling->set(1, 0, nonmortarPatch_[0], mortarPatch_[0], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, nullptr);
+    coupling->setWeakeningPatch(weakPatches_[0]);
+    coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<GlobalCoords>>(this->parset_.sub("boundary.friction"), *weakPatches_[0], TriangleGeometry::lengthScale()));
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
+    using LeafBoundaryCondition = BoundaryCondition<LeafGridView, dim>;
+
+    using Function = Dune::VirtualFunction<double, double>;
+    std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
+    std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>(-1.0*this->parset_.template get<double>("boundary.dirichlet.finalVelocity"), 0.25);
+
+    const double lengthScale = TriangleGeometry::lengthScale();
+
+    // body0
+    const auto& body0 = this->contactNetwork_.body(0);
+    const auto& leafVertexCount0 = body0->nVertices();
+    std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount0);
+    std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount0);
+
+    const auto& gridView0 = body0->gridView();
+    const auto& indexSet0 = gridView0.indexSet();
+
+    for (const auto& e : elements(gridView0)) {
+        for (const auto& is : intersections(gridView0, e)) {
+            if (!is.boundary())
+                continue;
+
+            auto isCenter = is.geometry().center();
+
+            // Dirichlet boundary (c)
+            if (leafFaces_[0]->c.contains(is)) {
+                isCenter -=  triangleGeometries_[0]->A();
+
+                if (isCenter.two_norm()/std::abs(triangleGeometries_[0]->length()) > 0.1) {
+                    const auto inside = is.inside();
+
+                    auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
+                    int n = refElement.size(is.indexInInside(),1,dim);
+
+                    for (int i=0; i<n; i++) {
+                        int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
+                        int globalIdx = indexSet0.template subIndex(inside,faceIdxi,dim);
+                        (*zeroDirichletNodes)[globalIdx][1] = true;
+                    }
+                }
+            }
+
+            // velocity Dirichlet boundary (a)
+            if (leafFaces_[0]->a.contains(is)) {
+                isCenter -=  triangleGeometries_[0]->C();
+
+                if (isCenter.two_norm()/std::abs(triangleGeometries_[0]->height()) > 0.1) {
+                    const auto inside = is.inside();
+
+                    auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
+                    int n = refElement.size(is.indexInInside(),1,dim);
+
+                    for (int i=0; i<n; i++) {
+                        int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
+                        int globalIdx = indexSet0.template subIndex(inside,faceIdxi,dim);
+                        (*velocityDirichletNodes)[globalIdx][0] = true;
+                    }
+                }
+            }
+        }
+    }
+
+    std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
+
+    std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+    zeroDirichletBoundary->setBoundaryPatch(body0->gridView(), zeroDirichletNodes);
+    zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
+    body0->addBoundaryCondition(zeroDirichletBoundary);
+
+    std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+    velocityDirichletBoundary->setBoundaryPatch(body0->gridView(), velocityDirichletNodes);
+    velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
+    body0->addBoundaryCondition(velocityDirichletBoundary);
+
+    // body1
+    const auto& body1 = this->contactNetwork_.body(1);
+    const auto& leafVertexCount1 = body1->nVertices();
+    std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes1 = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
+    std::shared_ptr<Dune::BitSetVector<dim>> loadNeumannNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
+
+    const auto& gridView1 = body1->gridView();
+    const auto& indexSet1 = gridView1.indexSet();
+
+    for (const auto& e : elements(gridView1)) {
+        for (const auto& is : intersections(gridView1, e)) {
+            if (!is.boundary())
+                continue;
+
+            auto isCenter = is.geometry().center();
+
+            // Neumann boundary, normal load (c)
+            if (leafFaces_[1]->c.contains(is)) {
+                isCenter -=  triangleGeometries_[1]->A();
+
+                if (isCenter.two_norm()/std::abs(triangleGeometries_[1]->length()) > 0.1) {
+                    const auto inside = is.inside();
+
+                    auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
+                    int n = refElement.size(is.indexInInside(),1,dim);
+
+                    for (int i=0; i<n; i++) {
+                        int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
+                        int globalIdx = indexSet1.template subIndex(inside,faceIdxi,dim);
+                        (*loadNeumannNodes)[globalIdx][1] = true;
+                    }
+                }
+            }
+
+            // Dirichlet boundary (a)
+            if (leafFaces_[1]->a.contains(is)) {
+                isCenter -=  triangleGeometries_[1]->C();
+
+                if (isCenter.two_norm()/std::abs(triangleGeometries_[1]->height()) > 0.1) {
+                    const auto inside = is.inside();
+
+                    auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
+                    int n = refElement.size(is.indexInInside(),1,dim);
+
+                    for (int i=0; i<n; i++) {
+                        int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
+                        int globalIdx = indexSet1.template subIndex(inside,faceIdxi,dim);
+                        (*zeroDirichletNodes1)[globalIdx][0] = true;
+                    }
+                }
+            }
+        }
+    }
+
+    std::shared_ptr<Function> constantFunction = std::make_shared<NeumannCondition>(this->parset_.template get<double>("boundary.neumann.sigmaN"));
+
+    auto loadNeumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
+    loadNeumannBoundary->setBoundaryPatch(body1->gridView(), loadNeumannNodes);
+    loadNeumannBoundary->setBoundaryFunction(constantFunction);
+    body1->addBoundaryCondition(loadNeumannBoundary);
+
+    std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary1 = std::make_shared<LeafBoundaryCondition>("dirichlet");
+    zeroDirichletBoundary1->setBoundaryPatch(body1->gridView(), zeroDirichletNodes1);
+    zeroDirichletBoundary1->setBoundaryFunction(zeroFunction);
+    body1->addBoundaryCondition(zeroDirichletBoundary1);
+
+    // body0, body1: natural boundary conditions
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& body = this->contactNetwork_.body(i);
+
+        // Neumann boundary
+        auto neumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
+        auto neumannPatch = std::make_shared<LeafBoundaryPatch>(body->gridView(), true);
+        neumannBoundary->setBoundaryPatch(neumannPatch);
+        neumannBoundary->setBoundaryFunction(neumannFunction);
+        body->addBoundaryCondition(neumannBoundary);
+    }
+}
+
+#include "strikeslipfactory_tmpl.cc"
--- a/dune/tectonic/factories/strikeslipfactory.hh
+++ b/dune/tectonic/factories/strikeslipfactory.hh
+#ifndef DUNE_TECTONIC_FACTORIES_STRIKESLIPFACTORY_HH
+#define DUNE_TECTONIC_FACTORIES_STRIKESLIPFACTORY_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include "contactnetworkfactory.hh"
+
+#include "../problem-data/strikeslipbody.hh"
+#include "../problem-data/grid/trianglegrids.hh"
+#include "../problem-data/grid/trianglegeometry.hh"
+
+template <class HostGridType, class VectorType> class StrikeSlipFactory : public ContactNetworkFactory<HostGridType, VectorType>{
+private:
+    using Base = ContactNetworkFactory<HostGridType, VectorType>;
+
+public:
+    using ContactNetwork = typename Base::ContactNetwork;
+
+private:
+    using GlobalCoords = typename ContactNetwork::LocalVector;
+
+    using LeafGridView = typename ContactNetwork::GridView;
+    using LevelGridView = typename ContactNetwork::GridType::LevelGridView;
+
+    using LevelBoundaryPatch = BoundaryPatch<LevelGridView>;
+    using LeafBoundaryPatch = BoundaryPatch<LeafGridView>;
+
+    using LeafFaces = MyFaces<LeafGridView>;
+    using LevelFaces = MyFaces<LevelGridView>;
+
+    using TriangleGeometry= TriangleGeometry<typename GlobalCoords::field_type>;
+    using WeakeningRegion = typename TriangleGeometry::WeakeningRegion;
+
+    static const int dim = ContactNetwork::dim;
+
+    std::vector<std::shared_ptr<MyBodyData<dim>>> bodyData_;     // material properties of bodies
+
+    std::unique_ptr<GridsConstructor<HostGridType>> gridConstructor_;
+
+    std::vector<std::shared_ptr<TriangleGeometry>> triangleGeometries_;
+
+    std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
+    std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
+
+    std::vector<std::shared_ptr<WeakeningRegion>> weakPatches_;
+
+    std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
+    std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
+
+public:
+    StrikeSlipFactory(const Dune::ParameterTree& parset) :
+        Base(parset, 2, 1),
+        bodyData_(2),
+        triangleGeometries_(2),
+        leafFaces_(2),
+        levelFaces_(2),
+        weakPatches_(2),
+        nonmortarPatch_(1),
+        mortarPatch_(1)
+    {}
+
+    void setBodies();
+    void setLevelBodies();
+    void setCouplings();
+    void setLevelCouplings() {}
+    void setBoundaryConditions();
+
+    auto& weakPatches() {
+        return weakPatches_;
+    }
+
+private:
+    static constexpr Dune::FieldVector<double, MY_DIM> zenith_() {
+        #if MY_DIM == 2
+            return {0, 0};
+        #else
+            return {0, 0, 1};
+        #endif
+    }
+};
+#endif
+
+
--- a/dune/tectonic/factories/strikeslipfactory_tmpl.cc
+++ b/dune/tectonic/factories/strikeslipfactory_tmpl.cc
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+template class StrikeSlipFactory<Grid, Vector>;
--- a/dune/tectonic/factories/threeblocksfactory.cc
+++ b/dune/tectonic/factories/threeblocksfactory.cc
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include "../problem-data/bc.hh"
+#include "../problem-data/myglobalfrictiondata.hh"
+
+#include "../utils/diameter.hh"
+
+#include "threeblocksfactory.hh"
+
+template <class HostGridType, class VectorTEMPLATE>
+void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setBodies() {
+    // set up cuboid geometries
+
+    std::array<double, 3> lengths = {2.0, 1.0, 1.0};
+    std::array<double, 3> heights = {1.0, 1.0, 1.0};
+
+    std::array<GlobalCoords, 3> origins;
+
+    const auto& subParset = this->parset_.sub("boundary.friction.weakening");
+
+#if MY_DIM == 3 // TODO: not implemented
+        double const depth = 0.60;
+
+        origins[0] = {0, 0, 0};
+        origins[1] = {0, origins[0][1] + heights[0], 0};
+        origins[2] = {origins[1][0] + lengths[1], origins[0][1] + heights[0], 0};
+
+        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lengths[0], heights[0], depth);
+        cuboidGeometries_[0]->addWeakeningPatch(subParset, {origins[0][0], origins[0][1]+ heights[0], 0}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0], 0});
+
+        cuboidGeometries_[1] = std::make_shared<CuboidGeometry>(origins[1], lengths[1], heights[1], depth);
+        cuboidGeometries_[1]->addWeakeningPatch(subParset, origins[1], {origins[1][0] + lengths[1], origins[1][1], 0});
+        cuboidGeometries_[1]->addWeakeningPatch(subParset, origins[1], {origins[1][0] + lengths[1], origins[1][1] + heights[1], 0});
+
+        cuboidGeometries_[2] = std::make_shared<CuboidGeometry>(origins[2], lengths[2], heights[2], depth);
+        cuboidGeometries_[2]->addWeakeningPatch(subParset, origins[2], {origins[2][0] + lengths[2], origins[2][1], 0});
+        cuboidGeometries_[2]->addWeakeningPatch(subParset, origins[2], {origins[2][0], origins[2][1] + widths[2], 0});
+#elif MY_DIM == 2
+        origins[0] = {0, 0};
+        origins[1] = {0, origins[0][1] + heights[0]};
+        origins[2] = {origins[1][0] + lengths[1], origins[0][1] + heights[0]};
+
+        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lengths[0], heights[0]);
+        cuboidGeometries_[0]->addWeakeningPatch(subParset, {origins[0][0], origins[0][1]+ heights[0]}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0]});
+
+        cuboidGeometries_[1] = std::make_shared<CuboidGeometry>(origins[1], lengths[1], heights[1]);
+        cuboidGeometries_[1]->addWeakeningPatch(subParset, origins[1], {origins[1][0] + lengths[1], origins[1][1]});
+        cuboidGeometries_[1]->addWeakeningPatch(subParset, {origins[1][0] + lengths[1], origins[1][1]}, {origins[1][0] + lengths[1], origins[1][1] + heights[1]});
+
+        cuboidGeometries_[2] = std::make_shared<CuboidGeometry>(origins[2], lengths[2], heights[2]);
+        cuboidGeometries_[2]->addWeakeningPatch(subParset, origins[2], {origins[2][0] + lengths[2], origins[2][1]});
+        cuboidGeometries_[2]->addWeakeningPatch(subParset, origins[2], {origins[2][0], origins[2][1] + heights[2]});
+#else
+#error CuboidGeometry only supports 2D and 3D!"
+#endif
+
+    // set up reference grids
+    gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(cuboidGeometries_);
+    auto& grids = gridConstructor_->getGrids();
+
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& cuboidGeometry = *cuboidGeometries_[i];
+
+        // define weak patch and refine grid
+        const auto& weakeningRegions = cuboidGeometry.weakeningRegions();
+        for (size_t j=0; j<weakeningRegions.size(); j++) {
+            refine(*grids[i], weakeningRegions[j], this->parset_.template get<double>("boundary.friction.smallestDiameter"), CuboidGeometry::lengthScale());
+        }
+
+        // determine minDiameter and maxDiameter
+        double minDiameter = std::numeric_limits<double>::infinity();
+        double maxDiameter = 0.0;
+        for (auto &&e : elements(grids[i]->leafGridView())) {
+          auto const geometry = e.geometry();
+          auto const diam = diameter(geometry);
+          minDiameter = std::min(minDiameter, diam);
+          maxDiameter = std::max(maxDiameter, diam);
+        }
+        std::cout << "Grid" << i << " min diameter: " << minDiameter << std::endl;
+        std::cout << "Grid" << i << " max diameter: " << maxDiameter << std::endl;
+    }
+
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        this->bodies_[i] = std::make_shared<typename Base::LeafBody>(bodyData_, grids[i]);
+    }
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setLevelBodies() {
+    const size_t maxLevel = std::max({this->bodies_[0]->grid()->maxLevel(), this->bodies_[1]->grid()->maxLevel(), this->bodies_[2]->grid()->maxLevel()});
+
+    for (size_t l=0; l<=maxLevel; l++) {
+        std::vector<size_t> bodyLevels(3, l);
+        this->contactNetwork_.addLevel(bodyLevels, l);
+    }
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setCouplings() {
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& cuboidGeometry = *cuboidGeometries_[i];
+        std::cout << this->bodies_[i]->nVertices() << std::endl;
+        std::cout << (leafFaces_[i] == nullptr) << std::endl;
+        leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), cuboidGeometry);
+        levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), cuboidGeometry);
+    }
+
+    auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
+    std::shared_ptr<typename Base::FrictionCouplingPair::BackEndType> backend = nullptr;
+
+    std::array<std::array<size_t, 2>, 3> couplingIndices;
+
+    nonmortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->lower);
+    mortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->upper);
+    weakPatches_[0] = std::make_shared<WeakeningRegion>(cuboidGeometries_[1]->weakeningRegions()[0]);
+    couplingIndices[0] = {1, 0};
+
+    nonmortarPatch_[1] = std::make_shared<LevelBoundaryPatch>(levelFaces_[2]->lower);
+    mortarPatch_[1] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->upper);
+    weakPatches_[1] = std::make_shared<WeakeningRegion>(cuboidGeometries_[2]->weakeningRegions()[0]);
+    couplingIndices[1] = {2, 0};
+
+    nonmortarPatch_[2] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->right);
+    mortarPatch_[2] = std::make_shared<LevelBoundaryPatch>(levelFaces_[2]->left);
+    weakPatches_[2] = std::make_shared<WeakeningRegion>(cuboidGeometries_[1]->weakeningRegions()[1]);
+    couplingIndices[2] = {1, 2};
+
+    for (size_t i=0; i<this->couplings_.size(); i++) {
+      auto& coupling = this->couplings_[i];
+      coupling = std::make_shared<typename Base::FrictionCouplingPair>();
+
+      coupling->set(couplingIndices[i][0], couplingIndices[i][1], nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
+      coupling->setWeakeningPatch(weakPatches_[i]);
+      coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<GlobalCoords>>(this->parset_.sub("boundary.friction"), *weakPatches_[i], CuboidGeometry::lengthScale()));
+    }
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void ThreeBlocksFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
+    using LeafBoundaryCondition = BoundaryCondition<LeafGridView, dim>;
+
+    using Function = Dune::VirtualFunction<double, double>;
+    std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
+    std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>(this->parset_.template get<double>("boundary.dirichlet.finalVelocity"));
+
+    const double lengthScale = CuboidGeometry::lengthScale();
+
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& body = this->contactNetwork_.body(i);
+        const auto& leafVertexCount = body->nVertices();
+
+        // Neumann boundary
+        std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->gridView()), neumannFunction, "neumann");
+        body->addBoundaryCondition(neumannBoundary);
+
+        // upper Dirichlet Boundary
+        // identify Dirichlet nodes on leaf level
+        if (i>0) {
+            std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount);
+            for (int j=0; j<leafVertexCount; j++) {
+                if (leafFaces_[i]->upper.containsVertex(j))
+                    (*velocityDirichletNodes)[j][0] = true;
+
+                #if MY_DIM == 3 //TODO: wrong, needs revision
+                if (leafFaces_[i]->front.containsVertex(j) || leafFaces_[i]->back.containsVertex(j))
+                    zeroDirichletNodes->at(j)[2] = true;
+                #endif
+            }
+
+            std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+
+            velocityDirichletBoundary->setBoundaryPatch(body->gridView(), velocityDirichletNodes);
+            velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
+            body->addBoundaryCondition(velocityDirichletBoundary);
+        }
+    }
+
+    // lower Dirichlet Boundary
+    const auto& firstBody = this->contactNetwork_.body(0);
+    const auto& firstLeafVertexCount = firstBody->nVertices();
+    std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(firstLeafVertexCount);
+    for (int j=0; j<firstLeafVertexCount; j++) {
+        if (leafFaces_[0]->lower.containsVertex(j)) {
+            for (size_t d=0; d<dim; d++) {
+              (*zeroDirichletNodes)[j][d] = true;
+            }
+        }
+    }
+
+    std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+    zeroDirichletBoundary->setBoundaryPatch(firstBody->gridView(), zeroDirichletNodes);
+
+    std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
+    zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
+    firstBody->addBoundaryCondition(zeroDirichletBoundary);
+}
+
+#include "threeblocksfactory_tmpl.cc"
--- a/dune/tectonic/factories/threeblocksfactory.hh
+++ b/dune/tectonic/factories/threeblocksfactory.hh
+#ifndef SRC_THREEBLOCKSFACTORY_HH
+#define SRC_THREEBLOCKSFACTORY_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include "contactnetworkfactory.hh"
+
+#include "../problem-data/mybody.hh"
+#include "../problem-data/grid/mygrids.hh"
+#include "../problem-data/grid/cuboidgeometry.hh"
+
+template <class HostGridType, class VectorType>
+class ThreeBlocksFactory : public ContactNetworkFactory<HostGridType, VectorType>{
+private:
+    using Base = ContactNetworkFactory<HostGridType, VectorType>;
+
+public:
+    using ContactNetwork = typename Base::ContactNetwork;
+
+private:
+    using GlobalCoords = typename ContactNetwork::LocalVector;
+
+    using LeafGridView = typename ContactNetwork::GridView;
+    using LevelGridView = typename ContactNetwork::GridType::LevelGridView;
+
+    using LevelBoundaryPatch = BoundaryPatch<LevelGridView>;
+    using LeafBoundaryPatch = BoundaryPatch<LeafGridView>;
+
+    using LeafFaces = MyFaces<LeafGridView>;
+    using LevelFaces = MyFaces<LevelGridView>;
+
+    using CuboidGeometry= CuboidGeometry<typename GlobalCoords::field_type>;
+    using WeakeningRegion = typename CuboidGeometry::WeakeningRegion;
+
+    static const int dim = ContactNetwork::dim;
+
+    const std::shared_ptr<MyBodyData<dim>> bodyData_;     // material properties of bodies
+
+    std::unique_ptr<GridsConstructor<HostGridType>> gridConstructor_;
+
+    std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries_;
+
+    std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
+    std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
+
+    std::vector<std::shared_ptr<WeakeningRegion>> weakPatches_;
+
+    std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
+    std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
+
+public:
+    ThreeBlocksFactory(const Dune::ParameterTree& parset) :
+        Base(parset, 3, 3),
+        bodyData_(std::make_shared<MyBodyData<dim>>(this->parset_.sub("body"), this->parset_.template get<double>("gravity"), zenith_())),
+        cuboidGeometries_(3),
+        leafFaces_(3),
+        levelFaces_(3),
+        weakPatches_(3),
+        nonmortarPatch_(3),
+        mortarPatch_(3)
+    {}
+
+    void setBodies();
+    void setLevelBodies();
+    void setCouplings();
+    void setLevelCouplings() {}
+    void setBoundaryConditions();
+
+    auto& weakPatches() {
+        return weakPatches_;
+    }
+private:
+    static constexpr Dune::FieldVector<double, MY_DIM> zenith_() {
+        #if MY_DIM == 2
+        return {0, 1};
+        #elif MY_DIM == 3
+        return {0, 1, 0};
+        #endif
+    }
+};
+#endif
+
+
--- a/dune/tectonic/factories/threeblocksfactory_tmpl.cc
+++ b/dune/tectonic/factories/threeblocksfactory_tmpl.cc
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+template class ThreeBlocksFactory<Grid, Vector>;
--- a/dune/tectonic/factories/twoblocksfactory.cc
+++ b/dune/tectonic/factories/twoblocksfactory.cc
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <string>
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include "../problem-data/bc.hh"
+#include "../problem-data/myglobalfrictiondata.hh"
+
+#include "../utils/diameter.hh"
+
+#include "twoblocksfactory.hh"
+
+template <class HostGridType, class VectorTEMPLATE>
+void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setBodies() {
+    // set up cuboid geometries
+
+    std::array<double, 2> lengths = {this->parset_.template get<double>("body0.length"), this->parset_.template get<double>("body1.length")};
+    std::array<double, 2> heights = {this->parset_.template get<double>("body0.height"), this->parset_.template get<double>("body1.height")};
+
+    std::array<GlobalCoords, 2> origins;
+
+    const auto& frictionParset = this->parset_.sub("boundary.friction");
+
+#if MY_DIM == 3 // TODO: not implemented
+        std::array<double, 2> depths = {this->parset_.template get<double>("body0.depth"), this->parset_.template get<double>("body1.depth")};
+
+        origins[0] = {0, 0, 0};
+        origins[1] = {lengths[0]/2.0, origins[0][1] + heights[0], 0};
+
+        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lengths[0], heights[0], depths[0]);
+        cuboidGeometries_[0]->addWeakeningPatch(frictionParset, {origins[0][0], origins[0][1]+ heights[0], 0}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0], 0});
+
+        cuboidGeometries_[1] = std::make_shared<CuboidGeometry>(origins[1], lengths[1], heights[1], depths[1]);
+        cuboidGeometries_[1]->addWeakeningPatch(frictionParset, origins[1], {origins[1][0] + lengths[1], origins[1][1], 0});
+#elif MY_DIM == 2
+        origins[0] = {-lengths[0]/2.0, 0};
+        origins[1] = {-lengths[1]/2.0, origins[0][1] + heights[0]};
+
+        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lengths[0], heights[0]);
+        cuboidGeometries_[0]->addWeakeningPatch(frictionParset, {origins[0][0], origins[0][1]+ heights[0]}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0]});
+
+        cuboidGeometries_[1] = std::make_shared<CuboidGeometry>(origins[1], lengths[1], heights[1]);
+        cuboidGeometries_[1]->addWeakeningPatch(frictionParset, origins[1], {origins[1][0] + lengths[1], origins[1][1]});
+#else
+#error CuboidGeometry only supports 2D and 3D!"
+#endif
+
+    const auto gravity = this->parset_.template get<double>("general.gravity");
+
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        std::string body_str = "body" + std::to_string(i);
+        const auto& cuboidGeometry = *cuboidGeometries_[i];
+
+        std::array<unsigned int, 2> initialElements = {
+            { this->parset_.template get<size_t>("body" + std::to_string(i) + ".initialXElements"),
+              this->parset_.template get<size_t>("body" + std::to_string(i) + ".initialYElements") } };
+
+        // set up reference grid
+        DefaultCuboidGridConstructor<HostGridType> gridConstructor(cuboidGeometry, initialElements);
+        gridConstructor.refine(this->parset_.template get<size_t>(body_str + ".refinements"));
+
+        bodyData_[i] = std::make_shared<MyBodyData<dim>>(this->parset_.sub(body_str), gravity, zenith_());
+        this->bodies_[i] = std::make_shared<typename Base::LeafBody>(bodyData_[i], gridConstructor.grid());
+    }
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setLevelBodies() {
+    const size_t maxLevel = std::max({this->bodies_[0]->grid()->maxLevel(), this->bodies_[1]->grid()->maxLevel()});
+
+    for (size_t l=0; l<=maxLevel; l++) {
+        std::vector<size_t> bodyLevels(2, l);
+        this->contactNetwork_.addLevel(bodyLevels, l);
+    }
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setCouplings() {
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& cuboidGeometry = *cuboidGeometries_[i];
+        leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), cuboidGeometry);
+        levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), cuboidGeometry);
+    }
+
+    auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
+
+    nonmortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->lower);
+    mortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->upper);
+    weakPatches_[0] = std::make_shared<WeakeningRegion>(cuboidGeometries_[1]->weakeningRegions()[0]);
+
+    auto& coupling = this->couplings_[0];
+    coupling = std::make_shared<typename Base::FrictionCouplingPair>();
+
+    coupling->set(1, 0, nonmortarPatch_[0], mortarPatch_[0], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, nullptr);
+    coupling->setWeakeningPatch(weakPatches_[0]);
+    coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<GlobalCoords>>(this->parset_.sub("boundary.friction"), *weakPatches_[0], CuboidGeometry::lengthScale()));
+}
+
+template <class HostGridType, class VectorTEMPLATE>
+void TwoBlocksFactory<HostGridType, VectorTEMPLATE>::setBoundaryConditions() {
+    using LeafBoundaryCondition = BoundaryCondition<LeafGridView, dim>;
+
+    using Function = Dune::VirtualFunction<double, double>;
+    std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
+    //std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletConditionLinearLoading>(this->parset_.template get<double>("boundary.dirichlet.finalVelocity"), 0.25);
+    std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>(this->parset_.template get<double>("boundary.dirichlet.finalVelocity"), 0.25);
+    //std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityStepDirichletCondition>(this->parset_.template get<double>("boundary.dirichlet.finalVelocity"), 10*this->parset_.template get<double>("boundary.dirichlet.finalVelocity"), 0.25, 0.5);
+
+    const double lengthScale = CuboidGeometry::lengthScale();
+
+    // body0: Dirichlet boundary (lower)
+    const auto& body0 = this->contactNetwork_.body(0);
+    const auto& leafVertexCount0 = body0->nVertices();
+    std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount0);
+    for (int j=0; j<leafVertexCount0; j++) {
+        if (leafFaces_[0]->lower.containsVertex(j)) {
+            for (size_t d=0; d<dim; d++) {
+              (*zeroDirichletNodes)[j][d] = true;
+            }
+        }
+    }
+
+    std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+    zeroDirichletBoundary->setBoundaryPatch(body0->gridView(), zeroDirichletNodes);
+
+    std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
+    zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
+    body0->addBoundaryCondition(zeroDirichletBoundary);
+
+    // body1: Dirichlet boundary (upper)
+    const auto& body1 = this->contactNetwork_.body(1);
+    const auto& leafVertexCount1 = body1->nVertices();
+    std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
+    for (size_t j=0; j<leafVertexCount1; j++) {
+        if (leafFaces_[1]->upper.containsVertex(j))
+            (*velocityDirichletNodes)[j][0] = true;
+
+        #if MY_DIM == 3 //TODO: wrong, needs revision
+        if (leafFaces_[1]->front.containsVertex(j) || leafFaces_[1]->back.containsVertex(j))
+            zeroDirichletNodes->at(j)[2] = true;
+        #endif
+    }
+
+    std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
+
+    velocityDirichletBoundary->setBoundaryPatch(body1->gridView(), velocityDirichletNodes);
+    velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
+    body1->addBoundaryCondition(velocityDirichletBoundary);
+
+    // body0, body1: natural boundary conditions
+    for (size_t i=0; i<this->bodyCount_; i++) {
+        const auto& body = this->contactNetwork_.body(i);
+
+        // Neumann boundary
+        auto neumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
+        auto neumannPatch = std::make_shared<LeafBoundaryPatch>(body->gridView(), true);
+        neumannBoundary->setBoundaryPatch(neumannPatch);
+        neumannBoundary->setBoundaryFunction(neumannFunction);
+        body->addBoundaryCondition(neumannBoundary);
+    }
+
+    // body1: Neumann boundary (upper)
+    // normal load
+    std::shared_ptr<Function> constantFunction = std::make_shared<NeumannCondition>(this->parset_.template get<double>("boundary.neumann.sigmaN"));
+
+    std::shared_ptr<Dune::BitSetVector<dim>> loadNeumannNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
+    for (size_t j=0; j<leafVertexCount1; j++) {
+        if (leafFaces_[1]->upper.containsVertex(j))
+            (*loadNeumannNodes)[j][1] = true;
+
+        #if MY_DIM == 3 //TODO: wrong, needs revision
+        if (leafFaces_[1]->front.containsVertex(j) || leafFaces_[1]->back.containsVertex(j))
+            zeroDirichletNodes->at(j)[2] = true;
+        #endif
+    }
+
+    auto loadNeumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
+    loadNeumannBoundary->setBoundaryPatch(body1->gridView(), loadNeumannNodes);
+    loadNeumannBoundary->setBoundaryFunction(constantFunction);
+    body1->addBoundaryCondition(loadNeumannBoundary);
+}
+
+#include "twoblocksfactory_tmpl.cc"
--- a/dune/tectonic/factories/twoblocksfactory.hh
+++ b/dune/tectonic/factories/twoblocksfactory.hh
+#ifndef DUNE_TECTONIC_FACTORIES_TWOBLOCKSFACTORY_HH
+#define DUNE_TECTONIC_FACTORIES_TWOBLOCKSFACTORY_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include "contactnetworkfactory.hh"
+
+#include "../problem-data/mybody.hh"
+#include "../problem-data/grid/cuboidgeometry.hh"
+#include "../problem-data/grid/cuboidgridconstructor.hh"
+
+template <class HostGridType, class VectorType> class TwoBlocksFactory : public ContactNetworkFactory<HostGridType, VectorType>{
+private:
+    using Base = ContactNetworkFactory<HostGridType, VectorType>;
+
+public:
+    using ContactNetwork = typename Base::ContactNetwork;
+
+private:
+    using GlobalCoords = typename ContactNetwork::LocalVector;
+
+    using LeafGridView = typename ContactNetwork::GridView;
+    using LevelGridView = typename ContactNetwork::GridType::LevelGridView;
+
+    using LevelBoundaryPatch = BoundaryPatch<LevelGridView>;
+    using LeafBoundaryPatch = BoundaryPatch<LeafGridView>;
+
+    using LeafFaces = MyFaces<LeafGridView>;
+    using LevelFaces = MyFaces<LevelGridView>;
+
+    using CuboidGeometry= CuboidGeometry<typename GlobalCoords::field_type>;
+    using WeakeningRegion = typename CuboidGeometry::WeakeningRegion;
+
+    static const int dim = ContactNetwork::dim;
+
+    std::vector<std::shared_ptr<MyBodyData<dim>>> bodyData_;     // material properties of bodies
+
+    std::vector<std::shared_ptr<CuboidGeometry>> cuboidGeometries_;
+
+    std::vector<std::shared_ptr<LeafFaces>> leafFaces_;
+    std::vector<std::shared_ptr<LevelFaces>> levelFaces_;
+
+    std::vector<std::shared_ptr<WeakeningRegion>> weakPatches_;
+
+    std::vector<std::shared_ptr<LevelBoundaryPatch>> nonmortarPatch_;
+    std::vector<std::shared_ptr<LevelBoundaryPatch>> mortarPatch_;
+
+public:
+    TwoBlocksFactory(const Dune::ParameterTree& parset) :
+        Base(parset, 2, 1),
+        bodyData_(2),
+        cuboidGeometries_(2),
+        leafFaces_(2),
+        levelFaces_(2),
+        weakPatches_(2),
+        nonmortarPatch_(1),
+        mortarPatch_(1)
+    {}
+
+    void setBodies();
+    void setLevelBodies();
+    void setCouplings();
+    void setLevelCouplings() {}
+    void setBoundaryConditions();
+
+    auto& weakPatches() {
+        return weakPatches_;
+    }
+
+private:
+    static constexpr Dune::FieldVector<double, MY_DIM> zenith_() {
+        #if MY_DIM == 2
+        return {0, 1};
+        #elif MY_DIM == 3
+        return {0, 1, 0};
+        #endif
+    }
+};
+#endif
+
+
--- a/dune/tectonic/factories/twoblocksfactory_tmpl.cc
+++ b/dune/tectonic/factories/twoblocksfactory_tmpl.cc
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+template class TwoBlocksFactory<Grid, Vector>;
--- a/dune/tectonic/globalratestatefriction.hh
+++ b/dune/tectonic/globalratestatefriction.hh
-#ifndef DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
-#define DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
-
-#include <vector>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fvector.hh>
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-
-#include <dune/tectonic/geocoordinate.hh>
-#include <dune/tectonic/globalfrictiondata.hh>
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/index-in-sorted-range.hh>
-
-template <class Matrix, class Vector, class ScalarFriction, class GridView>
-class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
-public:
-  using GlobalFriction<Matrix, Vector>::block_size;
-  using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;
-
-private:
-  using typename GlobalFriction<Matrix, Vector>::ScalarVector;
-
-public:
-  GlobalRateStateFriction(BoundaryPatch<GridView> const &frictionalBoundary,
-                          GlobalFrictionData<block_size> const &frictionInfo,
-                          ScalarVector const &weights,
-                          ScalarVector const &weightedNormalStress)
-      : restrictions(), localToGlobal(), zeroFriction() {
-    auto const gridView = frictionalBoundary.gridView();
-    Dune::MultipleCodimMultipleGeomTypeMapper<
-        GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
-    for (auto it = gridView.template begin<block_size>();
-         it != gridView.template end<block_size>(); ++it) {
-      auto const i = vertexMapper.index(*it);
-
-      if (not frictionalBoundary.containsVertex(i))
-        continue;
-
-      localToGlobal.emplace_back(i);
-      restrictions.emplace_back(weights[i], weightedNormalStress[i],
-                                frictionInfo(geoToPoint(it->geometry())));
-    }
-    assert(restrictions.size() == frictionalBoundary.numVertices());
-    assert(localToGlobal.size() == frictionalBoundary.numVertices());
-  }
-
-  void updateAlpha(ScalarVector const &alpha) override {
-    for (size_t j = 0; j < restrictions.size(); ++j)
-      restrictions[j].updateAlpha(alpha[localToGlobal[j]]);
-  }
-
-  /*
-    Return a restriction of the outer function to the i'th node.
-  */
-  LocalNonlinearity const &restriction(size_t i) const override {
-    auto const index = indexInSortedRange(localToGlobal, i);
-    if (index == localToGlobal.size())
-      return zeroFriction;
-    return restrictions[index];
-  }
-
-private:
-  std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions;
-  std::vector<size_t> localToGlobal;
-  WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction;
-};
-#endif
--- a/src/gridselector.hh
+++ b/src/gridselector.hh
--- a/dune/tectonic/io/CMakeLists.txt
+++ b/dune/tectonic/io/CMakeLists.txt
+add_subdirectory("hdf5")
+
+add_custom_target(tectonic_dune_io SOURCES
+  hdf5-bodywriter.hh
+  hdf5-writer.hh
+  io-handler.hh
+  uniform-grid-writer.cc
+  vtk-bodywriter.hh
+  vtk.hh
+)
+
+#install headers
+install(FILES
+  hdf5-bodywriter.hh
+  hdf5-levelwriter.hh
+  io-handler.hh
+  vtk-bodywriter.hh
+  vtk.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
--- a/dune/tectonic/io/hdf5-bodywriter.hh
+++ b/dune/tectonic/io/hdf5-bodywriter.hh
+#ifndef SRC_HDF5_BODYWRITER_HH
+#define SRC_HDF5_BODYWRITER_HH
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/fufem/hdf5/file.hh>
+
+#include "hdf5/frictionalboundary-writer.hh"
+//#include "hdf5/patchinfo-writer.hh"
+
+template <class ProgramState, class VertexBasis, class GridView>
+class HDF5BodyWriter {
+private:
+    using Vector = typename ProgramState::Vector;
+    using ScalarVector = typename ProgramState::ScalarVector;
+    //using PatchInfoWriter = PatchInfoWriter<ProgramState, VertexBasis, GridView>;
+
+public:
+    using VertexCoordinates = typename ProgramState::Vector;
+    using Patch = typename FrictionalBoundaryWriter<GridView>::Patch;
+    //using WeakPatches = std::vector<const ConvexPolyhedron<LocalVector>* >;
+
+    using LocalVector = typename VertexCoordinates::block_type;
+
+    //friend class HDF5LevelWriter<ProgramState, VertexBasis, GridView>;
+
+    HDF5BodyWriter(
+            HDF5::File& file,
+            const size_t bodyID,
+            const VertexCoordinates& vertexCoordinates,
+            const VertexBasis& vertexBasis,
+            const Patch& frictionPatch) :
+           // const WeakPatches& weakPatches) :
+        id_(bodyID),
+        /*#if MY_DIM == 3
+                patchInfoWriters_(patchCount_),
+        #endif*/
+        group_(file, "body"+std::to_string(id_)) {
+
+        frictionBoundaryWriter_ = std::make_unique<FrictionalBoundaryWriter<GridView>>(group_, vertexCoordinates, frictionPatch);
+
+        /*#if MY_DIM == 3
+            patchInfoWriters_[i] = std::make_unique<PatchInfoWriter>(file_, vertexBasis, *frictionPatches[i], *weakPatches[i], i);
+        #endif*/
+    }
+
+    template <class Friction>
+    void reportSolution(const ProgramState& programState, const Vector& v, const ScalarVector& frictionCoeff, const Friction& friction) {
+        frictionBoundaryWriter_->write(programState.timeStep, programState.u[id_], v, programState.alpha[id_], frictionCoeff, friction);
+
+        /*#if MY_DIM == 3
+        patchInfoWriters_[i]->write(programState);
+        #endif*/
+    }
+
+    void reportWeightedNormalStress(const ProgramState& programState) {
+        frictionBoundaryWriter_->writeWeightedNormalStress(programState.timeStep, programState.weightedNormalStress[id_], programState.weights[id_]);
+    }
+
+    auto id() const {
+        return id_;
+    }
+
+private:
+    const size_t id_;
+    HDF5::Group group_;
+
+/*#if MY_DIM == 3
+    std::vector<std::unique_ptr<PatchInfoWriter>> patchInfoWriters_;
+#endif*/
+
+    std::unique_ptr<FrictionalBoundaryWriter<GridView>> frictionBoundaryWriter_;
+};
+#endif
--- a/dune/tectonic/io/hdf5-writer.hh
+++ b/dune/tectonic/io/hdf5-writer.hh
+#ifndef SRC_HDF5_WRITER_HH
+#define SRC_HDF5_WRITER_HH
+
+#include <vector>
+
+#include "hdf5/iteration-writer.hh"
+#include "hdf5/time-writer.hh"
+#include "hdf5-bodywriter.hh"
+
+#include <dune/fufem/hdf5/file.hh>
+
+#include "../utils/debugutils.hh"
+
+template <class ProgramState, class VertexBasis, class GridView>
+class HDF5Writer {
+public:
+    using HDF5BodyWriter = HDF5BodyWriter<ProgramState, VertexBasis, GridView>;
+    using VertexCoordinates = std::vector<typename HDF5BodyWriter::VertexCoordinates>;
+    using VertexBases = std::vector<const VertexBasis*>;
+    using FrictionPatches = std::vector<const typename HDF5BodyWriter::Patch*>;
+    using ScalarVector = typename ProgramState::ScalarVector;
+    //using WeakPatches = std::vector<typename HDF5BodyWriter::WeakPatches>;
+
+    //friend class HDF5NetworkWriter<ProgramState, VertexBasis, GridView>;
+
+    HDF5Writer(HDF5::File& file,
+             const VertexCoordinates& vertexCoordinates,
+             const VertexBases& vertexBases,
+             const FrictionPatches& frictionPatches)
+             //const WeakPatches& weakPatches)
+      : file_(file),
+        frictionPatches_(frictionPatches),
+        iterationWriter_(file_),
+        timeWriter_(file_) {
+
+        for (size_t i=0; i<vertexCoordinates.size(); i++) {
+            if (frictionPatches_[i]->numVertices() > 0)
+                bodyWriters_.push_back(std::make_unique<HDF5BodyWriter>(file_, i, vertexCoordinates[i], *vertexBases[i], *frictionPatches_[i])); //, weakPatches[i]));
+        }
+    }
+
+    template <class ContactNetwork, class Friction>
+    void reportSolution(const ProgramState& programState, const ContactNetwork& contactNetwork, const Friction& friction) {
+        timeWriter_.write(programState);
+
+        // extract relative velocities
+        using Vector = typename ProgramState::Vector;
+        Vector mortarV;
+        contactNetwork.nBodyAssembler().nodalToTransformed(programState.v, mortarV);
+
+        //printRegularityTruncation(friction, mortarV);
+
+        std::vector<Vector> v_rel;
+        split(mortarV, v_rel);
+
+        using ScalarVector = typename ProgramState::ScalarVector;
+        const auto frictionCoeff = friction.coefficientOfFriction(mortarV);
+
+        std::vector<ScalarVector> splitCoeff;
+        split(frictionCoeff, splitCoeff);
+
+        for (size_t i=0; i<bodyWriters_.size(); i++) {
+            auto bodyID = bodyWriters_[i]->id();
+            bodyWriters_[i]->reportSolution(programState, v_rel[bodyID], splitCoeff[bodyID], friction);
+        }
+    }
+
+    void reportIterations(const ProgramState& programState, const IterationRegister& iterationCount) {
+        iterationWriter_.write(programState.timeStep, iterationCount);
+    }
+
+
+    void reportWeightedNormalStress(const ProgramState& programState) {
+        for (size_t i=0; i<bodyWriters_.size(); i++) {
+            auto bodyID = bodyWriters_[i]->id();
+            bodyWriters_[i]->reportWeightedNormalStress(programState);
+        }
+    }
+
+    template <class VectorType>
+    void split(const VectorType& v, std::vector<VectorType>& splitV) const {
+        const auto nBodies = frictionPatches_.size();
+
+        size_t globalIdx = 0;
+        splitV.resize(nBodies);
+
+        for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
+            const auto& bodyNodes = *frictionPatches_[bodyIdx]->getVertices();
+
+            auto& splitV_ = splitV[bodyIdx];
+            splitV_.resize(bodyNodes.size());
+
+            for (size_t i=0; i<splitV_.size(); i++) {
+                if (toBool(bodyNodes[i])) {
+                    splitV_[i] = v[globalIdx];
+                }
+                globalIdx++;
+            }
+        }
+    }
+
+private:
+  HDF5::File& file_;
+  const FrictionPatches& frictionPatches_;
+
+  IterationWriter iterationWriter_;
+  TimeWriter<ProgramState> timeWriter_;
+
+  std::vector<std::unique_ptr<HDF5BodyWriter>> bodyWriters_;
+};
+#endif
--- a/dune/tectonic/io/hdf5/CMakeLists.txt
+++ b/dune/tectonic/io/hdf5/CMakeLists.txt
+add_custom_target(tectonic_dune_io_hdf5 SOURCES
+  frictionalboundary-writer.hh
+  #frictionalboundary-writer.cc
+  iteration-writer.hh 
+  #iteration-writer.cc
+  #patchinfo-writer.hh
+  #patchinfo-writer.cc
+  restartbody-io.hh
+  restart-io.hh
+  restrict.hh
+  surface-writer.hh
+  time-writer.hh
+)
+
+#install headers
+install(FILES
+  frictionalboundary-writer.hh
+  iteration-writer.hh 
+  #patchinfo-writer.hh
+  restartbody-io.hh
+  restart-io.hh
+  restrict.hh
+  surface-writer.hh
+  time-writer.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
--- a/dune/tectonic/io/hdf5/frictionalboundary-writer.hh
+++ b/dune/tectonic/io/hdf5/frictionalboundary-writer.hh
+#ifndef SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
+#define SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/hdf5/file.hh>
+#include <dune/fufem/hdf5/sequenceio.hh>
+#include <dune/fufem/hdf5/singletonwriter.hh>
+
+#include "restrict.hh"
+
+template <class GridView> class FrictionalBoundaryWriter {
+public:
+  using Patch = BoundaryPatch<GridView>;
+
+  template <class Vector>
+  FrictionalBoundaryWriter(HDF5::Group& group, const Vector& vertexCoordinates,
+                           const Patch& frictionalBoundary) : group_(group),
+      frictionalBoundary_(frictionalBoundary),
+      frictionalBoundaryDisplacementWriter_(group_, "displacement",
+                                            frictionalBoundary.numVertices(),
+                                            Vector::block_type::dimension),
+      frictionalBoundaryVelocityWriter_(group_, "velocity",
+                                        frictionalBoundary.numVertices(),
+                                        Vector::block_type::dimension),
+      frictionalBoundaryStateWriter_(group_, "state",
+                                     frictionalBoundary.numVertices()),
+      frictionalBoundaryCoefficientWriter_(group_, "coefficient",
+                                           frictionalBoundary.numVertices()),
+      frictionalBoundaryWeightsWriter_(group_, "weights",
+                                           frictionalBoundary.numVertices()),
+      frictionalBoundaryWeightedNormalStressWriter_(group_, "weightedNormalStress",
+                                           frictionalBoundary.numVertices()){
+
+  auto const frictionalBoundaryCoordinates =
+      restrictToSurface(vertexCoordinates, frictionalBoundary);
+
+  HDF5::SingletonWriter<2> frictionalBoundaryCoordinateWriter(
+      group_, "coordinates", frictionalBoundaryCoordinates.size(),
+      Vector::block_type::dimension);
+  setEntry(frictionalBoundaryCoordinateWriter, frictionalBoundaryCoordinates);
+}
+
+  template <class Vector, class ScalarVector, class Friction>
+  void write(const size_t timeStep, const Vector& u, const Vector& v, const ScalarVector& alpha, const ScalarVector& frictionCoeff, const Friction& friction) {
+
+      auto const frictionalBoundaryDisplacements = restrictToSurface(u, frictionalBoundary_);
+      addEntry(frictionalBoundaryDisplacementWriter_, timeStep, frictionalBoundaryDisplacements);
+
+      auto const frictionalBoundaryVelocities = restrictToSurface(v, frictionalBoundary_);
+      addEntry(frictionalBoundaryVelocityWriter_, timeStep, frictionalBoundaryVelocities);
+
+      auto const frictionalBoundaryStates = restrictToSurface(alpha, frictionalBoundary_);
+      addEntry(frictionalBoundaryStateWriter_, timeStep, frictionalBoundaryStates);
+
+      auto const frictionalBoundaryCoefficient = restrictToSurface(frictionCoeff, frictionalBoundary_);
+      addEntry(frictionalBoundaryCoefficientWriter_, timeStep, frictionalBoundaryCoefficient);
+
+     // auto const regularityTruncation = printRegularityTruncation()
+     // addEntry(frictionalBoundaryRegularityTruncationWriter_, timeStep, regularityTruncation);
+    }
+
+  template <class ScalarVector>
+  void writeWeightedNormalStress(const size_t timeStep, const ScalarVector& weightedNormalStress, const ScalarVector& weights) {
+
+      auto const frictionalBoundaryWeights = restrictToSurface(weights, frictionalBoundary_);
+      addEntry(frictionalBoundaryWeightsWriter_, timeStep, frictionalBoundaryWeights);
+
+      auto const frictionalBoundaryWeightedNormalStress = restrictToSurface(weightedNormalStress, frictionalBoundary_);
+      addEntry(frictionalBoundaryWeightedNormalStressWriter_, timeStep, frictionalBoundaryWeightedNormalStress);
+  }
+
+private:
+  HDF5::Group& group_;
+
+  const Patch& frictionalBoundary_;
+
+  HDF5::SequenceIO<2> frictionalBoundaryDisplacementWriter_;
+  HDF5::SequenceIO<2> frictionalBoundaryVelocityWriter_;
+  HDF5::SequenceIO<1> frictionalBoundaryStateWriter_;
+  HDF5::SequenceIO<1> frictionalBoundaryCoefficientWriter_;
+  HDF5::SequenceIO<1> frictionalBoundaryWeightsWriter_;
+  HDF5::SequenceIO<1> frictionalBoundaryWeightedNormalStressWriter_;
+  //HDF5::SequenceIO<0> frictionalBoundaryRegularityTruncationWriter_;
+};
+
+#endif
No results found