dune-tectonic.pc.in

@@ -10,6 +10,6 @@ Name: @PACKAGE_NAME@
 Version: @VERSION@
 Description: dune-tectonic module
 URL: http://dune-project.org/
-Requires: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Requires: dune-common dune-contact dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
 Libs: -L${libdir}
 Cflags: -I${includedir}

dune.module

@@ -5,4 +5,4 @@
 Module: dune-tectonic
 Version: 2.5-1
 Maintainer: elias.pipping@fu-berlin.de
-Depends: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Depends: dune-common dune-contact dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid

dune/tectonic/CMakeLists.txt

 install(FILES
-  body.hh
+  bodydata.hh
   frictiondata.hh
   frictionpotential.hh
   globalfrictiondata.hh
@@ -12,4 +12,5 @@ install(FILES
   mydirectionalconvexfunction.hh
   quadraticenergy.hh
   tectonic.hh
+  transformedglobalratestatefriction.hh
   DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/body.hh → dune/tectonic/bodydata.hh

-#ifndef DUNE_TECTONIC_BODY_HH
-#define DUNE_TECTONIC_BODY_HH
+#ifndef DUNE_TECTONIC_BODY_DATA_HH
+#define DUNE_TECTONIC_BODY_DATA_HH
 
-template <int dimension> struct Body {
+template <int dimension> struct BodyData {
   using ScalarFunction =
       Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
                             Dune::FieldVector<double, 1>>;

dune/tectonic/globalfriction.hh

@@ -9,9 +9,12 @@
 
 #include <dune/solvers/common/interval.hh>
 
+#include <dune/tnnmg/functionals/nonsmoothconvexfunctional.hh>
+
 #include <dune/tectonic/localfriction.hh>
 
-template <class Matrix, class Vector> class GlobalFriction {
+template <class Matrix, class Vector>
+class GlobalFriction { //: public Dune::TNNMG::NonsmoothConvexFunctional<> {
 protected:
   using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
 
@@ -82,5 +85,6 @@ template <class Matrix, class Vector> class GlobalFriction {
   }
 
   void virtual updateAlpha(ScalarVector const &alpha) = 0;
+
 };
 #endif

dune/tectonic/globalratestatefriction.hh

@@ -32,7 +32,7 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
       : restrictions(), localToGlobal(), zeroFriction() {
     auto const gridView = frictionalBoundary.gridView();
     Dune::MultipleCodimMultipleGeomTypeMapper<
-        GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
+        GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
     for (auto it = gridView.template begin<block_size>();
          it != gridView.template end<block_size>(); ++it) {
       auto const i = vertexMapper.index(*it);
@@ -44,8 +44,8 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
       restrictions.emplace_back(weights[i], weightedNormalStress[i],
                                 frictionInfo(geoToPoint(it->geometry())));
     }
-    assert(restrictions.size() == frictionalBoundary.numVertices());
-    assert(localToGlobal.size() == frictionalBoundary.numVertices());
+    assert(restrictions.size() == (size_t) frictionalBoundary.numVertices());
+    assert(localToGlobal.size() == (size_t) frictionalBoundary.numVertices());
   }
 
   void updateAlpha(ScalarVector const &alpha) override {

dune/tectonic/myblockproblem.hh

@@ -57,7 +57,7 @@ class MyBlockProblem : /* not public */ BlockNonlinearGSProblem<ConvexProblem> {
   MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem)
       : Base(parset, problem),
         maxEigenvalues_(problem.f.size()),
-        localBisection(0.0, 1.0, 0.0, true, 0.0) {
+        localBisection(0.0, 1.0, 0.0, 0.0) {
     for (size_t i = 0; i < problem.f.size(); ++i) {
       LocalVectorType eigenvalues;
       Dune::FMatrixHelp::eigenValues(problem.A[i][i], eigenvalues);

dune/tectonic/mydirectionalconvexfunction.hh

-#ifndef DUNE_TECTONIC_MYDIRECTIONALCONVEXFUNCTION_HH
-#define DUNE_TECTONIC_MYDIRECTIONALCONVEXFUNCTION_HH
-
-// Copied from dune/tnnmg/problem-classes/directionalconvexfunction.hh
-// Allows phi to be const
-
-#include <dune/fufem/arithmetic.hh>
-#include <dune/solvers/common/interval.hh>
-
-template <class Nonlinearity> class MyDirectionalConvexFunction {
-public:
-  using Vector = typename Nonlinearity::VectorType;
-  using Matrix = typename Nonlinearity::MatrixType;
-
-  MyDirectionalConvexFunction(double A, double b, Nonlinearity const &phi,
-                              Vector const &u, Vector const &v)
-      : A(A), b(b), phi(phi), u(u), v(v) {
-    phi.directionalDomain(u, v, dom);
-  }
-
-  double quadraticPart() const { return A; }
-
-  double linearPart() const { return b; }
-
-  void subDiff(double x, Dune::Solvers::Interval<double> &D) const {
-    Vector uxv = u;
-    Arithmetic::addProduct(uxv, x, v);
-    phi.directionalSubDiff(uxv, v, D);
-    auto const Axmb = A * x - b;
-    D[0] += Axmb;
-    D[1] += Axmb;
-  }
-
-  void domain(Dune::Solvers::Interval<double> &domain) const {
-    domain[0] = this->dom[0];
-    domain[1] = this->dom[1];
-  }
-
-  double A;
-  double b;
-
-private:
-  Nonlinearity const &phi;
-  Vector const &u;
-  Vector const &v;
-
-  Dune::Solvers::Interval<double> dom;
-};
-
-/*
-  1/2 <A(u + hv),u + hv> - <b, u + hv>
-  = 1/2 <Av,v> h^2 - <b - Au, v> h + const.
-
-  localA = <Av,v>
-  localb = <b - Au, v>
-*/
-
-template <class Matrix, class Vector, class Nonlinearity>
-MyDirectionalConvexFunction<Nonlinearity> restrict(Matrix const &A,
-                                                   Vector const &b,
-                                                   Vector const &u,
-                                                   Vector const &v,
-                                                   Nonlinearity const &phi) {
-  return MyDirectionalConvexFunction<Nonlinearity>(
-      Arithmetic::Axy(A, v, v), Arithmetic::bmAxy(A, b, u, v), phi, u, v);
-}
-
-#endif

dune/tectonic/tectonic.hh

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

dune/tectonic/transformedglobalratestatefriction.hh

+#ifndef DUNE_TECTONIC_TRANSFORMED_GLOBALRATESTATEFRICTION_HH
+#define DUNE_TECTONIC_TRANSFORMED_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/contact/assemblers/nbodyassembler.hh>
+
+//#include <dune/tectonic/geocoordinate.hh>
+//#include <dune/tectonic/globalfrictiondata.hh>
+#include <dune/tectonic/globalratestatefriction.hh>
+//#include <dune/tectonic/index-in-sorted-range.hh>
+
+template <class Matrix, class Vector, class ScalarFriction, class GridView>
+class TransformedGlobalRateStateFriction : public GlobalRateStateFriction<Matrix, Vector, ScalarFriction, GridView> {
+public:
+  using NBodyAssembler = Dune::Contact::NBodyAssembler<typename GridView::Grid, Vector>;
+
+  using GlobalFriction<Matrix, Vector>::block_size;
+  using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;
+
+private:
+  using typename GlobalFriction<Matrix, Vector>::ScalarVector;
+
+public:
+  TransformedGlobalRateStateFriction(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, Dune::mcmgVertexLayout());
+    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() == (size_t) frictionalBoundary.numVertices());
+    assert(localToGlobal.size() == (size_t) frictionalBoundary.numVertices());
+  }
+
+  double operator()(Vector const &x) const {
+    Vector nodalX(x.size());
+    nBodyAssembler_.
+    double tmp = 0;
+    for (size_t i = 0; i < x.size(); ++i) {
+      tmp += restriction(i)(x[i]);
+    }
+    return tmp;
+  }
+
+  /*
+    Return a restriction of the outer function to the i'th node.
+  */
+  LocalNonlinearity const virtual &restriction(size_t i) const = 0;
+
+  void addHessian(Vector const &v, Matrix &hessian) const {
+    for (size_t i = 0; i < v.size(); ++i)
+      restriction(i).addHessian(v[i], hessian[i][i]);
+  }
+
+  void directionalDomain(Vector const &, Vector const &,
+                         Dune::Solvers::Interval<double> &dom) const {
+    dom[0] = -std::numeric_limits<double>::max();
+    dom[1] = std::numeric_limits<double>::max();
+  }
+
+  void directionalSubDiff(
+      Vector const &u, Vector const &v,
+      Dune::Solvers::Interval<double> &subdifferential) const {
+    subdifferential[0] = subdifferential[1] = 0;
+    for (size_t i = 0; i < u.size(); ++i) {
+      Dune::Solvers::Interval<double> D;
+      restriction(i).directionalSubDiff(u[i], v[i], D);
+      subdifferential[0] += D[0];
+      subdifferential[1] += D[1];
+    }
+  }
+
+  void addHessianIndices(Dune::MatrixIndexSet &indices) const {
+    for (size_t i = 0; i < indices.rows(); ++i)
+      indices.add(i, i);
+  }
+
+  void addGradient(Vector const &v, Vector &gradient) const {
+    for (size_t i = 0; i < v.size(); ++i)
+      restriction(i).addGradient(v[i], gradient[i]);
+  }
+
+  double regularity(size_t i, typename Vector::block_type const &x) const {
+    return restriction(i).regularity(x);
+  }
+
+  ScalarVector coefficientOfFriction(Vector const &x) const {
+    ScalarVector ret(x.size());
+    for (size_t i = 0; i < x.size(); ++i)
+      ret[i] = restriction(i).coefficientOfFriction(x[i]);
+    return ret;
+  }
+
+  /*
+    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:
+  const NBodyAssembler& nBodyAssembler_;
+};
+#endif

program_structure.txt

+--------------------------
+--  Program structure:  --
+--------------------------
+
+
+1. build n-body system 
+    - contains grids, couplings, gridViews, assembler 
+    
+    Data structure: LevelContactNetwork
+    Factories:      StackedBlocksFactory
+    
+2. initialize/set up program state
+    - holds bodyStates, u, v, a, alpha for each body
+    - defines time, timeStep
+    - computes initial conditions
+    
+    Data structure: ProgramState
+
+-- tested until here
+
+3. assemble RSD friction
+
+4. set up TNNMG solver
+    - rate updater
+    - state updater
+    
+5. adaptive time stepper
+
+ 

src/CMakeLists.txt

+add_subdirectory("tests")
+
 set(SW_SOURCE_FILES
   assemblers.cc
-  enumparser.cc
-  hdf5/frictionalboundary-writer.cc
-  hdf5/iteration-writer.cc
-  hdf5/patchinfo-writer.cc
-  hdf5/restart-io.cc
-  hdf5/surface-writer.cc
-  hdf5/time-writer.cc
-  one-body-problem-data/mygeometry.cc
-  one-body-problem-data/mygrid.cc
-  one-body-problem.cc
+  data-structures/body.cc
+  data-structures/enumparser.cc
+  io/vtk.cc
+  io/hdf5/frictionalboundary-writer.cc
+  io/hdf5/iteration-writer.cc
+  io/hdf5/patchinfo-writer.cc
+  io/hdf5/restart-io.cc
+  io/hdf5/surface-writer.cc
+  io/hdf5/time-writer.cc
+#  one-body-problem-data/mygeometry.cc
+#  one-body-problem-data/mygrid.cc
+#  one-body-problem.cc
   spatial-solving/fixedpointiterator.cc
   spatial-solving/solverfactory.cc
   time-stepping/adaptivetimestepper.cc
@@ -17,26 +21,63 @@ set(SW_SOURCE_FILES
   time-stepping/rate.cc
   time-stepping/rate/rateupdater.cc
   time-stepping/state.cc
-  vtk.cc
 )
+
+set(MSW_SOURCE_FILES
+  assemblers.cc
+  nodalweights.cc
+  data-structures/body.cc
+  #data-structures/contactnetwork.cc
+  data-structures/enumparser.cc
+  data-structures/levelcontactnetwork.cc
+  #factories/cantorfactory.cc
+  factories/stackedblocksfactory.cc
+  io/vtk.cc
+  io/hdf5/frictionalboundary-writer.cc
+  io/hdf5/iteration-writer.cc
+  io/hdf5/patchinfo-writer.cc
+  io/hdf5/restart-io.cc
+  io/hdf5/surface-writer.cc
+  io/hdf5/time-writer.cc
+  #multi-body-problem-data/grid/cube.cc
+  #multi-body-problem-data/grid/cubefaces.cc
+  multi-body-problem-data/grid/cuboidgeometry.cc
+  multi-body-problem-data/grid/mygrids.cc
+  multi-body-problem-data/grid/simplexmanager.cc
+  multi-body-problem.cc
+  spatial-solving/fixedpointiterator.cc
+  spatial-solving/solverfactory.cc
+  spatial-solving/solverfactory_old.cc
+  time-stepping/adaptivetimestepper.cc
+  time-stepping/coupledtimestepper.cc
+  time-stepping/rate.cc
+  time-stepping/rate/rateupdater.cc
+  time-stepping/state.cc
+)
+
 set(UGW_SOURCE_FILES
   assemblers.cc # FIXME
+  io/uniform-grid-writer.cc
+  io/vtk.cc
   one-body-problem-data/mygrid.cc
-  uniform-grid-writer.cc
-  vtk.cc
 )
 
 foreach(_dim 2 3)
   set(_sw_target one-body-problem-${_dim}D)
+  set(_msw_target multi-body-problem-${_dim}D)
   set(_ugw_target uniform-grid-writer-${_dim}D)
 
   add_executable(${_sw_target} ${SW_SOURCE_FILES})
+  add_executable(${_msw_target} ${MSW_SOURCE_FILES})
   add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
   add_dune_ug_flags(${_sw_target})
+  add_dune_ug_flags(${_msw_target})
   add_dune_ug_flags(${_ugw_target})
 
   add_dune_hdf5_flags(${_sw_target})
+  add_dune_hdf5_flags(${_msw_target})
 
   set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  set_property(TARGET ${_msw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
   set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
 endforeach()

src/assemblers.cc

@@ -24,8 +24,8 @@
 #include "assemblers.hh"
 
 template <class GridView, int dimension>
-MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
-    : cellBasis(_gridView),
+MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView) :
+      cellBasis(_gridView),
       vertexBasis(_gridView),
       gridView(_gridView),
       cellAssembler(cellBasis, cellBasis),
@@ -33,8 +33,9 @@ MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
-    BoundaryPatch<GridView> const &frictionalBoundary,
-    ScalarMatrix &frictionalBoundaryMass) const {
+        BoundaryPatch<GridView> const &frictionalBoundary,
+        ScalarMatrix &frictionalBoundaryMass) const {
+
   BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
                         LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
       frictionalBoundaryMassAssembler(frictionalBoundary);
@@ -44,9 +45,9 @@ void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleMass(
-    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
-        densityFunction,
-    Matrix &M) const {
+        Dune::VirtualFunction<LocalVector, LocalScalarVector> const & densityFunction,
+        Matrix &M) const {
+
   // NOTE: We treat the weight as a constant function
   QuadratureRuleKey quadKey(dimension, 0);
 
@@ -58,8 +59,11 @@ void MyAssembler<GridView, dimension>::assembleMass(
 }
 
 template <class GridView, int dimension>
-void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
-                                                          Matrix &A) const {
+void MyAssembler<GridView, dimension>::assembleElasticity(
+        double E,
+        double nu,
+        Matrix &A) const {
+
   StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
       localStiffness(E, nu);
   vertexAssembler.assembleOperator(localStiffness, A);
@@ -67,9 +71,10 @@ void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleViscosity(
-    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
-    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
-    Matrix &C) const {
+        Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
+        Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
+        Matrix &C) const {
+
   // NOTE: We treat the weights as constant functions
   QuadratureRuleKey shearViscosityKey(dimension, 0);
   QuadratureRuleKey bulkViscosityKey(dimension, 0);
@@ -83,8 +88,9 @@ void MyAssembler<GridView, dimension>::assembleViscosity(
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleBodyForce(
-    Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
-    Vector &f) const {
+        Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
+        Vector &f) const {
+
   L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
       gravityFunctionalAssembler(gravityField);
   vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f);
@@ -92,9 +98,11 @@ void MyAssembler<GridView, dimension>::assembleBodyForce(
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleNeumann(
-    BoundaryPatch<GridView> const &neumannBoundary, Vector &f,
-    Dune::VirtualFunction<double, double> const &neumann,
-    double relativeTime) const {
+        BoundaryPatch<GridView> const &neumannBoundary,
+        Vector &f,
+        Dune::VirtualFunction<double, double> const &neumann,
+        double relativeTime) const {
+
   LocalVector localNeumann(0);
   neumann.evaluate(relativeTime, localNeumann[0]);
   NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
@@ -106,9 +114,11 @@ void MyAssembler<GridView, dimension>::assembleNeumann(
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
-    BoundaryPatch<GridView> const &frictionalBoundary,
-    ScalarVector &weightedNormalStress, double youngModulus,
-    double poissonRatio, Vector const &displacement) const {
+        BoundaryPatch<GridView> const &frictionalBoundary,
+        ScalarVector &weightedNormalStress,
+        double youngModulus,
+        double poissonRatio,
+        Vector const &displacement) const {
 
   BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
                                                               displacement);
@@ -138,11 +148,12 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
 
 template <class GridView, int dimension>
 auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
-    Config::FrictionModel frictionModel,
-    BoundaryPatch<GridView> const &frictionalBoundary,
-    GlobalFrictionData<dimension> const &frictionInfo,
-    ScalarVector const &weightedNormalStress) const
-    -> std::shared_ptr<GlobalFriction<Matrix, Vector>> {
+        Config::FrictionModel frictionModel,
+        BoundaryPatch<GridView> const &frictionalBoundary,
+        GlobalFrictionData<dimension> const &frictionInfo,
+        ScalarVector const &weightedNormalStress) const
+-> std::shared_ptr<GlobalFriction<Matrix, Vector>> {
+
   // Lumping of the nonlinearity
   ScalarVector weights;
   {
@@ -169,8 +180,11 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleVonMisesStress(
-    double youngModulus, double poissonRatio, Vector const &u,
-    ScalarVector &stress) const {
+        double youngModulus,
+        double poissonRatio,
+        Vector const &u,
+        ScalarVector &stress) const {
+
   auto const gridDisplacement =
       std::make_shared<BasisGridFunction<VertexBasis, Vector> const>(
           vertexBasis, u);

src/assemblers.hh

@@ -16,75 +16,85 @@
 #include <dune/tectonic/globalfriction.hh>
 #include <dune/tectonic/globalfrictiondata.hh>
 
-#include "enums.hh"
+#include "data-structures/enums.hh"
 
-template <class GridView, int dimension> class MyAssembler {
+template <class GridView, int dimension>
+class MyAssembler {
 public:
-  using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
-  using Matrix =
+    using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using Matrix =
       Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>;
-  using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
-  using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
+    using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+    using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
 
-  using CellBasis = P0Basis<GridView, double>;
-  using VertexBasis = P1NodalBasis<GridView, double>;
+    using CellBasis = P0Basis<GridView, double>;
+    using VertexBasis = P1NodalBasis<GridView, double>;
 
-  CellBasis const cellBasis;
-  VertexBasis const vertexBasis;
-  GridView const &gridView;
+    CellBasis const cellBasis;
+    VertexBasis const vertexBasis;
+    GridView const &gridView;
 
 private:
-  using Grid = typename GridView::Grid;
-  using LocalVector = typename Vector::block_type;
-  using LocalScalarVector = typename ScalarVector::block_type;
+    using Grid = typename GridView::Grid;
+    using LocalVector = typename Vector::block_type;
+    using LocalScalarVector = typename ScalarVector::block_type;
 
-  using LocalCellBasis = typename CellBasis::LocalFiniteElement;
-  using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
+    using LocalCellBasis = typename CellBasis::LocalFiniteElement;
+    using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
 
-  Assembler<CellBasis, CellBasis> cellAssembler;
-  Assembler<VertexBasis, VertexBasis> vertexAssembler;
+    Assembler<CellBasis, CellBasis> cellAssembler;
+    Assembler<VertexBasis, VertexBasis> vertexAssembler;
 
 public:
-  MyAssembler(GridView const &gridView);
-
-  void assembleFrictionalBoundaryMass(
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      ScalarMatrix &frictionalBoundaryMass) const;
-
-  void assembleMass(Dune::VirtualFunction<
-                        LocalVector, LocalScalarVector> const &densityFunction,
-                    Matrix &M) const;
-
-  void assembleElasticity(double E, double nu, Matrix &A) const;
-
-  void assembleViscosity(
-      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
-          shearViscosity,
-      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
-          bulkViscosity,
-      Matrix &C) const;
-
-  void assembleBodyForce(
-      Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
-      Vector &f) const;
-
-  void assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
-                       Vector &f,
-                       Dune::VirtualFunction<double, double> const &neumann,
-                       double relativeTime) const;
-
-  void assembleWeightedNormalStress(
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      ScalarVector &weightedNormalStress, double youngModulus,
-      double poissonRatio, Vector const &displacement) const;
-
-  std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity(
-      Config::FrictionModel frictionModel,
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      GlobalFrictionData<dimension> const &frictionInfo,
-      ScalarVector const &weightedNormalStress) const;
-
-  void assembleVonMisesStress(double youngModulus, double poissonRatio,
-                              Vector const &u, ScalarVector &stress) const;
+    MyAssembler(GridView const &gridView);
+
+    void assembleFrictionalBoundaryMass(
+            BoundaryPatch<GridView> const &frictionalBoundary,
+            ScalarMatrix &frictionalBoundaryMass) const;
+
+    void assembleMass(
+            Dune::VirtualFunction<
+            LocalVector, LocalScalarVector> const &densityFunction,
+            Matrix &M) const;
+
+    void assembleElasticity(
+            double E,
+            double nu,
+            Matrix &A) const;
+
+    void assembleViscosity(
+            Dune::VirtualFunction<LocalVector, LocalScalarVector> const & shearViscosity,
+            Dune::VirtualFunction<LocalVector, LocalScalarVector> const & bulkViscosity,
+            Matrix &C) const;
+
+    void assembleBodyForce(
+            Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
+            Vector &f) const;
+
+    void assembleNeumann(
+            BoundaryPatch<GridView> const &neumannBoundary,
+            Vector &f,
+            Dune::VirtualFunction<double, double> const &neumann,
+            double relativeTime) const;
+
+    void assembleWeightedNormalStress(
+            BoundaryPatch<GridView> const &frictionalBoundary,
+            ScalarVector &weightedNormalStress,
+            double youngModulus,
+            double poissonRatio,
+            Vector const &displacement) const;
+
+    auto assembleFrictionNonlinearity(
+            Config::FrictionModel frictionModel,
+            BoundaryPatch<GridView> const &frictionalBoundary,
+            GlobalFrictionData<dimension> const &frictionInfo,
+            ScalarVector const &weightedNormalStress) const
+    -> std::shared_ptr<GlobalFriction<Matrix, Vector>>;
+
+    void assembleVonMisesStress(
+            double youngModulus,
+            double poissonRatio,
+            Vector const &u,
+            ScalarVector &stress) const;
 };
 #endif

src/assemblers_tmpl.cc

@@ -3,5 +3,8 @@
 #endif
 
 #include "explicitgrid.hh"
+#include "explicitvectors.hh"
 
-template class MyAssembler<GridView, MY_DIM>;
+#include "data-structures/body.hh"
+
+template class MyAssembler<Body<Grid, Vector, MY_DIM>::DeformedLeafGridView, MY_DIM>;

src/boundarycondition.hh

+#ifndef SRC_BOUNDARYCONDITION_HH
+#define SRC_BOUNDARYCONDITION_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+
+template <class GridView, int dims>
+class BoundaryCondition {
+public:
+  using BoundaryPatch = BoundaryPatch<GridView>;
+
+private:
+  using Function = Dune::VirtualFunction<double, double>;
+
+  const std::string tag_; // store type of boundary condition, e.g. dirichlet, neumann, friction, etc
+
+  std::shared_ptr<BoundaryPatch> boundaryPatch_;
+  std::shared_ptr<Dune::BitSetVector<dims>> boundaryNodes_;
+  std::shared_ptr<Function> boundaryFunction_;
+
+public:
+  BoundaryCondition(const std::string& tag = "") :
+      tag_(tag)
+  {}
+
+  BoundaryCondition(std::shared_ptr<BoundaryPatch> patch, std::shared_ptr<Function> function, const std::string& tag = "") :
+      tag_(tag),
+      boundaryPatch_(patch),
+      boundaryFunction_(function)
+  {}
+
+  void setBoundaryPatch(const GridView& gridView, std::shared_ptr<Dune::BitSetVector<dims>> nodes) {
+      boundaryNodes_ = nodes;
+      boundaryPatch_ = std::make_shared<BoundaryPatch>(gridView, *nodes);
+  }
+
+  void setBoundaryPatch(std::shared_ptr<BoundaryPatch> patch) {
+      boundaryPatch_ = patch;
+  }
+
+  void setBoundaryPatch(const BoundaryPatch& patch) {
+      boundaryPatch_ = std::make_shared<BoundaryPatch>(patch);
+  }
+
+  void setBoundaryFunction(std::shared_ptr<Function> function) {
+      boundaryFunction_ = function;
+  }
+
+  void set(std::shared_ptr<BoundaryPatch> patch, std::shared_ptr<Function> function) {
+      boundaryPatch_ = patch;
+      boundaryFunction_ = function;
+  }
+
+  const std::string& tag() const {
+      return tag_;
+  }
+
+  const std::shared_ptr<BoundaryPatch>& boundaryPatch() const {
+      return boundaryPatch_;
+  }
+
+  const std::shared_ptr<Dune::BitSetVector<dims>>& boundaryNodes() const {
+      return boundaryNodes_;
+  }
+
+  const std::shared_ptr<Function>& boundaryFunction() const {
+      return boundaryFunction_;
+  }
+};
+
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include <dune/tectonic/globalfrictiondata.hh>
+
+template <class GridView, class LocalVectorType, int dims>
+class FrictionBoundaryCondition : public BoundaryCondition<GridView, dims>{
+private:
+    using Base = BoundaryCondition<GridView, dims>;
+    using LocalVector = LocalVectorType;
+
+    // friction data
+    std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch_;
+    std::shared_ptr<GlobalFrictionData<dims>> frictionData_;
+
+public:
+  FrictionBoundaryCondition() :
+      Base("friction")
+  {}
+
+  FrictionBoundaryCondition(std::shared_ptr<typename Base::BoundaryPatch> patch, std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch, std::shared_ptr<GlobalFrictionData<dims>> frictionData) :
+      Base(patch, nullptr, "friction"),
+      weakeningPatch_(weakeningPatch),
+      frictionData_(frictionData)
+  {}
+
+  void setWeakeningPatch(std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch) {
+      weakeningPatch_ = weakeningPatch;
+  }
+
+  void setFrictionData(std::shared_ptr<GlobalFrictionData<dims>> frictionData) {
+      frictionData_ = frictionData;
+  }
+
+  const std::shared_ptr<ConvexPolyhedron<LocalVector>>& weakeningPatch() const {
+      return weakeningPatch_;
+  }
+
+  const std::shared_ptr<GlobalFrictionData<dims>>& frictionData() const {
+      return frictionData_;
+  }
+};
+#endif

src/data-structures/.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+
+template class ContactNetwork<Grid, MY_DIM>;

src/data-structures/body.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+
+#include "body.hh"
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+Body<GridTEMPLATE, VectorTEMPLATE, dims>::Body(
+        const std::shared_ptr<BodyData<dims>>& bodyData,
+        const std::shared_ptr<GridType>& grid) :
+    bodyData_(bodyData),
+    grid_(grid),
+    deformation_(std::make_shared<DeformationFunction>(grid_->leafGridView())),
+    deformedGrid_(std::make_shared<DeformedGridType>(*grid_, *deformation_)),
+    leafView_(deformedGrid_->leafGridView()),
+    leafVertexCount_(leafView_.size(dims)),
+    assembler_(std::make_shared<Assembler>(leafView_)),
+    matrices_()
+{
+    leafBoundaryConditions_.resize(0);
+
+    levelBoundaryConditions_.resize(grid_->maxLevel()+1);
+    for (size_t i=0; i<levelBoundaryConditions_.size(); i++) {
+        levelBoundaryConditions_[i].resize(0);
+    }
+
+    externalForce_ = [&](double relativeTime, VectorType &ell) {
+        // Problem formulation: right-hand side
+        std::vector<std::shared_ptr<LeafBoundaryCondition>> leafNeumannConditions;
+        leafBoundaryConditions("neumann", leafNeumannConditions);
+        ell.resize(gravityFunctional_.size());
+        ell = gravityFunctional_;
+
+        for (size_t i=0; i<leafNeumannConditions.size(); i++) {
+            const auto& leafNeumannCondition = leafNeumannConditions[i];
+
+            VectorType b;
+            assembler_->assembleNeumann(*leafNeumannCondition->boundaryPatch(), b, *leafNeumannCondition->boundaryFunction(),
+                                    relativeTime);
+            ell += b;
+        }
+    };
+}
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::assemble() {
+
+    // assemble matrices
+    assembler_->assembleElasticity(bodyData_->getYoungModulus(), bodyData_->getPoissonRatio(), *matrices_.elasticity);
+    assembler_->assembleViscosity(bodyData_->getShearViscosityField(), bodyData_->getBulkViscosityField(), *matrices_.damping);
+    assembler_->assembleMass(bodyData_->getDensityField(), *matrices_.mass);
+
+    // assemble state energy norm
+    typename FrictionBoundaryCondition::BoundaryPatch frictionBoundary(this->leafView(), false); // boundary patch containing all friction boundaries (possibly disconnected)
+    for (size_t i=0; i<frictionBoundaryConditions_.size(); i++) {
+        const auto& frictionCondition = frictionBoundaryConditions_[i];
+        //std::cout << "grid Pointer " << frictionCondition->boundaryPatch()-> << std::endl;
+        frictionBoundary.addPatch(*frictionCondition->boundaryPatch());
+    }
+
+    ScalarMatrix relativeFrictionBoundaryMass;
+    assembler_->assembleFrictionalBoundaryMass(frictionBoundary, relativeFrictionBoundaryMass);
+    relativeFrictionBoundaryMass /= frictionBoundary.area(); // TODO: weight by individual friction patches?
+    stateEnergyNorm_ = std::make_shared<const EnergyNorm<ScalarMatrix, ScalarVector>>(relativeFrictionBoundaryMass);
+
+    // assemble forces
+    assembler_->assembleBodyForce(bodyData_->getGravityField(), gravityFunctional_);
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::assembleFriction(
+        const Config::FrictionModel& frictionModel,
+        const ScalarVector& weightedNormalStress) {
+
+    globalFriction_.resize({frictionBoundaryConditions_.size()});
+
+    for (size_t i=0; i<globalFriction_.size(); i++) {
+        const auto& frictionBoundaryCondition = frictionBoundaryConditions_[i];
+        globalFriction_[i] = assembler_->assembleFrictionNonlinearity(frictionModel, *frictionBoundaryCondition->boundaryPatch(), *frictionBoundaryCondition->frictionData(), weightedNormalStress);
+    }
+}
+
+// setter and getter
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::data() const
+-> const std::shared_ptr<BodyData<dims>>& {
+
+    return bodyData_;
+}
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::setDeformation(const VectorType& def) {
+    deformation_->setDeformation(def);
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::deformation() const
+-> DeformationFunction& {
+
+    return *deformation_;
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::deformedGrid() const
+-> const DeformedGridType* {
+
+    return deformedGrid_.get();
+}
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::leafView() const
+-> const DeformedLeafGridView& {
+
+    return leafView_;
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::leafVertexCount() const
+-> int {
+
+    return leafVertexCount_;
+}
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::assembler() const
+-> const std::shared_ptr<Assembler>& {
+
+    return assembler_;
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::matrices() const
+-> const Matrices<Matrix,1>& {
+
+    return matrices_;
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::externalForce() const
+-> const ExternalForce& {
+
+    return externalForce_;
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::stateEnergyNorm() const
+-> const std::shared_ptr<const StateEnergyNorm>& {
+
+    return stateEnergyNorm_;
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::addBoundaryCondition(std::shared_ptr<LeafBoundaryCondition> leafBoundaryCondition) {
+    leafBoundaryConditions_.emplace_back(leafBoundaryCondition);
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::leafBoundaryConditions(
+        const std::string& tag,
+        std::vector<std::shared_ptr<LeafBoundaryCondition>>& selectedConditions) const {
+
+    if (tag == "friction") {
+        selectedConditions.resize(frictionBoundaryConditions_.size());
+        for (size_t i=0; i<frictionBoundaryConditions_.size(); i++) {
+                    selectedConditions[i] = frictionBoundaryConditions_[i];
+        }
+
+    } else {
+        selectedConditions.resize(0);
+        for (size_t i=0; i<leafBoundaryConditions_.size(); i++) {
+            if (leafBoundaryConditions_[i]->tag() == tag)
+                selectedConditions.emplace_back(leafBoundaryConditions_[i]);
+        }
+    }
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::leafBoundaryConditions() const
+-> const std::vector<std::shared_ptr<LeafBoundaryCondition>>& {
+
+    return leafBoundaryConditions_;
+}
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::addBoundaryCondition(
+        std::shared_ptr<LevelBoundaryCondition> levelBoundaryCondition,
+        int level) {
+
+    levelBoundaryConditions_[level].push_back(levelBoundaryCondition);
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::levelBoundaryConditions(
+        const std::string& tag,
+        std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>& selectedConditions) const {
+
+    selectedConditions.resize(levelBoundaryConditions_.size());
+    for (size_t level=0; level<levelBoundaryConditions_.size(); level++) {
+        const std::vector<std::shared_ptr<LevelBoundaryCondition>>& levelBoundaryConditions = levelBoundaryConditions_[level];
+        for (size_t i=0; i<levelBoundaryConditions.size(); i++) {
+            if (levelBoundaryConditions[i]->tag() == tag)
+                selectedConditions[level].push_back(levelBoundaryConditions[i]);
+        }
+    }
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::levelBoundaryConditions(
+        const std::string& tag,
+        int level,
+        std::vector<std::shared_ptr<LevelBoundaryCondition>>& selectedConditions) const {
+
+    selectedConditions.resize(0);
+    const std::vector<std::shared_ptr<LevelBoundaryCondition>>& levelBoundaryConditions = levelBoundaryConditions_[level];
+    for (size_t i=0; i<levelBoundaryConditions.size(); i++) {
+        if (levelBoundaryConditions[i]->tag() == tag)
+            selectedConditions.push_back(levelBoundaryConditions[i]);
+    }
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::levelBoundaryConditions(int level) const
+-> const std::vector<std::shared_ptr<LevelBoundaryCondition>>& {
+
+    return levelBoundaryConditions_[level];
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::levelBoundaryConditions() const
+-> const std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>& {
+
+    return levelBoundaryConditions_;
+}
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::addBoundaryCondition(std::shared_ptr<FrictionBoundaryCondition> frictionBoundaryCondition) {
+    frictionBoundaryConditions_.emplace_back(frictionBoundaryCondition);
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::frictionBoundaryConditions() const
+-> const std::vector<std::shared_ptr<FrictionBoundaryCondition>>& {
+
+    return frictionBoundaryConditions_;
+}
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+auto Body<GridTEMPLATE, VectorTEMPLATE, dims>::globalFriction()
+-> GlobalFrictionContainer& {
+    return globalFriction_;
+}
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::boundaryPatchNodes(
+        const std::string& tag,
+        BoundaryPatchNodes& nodes) const {
+
+    std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
+    this->leafBoundaryConditions(tag, boundaryConditions);
+
+    nodes.resize(boundaryConditions.size());
+
+    for (size_t i=0; i<nodes.size(); i++) {
+        nodes[i] = boundaryConditions[i]->boundaryPatch()->getVertices();
+    }
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::boundaryNodes(
+        const std::string& tag,
+        BoundaryNodes& nodes) const {
+
+    std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
+    this->leafBoundaryConditions(tag, boundaryConditions);
+
+    nodes.resize(boundaryConditions.size());
+
+    for (size_t i=0; i<nodes.size(); i++) {
+        nodes[i] = boundaryConditions[i]->boundaryNodes().get();
+    }
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::boundaryFunctions(
+        const std::string& tag,
+        BoundaryFunctions& functions) const {
+
+    std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
+    this->leafBoundaryConditions(tag, boundaryConditions);
+
+    functions.resize(boundaryConditions.size());
+
+    for (size_t i=0; i<functions.size(); i++) {
+        functions[i] = boundaryConditions[i]->boundaryFunction().get();
+    }
+}
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+void Body<GridTEMPLATE, VectorTEMPLATE, dims>::boundaryPatches(
+        const std::string& tag,
+        BoundaryPatches& patches) const {
+
+    std::vector<std::shared_ptr<typename Body::LeafBoundaryCondition>> boundaryConditions;
+    this->leafBoundaryConditions(tag, boundaryConditions);
+
+    patches.resize(boundaryConditions.size());
+
+    for (size_t i=0; i<patches.size(); i++) {
+        patches[i] = boundaryConditions[i]->boundaryPatch().get();
+    }
+}
+
+#include "body_tmpl.cc"

src/data-structures/body.hh

+#ifndef SRC_DATA_STRUCTURES_BODY_HH
+#define SRC_DATA_STRUCTURES_BODY_HH
+
+#include <functional>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+//#include <dune/fufem/assemblers/assembler.hh>
+//#pragma clang diagnostic push
+//#pragma clang diagnostic ignored "-Wsign-compare"
+//#include <dune/fufem/functionspacebases/p0basis.hh>
+//#pragma clang diagnostic pop
+//#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+
+#include <dune/grid/geometrygrid/grid.hh>
+
+#include <dune/fufem/functions/deformationfunction.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+
+#include <dune/tectonic/bodydata.hh>
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "../assemblers.hh"
+#include "../boundarycondition.hh"
+#include "globalfrictioncontainer.hh"
+#include "enums.hh"
+#include "matrices.hh"
+
+
+
+template <class GridTEMPLATE, class VectorTEMPLATE, int dims>
+class Body {
+public:
+    using GridType = GridTEMPLATE;
+    using VectorType = VectorTEMPLATE;
+
+    using DeformationFunction = DeformationFunction<typename GridType::LeafGridView, VectorType>;
+    using DeformedGridType = Dune::GeometryGrid<GridType, DeformationFunction>;
+    using DeformedLeafGridView = typename DeformedGridType::LeafGridView;
+    using DeformedLevelGridView = typename DeformedGridType::LevelGridView;
+
+    using Assembler = MyAssembler<DeformedLeafGridView, dims>;
+    using Matrix = typename Assembler::Matrix;
+    using LocalVector = typename VectorType::block_type;
+
+    using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+    using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+
+    using Function = Dune::VirtualFunction<double, double>;
+    using ExternalForce = std::function<void(double, VectorType &)>;
+
+    using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dims>;
+    using FrictionBoundaryCondition = FrictionBoundaryCondition<DeformedLeafGridView, LocalVector, dims>;
+    using LevelBoundaryCondition = BoundaryCondition<DeformedLevelGridView, dims>;
+
+    using BoundaryFunctions = std::vector<const Function* >;
+    using BoundaryNodes = std::vector<const Dune::BitSetVector<dims>* >;
+    using BoundaryPatchNodes = std::vector<const Dune::BitSetVector<1>* >;
+    using BoundaryPatches = std::vector<const typename LeafBoundaryCondition::BoundaryPatch* >;
+
+    using GlobalFriction = GlobalFriction<Matrix, VectorType>;
+    using GlobalFrictionContainer = GlobalFrictionContainer<GlobalFriction, 1>;
+
+    using StateEnergyNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
+
+private:
+    std::shared_ptr<Body<GridType, VectorType, dims>> parent_;
+    std::vector<std::shared_ptr<Body<GridType, VectorType, dims>>> children_;
+
+    std::shared_ptr<BodyData<dims>> bodyData_;
+
+    std::shared_ptr<GridType> grid_;
+    std::shared_ptr<DeformationFunction> deformation_;
+    std::shared_ptr<DeformedGridType> deformedGrid_;
+
+    DeformedLeafGridView leafView_;
+    int leafVertexCount_;
+
+    std::shared_ptr<Assembler> assembler_;
+    Matrices<Matrix,1> matrices_;
+
+    ExternalForce externalForce_;
+
+    VectorType gravityFunctional_;
+    std::shared_ptr<const StateEnergyNorm> stateEnergyNorm_;
+
+    // boundary conditions
+    std::vector<std::shared_ptr<LeafBoundaryCondition>> leafBoundaryConditions_;
+    std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>> levelBoundaryConditions_; // first index: level, second index: bc on lvl
+
+    // friction boundary conditions
+    std::vector<std::shared_ptr<FrictionBoundaryCondition>> frictionBoundaryConditions_;
+    GlobalFrictionContainer globalFriction_;
+
+public:
+    Body(
+            const std::shared_ptr<BodyData<dims>>& bodyData,
+            const std::shared_ptr<GridType>& grid);
+
+    void assemble();
+    void assembleFriction(
+            const Config::FrictionModel& frictionModel,
+            const ScalarVector& weightedNormalStress);
+
+    // setter and getter
+    auto data() const -> const std::shared_ptr<BodyData<dims>>&;
+
+    void setDeformation(const VectorType& def);
+    auto deformation() const -> DeformationFunction&;
+    auto deformedGrid() const -> const DeformedGridType*;
+
+    auto leafView() const -> const DeformedLeafGridView&;
+    auto leafVertexCount() const -> int;
+
+    auto assembler() const -> const std::shared_ptr<Assembler>&;
+    auto matrices() const -> const Matrices<Matrix,1>&;
+
+    auto externalForce() const -> const ExternalForce&;
+
+    auto stateEnergyNorm() const -> const std::shared_ptr<const StateEnergyNorm>&;
+
+    void addBoundaryCondition(std::shared_ptr<LeafBoundaryCondition> leafBoundaryCondition);
+    void leafBoundaryConditions(
+            const std::string& tag,
+            std::vector<std::shared_ptr<LeafBoundaryCondition>>& selectedConditions) const;
+    auto leafBoundaryConditions() const -> const std::vector<std::shared_ptr<LeafBoundaryCondition>>&;
+
+    void addBoundaryCondition(std::shared_ptr<LevelBoundaryCondition> levelBoundaryCondition, int level);
+    void levelBoundaryConditions(
+            const std::string& tag,
+            std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>& selectedConditions) const;
+    void levelBoundaryConditions(
+            const std::string& tag,
+            int level,
+            std::vector<std::shared_ptr<LevelBoundaryCondition>>& selectedConditions) const;
+    auto levelBoundaryConditions(int level) const -> const std::vector<std::shared_ptr<LevelBoundaryCondition>>&;
+    auto levelBoundaryConditions() const -> const std::vector<std::vector<std::shared_ptr<LevelBoundaryCondition>>>&;
+
+    void addBoundaryCondition(std::shared_ptr<FrictionBoundaryCondition> frictionBoundaryCondition);
+    auto frictionBoundaryConditions() const -> const std::vector<std::shared_ptr<FrictionBoundaryCondition>>&;
+
+    auto globalFriction() -> GlobalFrictionContainer&;
+
+    void boundaryPatchNodes(
+            const std::string& tag,
+            BoundaryPatchNodes& nodes) const;
+    void boundaryNodes(
+            const std::string& tag,
+            BoundaryNodes& nodes) const;
+    void boundaryFunctions(
+            const std::string& tag,
+            BoundaryFunctions& functions) const;
+    void boundaryPatches(
+            const std::string& tag,
+            BoundaryPatches& patches) const;
+};
+#endif

src/hdf5/surface-writer_tmpl.cc → src/data-structures/body_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 Body<Grid, Vector, MY_DIM>;