Compare revisions

fbd31be5 · fbd31be5 · fbd31be5 · fbd31be5 · fbd31be5 · fbd31be5
--- a/src/factories/contactnetworkfactory.hh
+++ b/src/factories/contactnetworkfactory.hh
+#ifndef SRC_CONTACTNETWORKFACTORY_HH
+#define SRC_CONTACTNETWORKFACTORY_HH
+
+#include <dune/common/parametertree.hh>
+
+#include "../data-structures/contactnetwork.hh"
+
+template <class HostGridType, class VectorType>
+class ContactNetworkFactory {
+public:
+    using ContactNetwork = ContactNetwork<HostGridType, VectorType>;
+
+protected:
+    using LeafBody = typename ContactNetwork::LeafBody;
+    using FrictionCouplingPair = typename ContactNetwork::FrictionCouplingPair;
+
+    const Dune::ParameterTree& parset_;
+    const size_t bodyCount_;
+    const size_t couplingCount_;
+
+    std::vector<std::shared_ptr<LeafBody>> bodies_;
+    std::vector<std::shared_ptr<FrictionCouplingPair>> couplings_;
+
+    ContactNetwork contactNetwork_;
+
+private:
+    virtual void setBodies() = 0;
+    virtual void setLevelBodies() = 0;
+    virtual void setCouplings() = 0;
+    virtual void setLevelCouplings() = 0;
+    virtual void setBoundaryConditions() = 0;
+
+public:
+    ContactNetworkFactory(const Dune::ParameterTree& parset, size_t bodyCount, size_t couplingCount) :
+        parset_(parset),
+        bodyCount_(bodyCount),
+        couplingCount_(couplingCount),
+        bodies_(bodyCount_),
+        couplings_(couplingCount_),
+        contactNetwork_(bodyCount_, couplingCount_) {}
+
+    void build() {
+        setBodies();
+        contactNetwork_.setBodies(bodies_);
+
+        setCouplings();
+        contactNetwork_.setCouplings(couplings_);
+
+        setLevelBodies();
+        setLevelCouplings();
+
+        setBoundaryConditions();      
+    }
+
+    ContactNetwork& contactNetwork() {
+        return contactNetwork_;
+    }
+};
+#endif
--- a/src/factories/levelcontactnetworkfactory.hh
+++ b/src/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/src/factories/stackedblocksfactory.cc
+++ b/src/factories/stackedblocksfactory.cc
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include "../multi-body-problem-data/bc.hh"
+#include "../multi-body-problem-data/grid/cuboidgeometry.hh"
+#include "../multi-body-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
+    double const length = 1.00;
+    double const height = 0.27;
+    double const weakLength = 0.60;
+
+    std::vector<GlobalCoords> origins(this->bodyCount_);
+
+    const auto& subParset = this->parset_.sub("boundary.friction.weakening");
+
+#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, origin[1]};
+            return h;
+        };
+
+        origins[0] = {0,0};
+        GlobalCoords lowerWeakOrigin = {0.2, 0};
+        GlobalCoords upperWeakOrigin = {0.2, 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();
+            GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
+            GlobalCoords upperWeakOrigin_i = upperWeakOrigin + origins[i];
+
+            cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height);
+
+
+            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);
+        cuboidGeometries_[idx]->addWeakeningPatch(subParset, lowerWeakOrigin, weakBound(lowerWeakOrigin));
+#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 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]->upper);
+      mortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i+1]->lower);
+      weakPatches_[i] = std::make_shared<WeakeningRegion>(cuboidGeometries_[i]->weakeningRegions()[0]);
+
+      coupling->set(i, 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 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]--;
+        }
+    }
+}
+
+
+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>();
+
+    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;
+
+                #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);
+        } 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;
+            }
+        }
+    }
+
+    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/src/factories/stackedblocksfactory.hh
+++ b/src/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 "../multi-body-problem-data/mybody.hh"
+#include "../multi-body-problem-data/grid/mygrids.hh"
+#include "../multi-body-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::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:
+    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_, 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();
+
+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/src/factories/threeblocksfactory.cc
+++ b/src/factories/threeblocksfactory.cc
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include <dune/contact/projections/normalprojection.hh>
+
+#include "../multi-body-problem-data/bc.hh"
+#include "../multi-body-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 = {3.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] + widths[0], 0};
+        origins[2] = {origins[1][0] + lengths[1], origins[0][1] + widths[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_[0]->upper);
+    mortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->lower);
+    weakPatches_[0] = std::make_shared<WeakeningRegion>(cuboidGeometries_[0]->weakeningRegions()[0]);
+    couplingIndices[0] = {0, 1};
+
+    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>();
+
+    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/src/factories/threeblocksfactory.hh
+++ b/src/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 "../multi-body-problem-data/mybody.hh"
+#include "../multi-body-problem-data/grid/mygrids.hh"
+#include "../multi-body-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_, zenith_())),
+        cuboidGeometries_(3),
+        leafFaces_(3),
+        levelFaces_(3),
+        weakPatches_(3),
+        nonmortarPatch_(3),
+        mortarPatch_(3)
+    {}
+
+    void setBodies();
+    void setLevelBodies();
+    void setCouplings();
+    void setLevelCouplings() {}
+    void setBoundaryConditions();
+
+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/factories/threeblocksfactory_tmpl.cc
+++ b/src/factories/threeblocksfactory_tmpl.cc
+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+template class ThreeBlocksFactory<Grid, Vector>;
--- a/src/frictioncouplingpair.hh
+++ b/src/frictioncouplingpair.hh
+#ifndef SRC_FRICTIONCOUPLINGPAIR_HH
+#define SRC_FRICTIONCOUPLINGPAIR_HH
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/tectonic/globalfrictiondata.hh>
+
+template <class GridType, class LocalVectorType, class field_type = double>
+class FrictionCouplingPair : public Dune::Contact::CouplingPair<GridType, GridType, field_type>{
+private:
+    static const int dim = GridType::dimensionworld;
+
+    using Base = Dune::Contact::CouplingPair<GridType,GridType,field_type>;
+    using LocalVector = LocalVectorType;
+
+    // friction data
+    std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch_;
+    std::shared_ptr<GlobalFrictionData<dim>> frictionData_;
+
+public:
+  void setWeakeningPatch(std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch) {
+      weakeningPatch_ = weakeningPatch;
+  }
+
+  void setFrictionData(std::shared_ptr<GlobalFrictionData<dim>> frictionData) {
+      frictionData_ = frictionData;
+  }
+
+  const auto& weakeningPatch() const {
+      return *weakeningPatch_;
+  }
+
+  const auto& frictionData() const {
+      return *frictionData_;
+  }
+};
+#endif
--- a/src/hdf5-writer.hh
+++ b/src/hdf5-writer.hh
-#ifndef SRC_HDF5_WRITER_HH
-#define SRC_HDF5_WRITER_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/iteration-writer.hh"
-#include "hdf5/patchinfo-writer.hh"
-#include "hdf5/surface-writer.hh"
-#include "hdf5/time-writer.hh"
-
-template <class ProgramState, class VertexBasis, class GridView>
-class HDF5Writer {
-private:
-  using Vector = typename ProgramState::Vector;
-  using Patch = BoundaryPatch<GridView>;
-
-  using LocalVector = typename Vector::block_type;
-
-public:
-  HDF5Writer(HDF5::Grouplike &file, Vector const &vertexCoordinates,
-             VertexBasis const &vertexBasis, Patch const &surface,
-             Patch const &frictionalBoundary,
-             ConvexPolyhedron<LocalVector> const &weakPatch)
-      : file_(file),
-        iterationWriter_(file_),
-        timeWriter_(file_),
-#if MY_DIM == 3
-        patchInfoWriter_(file_, vertexBasis, frictionalBoundary, weakPatch),
-#endif
-        surfaceWriter_(file_, vertexCoordinates, surface),
-        frictionalBoundaryWriter_(file_, vertexCoordinates,
-                                  frictionalBoundary) {
-  }
-
-  template <class Friction>
-  void reportSolution(ProgramState const &programState,
-                      // for the friction coefficient
-                      Friction &friction) {
-    timeWriter_.write(programState);
-#if MY_DIM == 3
-    patchInfoWriter_.write(programState);
-#endif
-    surfaceWriter_.write(programState);
-    frictionalBoundaryWriter_.write(programState, friction);
-  }
-
-  void reportIterations(ProgramState const &programState,
-                        IterationRegister const &iterationCount) {
-    iterationWriter_.write(programState.timeStep, iterationCount);
-  }
-
-private:
-  HDF5::Grouplike &file_;
-
-  IterationWriter iterationWriter_;
-  TimeWriter<ProgramState> timeWriter_;
-#if MY_DIM == 3
-  PatchInfoWriter<ProgramState, VertexBasis, GridView> patchInfoWriter_;
-#endif
-  SurfaceWriter<ProgramState, GridView> surfaceWriter_;
-  FrictionalBoundaryWriter<ProgramState, GridView> frictionalBoundaryWriter_;
-};
-#endif
--- a/src/hdf5/patchinfo-writer_tmpl.cc
+++ b/src/hdf5/patchinfo-writer_tmpl.cc
-#include "../explicitvectors.hh"
-#include "../explicitgrid.hh"
-
-#include "../program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-using P1Basis = P1NodalBasis<GridView, double>;
-using MyFunction = BasisGridFunction<P1Basis, Vector>;
-
-template class GridEvaluator<LocalVector, GridView>;
-template Dune::Matrix<typename MyFunction::RangeType>
-GridEvaluator<LocalVector, GridView>::evaluate(MyFunction const &f) const;
-
-template class PatchInfoWriter<MyProgramState, P1Basis, GridView>;
--- a/src/io/hdf5-bodywriter.hh
+++ b/src/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 VertexCoordinates = typename ProgramState::Vector;
+    using FrictionPatches = std::vector<const BoundaryPatch<GridView>* >;
+    using WeakPatches = std::vector<const ConvexPolyhedron<LocalVector>* >;
+
+    using LocalVector = typename VertexCoordinates::block_type;
+
+    friend class HDF5LevelWriter<ProgramState, VertexBasis, GridView>;
+
+public:
+    HDF5BodyWriter(
+            HDF5::Grouplike &file,
+            const VertexCoordinates& vertexCoordinates,
+            const VertexBasis& vertexBasis,
+            const FrictionPatches& frictionPatches,
+            const WeakPatches& weakPatches) :
+        patchCount_(frictionPatches.size()),
+        file_(file),
+#if MY_DIM == 3
+        patchInfoWriters_(patchCount_),
+#endif
+        frictionBoundaryWriters_(patchCount_) {
+
+        assert(patchCount_ == weakPatches.size());
+
+        for (size_t i=0; i<patchCount_; i++) {
+#if MY_DIM == 3
+            patchInfoWriters_[i] = new PatchInfoWriter<ProgramState, VertexBasis, GridView>(file_, vertexBasis, *frictionPatches[i], *weakPatches[i], i);
+#endif
+            frictionBoundaryWriters_[i] = new FrictionalBoundaryWriter<ProgramState, GridView>(file_, vertexCoordinates, *frictionPatches[i], i);
+        }
+    }
+
+    ~HDF5BodyWriter() {
+        for (size_t i=0; i<patchCount_; i++) {
+#if MY_DIM == 3
+            delete patchInfoWriters_[i];
+#endif
+            delete frictionBoundaryWriters_[i];
+        }
+    }
+
+    template <class Friction>
+    void reportSolution(ProgramState const &programState,
+                      // for the friction coefficient
+                      std::vector<std::shared_ptr<Friction>>& friction) {
+
+        for (size_t i=0; i<patchCount_; i++) {
+#if MY_DIM == 3
+            patchInfoWriters_[i]->write(programState);
+#endif
+            frictionBoundaryWriters_[i]->write(programState, *friction[i]);
+        }
+    }
+
+private:
+    const size_t patchCount_;
+
+    HDF5::Grouplike &file_;
+
+#if MY_DIM == 3
+    std::vector<PatchInfoWriter<ProgramState, VertexBasis, GridView>* > patchInfoWriters_;
+#endif
+
+    std::vector<FrictionalBoundaryWriter<ProgramState, GridView>* > frictionBoundaryWriters_;
+};
+#endif
--- a/src/io/hdf5-levelwriter.hh
+++ b/src/io/hdf5-levelwriter.hh
+#ifndef SRC_HDF5_LEVELWRITER_HH
+#define SRC_HDF5_LEVELWRITER_HH
+
+#include <dune/fufem/hdf5/file.hh>
+
+#include "hdf5/iteration-writer.hh"
+#include "hdf5/time-writer.hh"
+#include "hdf5-bodywriter.hh"
+
+template <class ProgramState, class VertexBasis, class GridView>
+class HDF5LevelWriter {
+private:
+    using HDF5BodyWriter = HDF5BodyWriter<ProgramState, VertexBasis, GridView>;
+    using VertexCoordinates = std::vector<typename HDF5BodyWriter::VertexCoordinates>;
+    using VertexBases = std::vector<VertexBasis>;
+    using FrictionPatches = std::vector<typename HDF5BodyWriter::FrictionPatches>;
+    using WeakPatches = std::vector<typename HDF5BodyWriter::WeakPatches>;
+
+    friend class HDF5NetworkWriter<ProgramState, VertexBasis, GridView>;
+
+public:
+
+    HDF5LevelWriter(HDF5::Grouplike &file,
+             const VertexCoordinates& vertexCoordinates,
+             const VertexBases& vertexBases,
+             const FrictionPatches& frictionBoundaries,
+             const WeakPatches& weakPatches)
+      : bodyCount_(vertexCoordinates.size()),
+        file_(file),
+        iterationWriter_(file_),
+        timeWriter_(file_),
+        bodyWriters_(bodyCount_) {
+
+        for (size_t i=0; i<bodyCount_; i++) {
+            bodyWriters_[i] = new HDF5BodyWriter<ProgramState, VertexBasis, GridView>(file_, vertexCoordinates[i], vertexBases[i], frictionBoundaries[i], weakPatches[i]);
+        }
+    }
+
+    ~HDF5LevelWriter() {
+        for (size_t i=0; i<bodyCount_; i++) {
+            delete bodyWriters_[i];
+        }
+    }
+
+    template <class Friction>
+    void reportSolution(
+            ProgramState const &programState,
+            // for the friction coefficient
+            std::vector<std::shared_ptr<Friction>>& friction) {
+        timeWriter_.write(programState);
+
+        for (size_t i=0; i<bodyCount_; i++) {
+            bodyWriters_[i]->reportSolution(programState, *friction[i]); //TODO
+        }
+    }
+
+    void reportIterations(
+            ProgramState const &programState,
+            IterationRegister const &iterationCount) {
+        iterationWriter_.write(programState.timeStep, iterationCount);
+  }
+
+private:
+  const size_t bodyCount_;
+
+  HDF5::Grouplike &file_;
+
+  IterationWriter iterationWriter_;
+  TimeWriter<ProgramState> timeWriter_;
+
+  std::vector<HDF5BodyWriter* > bodyWriters_;
+};
+#endif
--- a/src/hdf5/frictionalboundary-writer.cc
+++ b/src/hdf5/frictionalboundary-writer.cc
@@ -10,8 +10,9 @@
 template <class ProgramState, class GridView>
 FrictionalBoundaryWriter<ProgramState, GridView>::FrictionalBoundaryWriter(
    HDF5::Grouplike &file, Vector const &vertexCoordinates,
-    Patch const &frictionalBoundary)
-    : group_(file, "frictionalBoundary"),
+    Patch const &frictionalBoundary, size_t const id)
+    : id_(id),
+      group_(file, "frictionalBoundary"+std::to_string(id_)),
      frictionalBoundary_(frictionalBoundary),
      frictionalBoundaryDisplacementWriter_(group_, "displacement",
                                            frictionalBoundary.numVertices(),
@@ -35,23 +36,24 @@ template <class ProgramState, class GridView>
 template <class Friction>
 void FrictionalBoundaryWriter<ProgramState, GridView>::write(
    ProgramState const &programState, Friction &friction) {
+
  auto const frictionalBoundaryDisplacements =
-      restrictToSurface(programState.u, frictionalBoundary_);
+      restrictToSurface(programState.u[id_], frictionalBoundary_);
  addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep,
           frictionalBoundaryDisplacements);

  auto const frictionalBoundaryVelocities =
-      restrictToSurface(programState.v, frictionalBoundary_);
+      restrictToSurface(programState.v[id_], frictionalBoundary_);
  addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep,
           frictionalBoundaryVelocities);

  auto const frictionalBoundaryStates =
-      restrictToSurface(programState.alpha, frictionalBoundary_);
+      restrictToSurface(programState.alpha[id_], frictionalBoundary_);
  addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
           frictionalBoundaryStates);

  friction.updateAlpha(programState.alpha);
-  auto const c = friction.coefficientOfFriction(programState.v);
+  auto const c = friction.coefficientOfFriction(programState.v[id_]);
  auto const frictionalBoundaryCoefficient =
      restrictToSurface(c, frictionalBoundary_);
  addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep,

--- a/src/hdf5/frictionalboundary-writer.hh
+++ b/src/hdf5/frictionalboundary-writer.hh
@@ -11,13 +11,15 @@ template <class ProgramState, class GridView> class FrictionalBoundaryWriter {
  using Patch = BoundaryPatch<GridView>;

 public:
-  FrictionalBoundaryWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
-                           Patch const &frictionalBoundary);
+  FrictionalBoundaryWriter(HDF5::Grouplike& file, Vector const &vertexCoordinates,
+                           Patch const &frictionalBoundary, size_t const id);

  template <class Friction>
  void write(ProgramState const &programState, Friction &friction);

 private:
+  size_t const id_;
+
  HDF5::Group group_;

  Patch const &frictionalBoundary_;

--- a/src/hdf5/frictionalboundary-writer_tmpl.cc
+++ b/src/hdf5/frictionalboundary-writer_tmpl.cc
-#include "../explicitvectors.hh"
-#include "../explicitgrid.hh"
+#include "../../explicitvectors.hh"
+#include "../../explicitgrid.hh"

-#include "../program_state.hh"
+#include "../../data-structures/program_state.hh"

 using MyProgramState = ProgramState<Vector, ScalarVector>;
 using MyFriction = GlobalFriction<Matrix, Vector>;

-template class FrictionalBoundaryWriter<MyProgramState, GridView>;
-template void FrictionalBoundaryWriter<MyProgramState, GridView>::write(
+template class FrictionalBoundaryWriter<MyProgramState, DeformedGrid::LeafGridView>;
+template void FrictionalBoundaryWriter<MyProgramState, DeformedGrid::LeafGridView>::write(
    MyProgramState const &programState, MyFriction &friction);
+
+
--- a/src/hdf5/iteration-writer.cc
+++ b/src/hdf5/iteration-writer.cc
--- a/src/hdf5/iteration-writer.hh
+++ b/src/hdf5/iteration-writer.hh
@@ -3,7 +3,7 @@

 #include <dune/fufem/hdf5/sequenceio.hh>

-#include "../time-stepping/adaptivetimestepper.hh"
+#include "../../time-stepping/adaptivetimestepper.hh"

 class IterationWriter {
 public:

--- a/src/hdf5/patchinfo-writer.cc
+++ b/src/hdf5/patchinfo-writer.cc
@@ -68,8 +68,10 @@ template <class ProgramState, class VertexBasis, class GridView>
 PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter(
    HDF5::Grouplike &file, VertexBasis const &vertexBasis,
    Patch const &frictionalBoundary,
-    ConvexPolyhedron<LocalVector> const &weakPatch)
-    : group_(file, "weakPatchGrid"),
+    ConvexPolyhedron<LocalVector> const &weakPatch,
+    size_t const id)
+    : id_(id),
+      group_(file, "weakPatchGrid"+std::to_string(id_)),
      vertexBasis_(vertexBasis),
      gridEvaluator_(weakPatch, frictionalBoundary.gridView()),
      weakPatchGridVelocityWriter_(
@@ -87,7 +89,8 @@ PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter(
 template <class ProgramState, class VertexBasis, class GridView>
 void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write(
    ProgramState const &programState) {
-  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v);
+    // TODO
+  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v[id_]);
  auto const gridVelocity = gridEvaluator_.evaluate(velocity);
  addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
 }

--- a/src/hdf5/patchinfo-writer.hh
+++ b/src/hdf5/patchinfo-writer.hh
@@ -8,7 +8,7 @@
 #include <dune/fufem/geometry/convexpolyhedron.hh>
 #include <dune/fufem/hdf5/sequenceio.hh>

-#include "../one-body-problem-data/mygeometry.hh"
+#include "../../one-body-problem-data/mygeometry.hh"

 template <class LocalVector, class GridView> class GridEvaluator {
  using Element = typename GridView::Grid::template Codim<0>::Entity;
@@ -36,11 +36,14 @@ class PatchInfoWriter {
 public:
  PatchInfoWriter(HDF5::Grouplike &file, VertexBasis const &vertexBasis,
                  Patch const &frictionalBoundary,
-                  ConvexPolyhedron<LocalVector> const &weakPatch);
+                  ConvexPolyhedron<LocalVector> const &weakPatch,
+                  size_t const id);

  void write(ProgramState const &programState);

 private:
+  const size_t id_;
+
  HDF5::Group group_;

  VertexBasis const &vertexBasis_;
No results found