data/tools/support/write_contours.py

+writeContours <- function (contours, level, file) {
+  file.create(file)
+  for (cl in contours[sapply(contours, function(x) x$level) == level]) {
+    conn <- file(file, 'a')
+    write.table(cbind(cl$x, cl$y, level),
+                file = file, row.names = FALSE, col.names = FALSE,
+                append = TRUE)
+    writeLines("\n", conn)
+    close(conn)
+  }
+}

debugging.m

+% householder reflection
+
+v0 = [1; 0];
+v0 = v0./norm(v0);
+
+v1 = [1; -1];
+v1 = v1./norm(v1);
+
+v = 1/2.*(v1-v0);
+vnorm = norm(v);
+
+H = eye(2) - 2/(vnorm^2).*(v*v')
+
+A = [...
+];
+
+f = [...
+0;
+-4.33681e-19;
+0;
+1.0842e-19;
+-2.66671e-06;
+-4.94914e-07;
+0;
+-2.71051e-20;
+-4.19923e-06;
+-4.21374e-07;
+-3.92235e-06;
+-8.64817e-07;
+0;
+-5.42101e-20;
+-4.96301e-06;
+-4.13502e-06;
+0;
+0;
+-7.03515e-06;
+-5.84881e-06;
+-4.93309e-06;
+-2.37804e-06;
+0;
+0;
+-2.90755e-06;
+-1.23334e-06;
+0;
+0;
+-2.2449e-06;
+-8.21563e-07;
+0;
+0;
+0;
+0;
+0;
+0;
+-7.24652e-06;
+-8.25293e-07;
+0;
+2.4104e-06;
+-8.19294e-06;
+-1.18617e-06;
+-5.58223e-06;
+-1.69348e-06;
+0;
+1.88001e-06;
+-1.03387e-05;
+-8.22886e-06;
+0;
+0;
+-1.07258e-05;
+-8.62344e-06;
+-9.76077e-06;
+-4.93954e-06;
+-9.10261e-06;
+-7.35307e-06;
+-4.01973e-06;
+-2.38351e-06;
+-5.92659e-06;
+-4.23858e-07;
+-7.9911e-06;
+-4.04139e-06;
+-3.80594e-06;
+-3.06903e-06;
+-3.92924e-06;
+-2.71681e-06;
+-5.12975e-06;
+-1.37642e-06;
+-4.79892e-06;
+-1.14935e-06;
+-7.19726e-06;
+-3.54967e-06;
+-6.62755e-06;
+-6.1499e-07;
+0;
+1.28701e-05;
+0;
+0;
+-1.09944e-05;
+-1.92221e-06;
+0;
+1.17853e-05;
+-1.21456e-05;
+-2.20784e-06;
+-8.3682e-06;
+-2.01449e-06;
+0;
+0;
+-1.43122e-05;
+-1.00938e-05;
+0;
+0;
+-1.48374e-05;
+-1.03947e-05;
+-1.3684e-05;
+-6.31825e-06;
+-1.35318e-05;
+-9.3703e-06;
+-5.76679e-06;
+-1.88878e-06;
+-9.18524e-06;
+-1.51544e-06;
+-1.16296e-05;
+-5.64398e-06;
+-5.20301e-06;
+-2.13559e-06;
+-5.55432e-06;
+-2.017e-06;
+-7.96841e-06;
+-1.7895e-06;
+-7.95375e-06;
+-1.76175e-06;
+-1.0596e-05;
+-5.24553e-06;
+-1.02345e-05;
+-1.74422e-06;
+-2.16064e-05;
+-1.85059e-06;
+-2.24752e-05;
+-3.34957e-07;
+-1.34774e-05;
+-2.51353e-06;
+-2.18312e-05;
+-1.04202e-06;
+-1.99898e-05;
+-1.463e-06;
+-1.94256e-05;
+-2.11821e-06;
+-5e-05;
+5.03269e-07;
+-2.00631e-05;
+-2.41829e-06;
+-2.11637e-05;
+-1.12298e-06;
+-1.61508e-05;
+-1.11752e-05;
+-1.7267e-05;
+-6.9337e-06;
+-1.43931e-05;
+-6.46903e-06;
+-1.70354e-05;
+-2.45177e-06;
+-5e-05;
+-2.50465e-06;
+-5e-05;
+-1.21588e-05;
+-2.01431e-05;
+-1.19797e-05;
+];
+
+d = 0.00075745;
+A2 = [...
+];
+
+f2 = [...
+0;
+0;
+0;
+0;
+-3.01994e-06;
+-1.76298e-07;
+0;
+0;
+-4.78204e-06;
+-2.0995e-07;
+-5.01951e-06;
+1.22653e-06;
+0;
+0;
+-4.98476e-06;
+-3.57931e-06;
+0;
+0;
+-7.06572e-06;
+-4.97786e-06;
+-5.40025e-06;
+-1.91591e-06;
+0;
+0;
+-3.26344e-06;
+-9.4196e-07;
+0;
+0;
+-3.47107e-06;
+1.69693e-06;
+0;
+0;
+0;
+0;
+0;
+0;
+-8.36624e-06;
+-3.21786e-07;
+0;
+0;
+-9.19771e-06;
+-4.8646e-07;
+-9.46319e-06;
+2.13218e-06;
+0;
+0;
+-1.03261e-05;
+-6.78209e-06;
+0;
+0;
+-1.09634e-05;
+-7.09926e-06;
+-1.01615e-05;
+-3.77678e-06;
+-9.09151e-06;
+-6.13484e-06;
+-9.06507e-06;
+4.28268e-06;
+-6.83374e-06;
+-2.4706e-07;
+-8.51881e-06;
+-3.22169e-06;
+-7.83369e-06;
+3.49294e-06;
+-8.45665e-06;
+3.90404e-06;
+-7.913e-06;
+1.91502e-06;
+-7.13457e-06;
+1.77604e-06;
+-7.71821e-06;
+-2.87949e-06;
+-7.58455e-06;
+-3.19117e-07;
+0;
+0;
+0;
+0;
+-1.17905e-05;
+-6.57623e-07;
+0;
+0;
+-1.26738e-05;
+-7.75536e-07;
+-1.28663e-05;
+2.59111e-06;
+0;
+0;
+-1.33726e-05;
+-8.23012e-06;
+0;
+0;
+-1.40382e-05;
+-8.4827e-06;
+-1.3445e-05;
+-4.60248e-06;
+-1.21277e-05;
+-7.7075e-06;
+-1.20933e-05;
+5.59837e-06;
+-1.00511e-05;
+-5.83173e-07;
+-1.16989e-05;
+-4.19969e-06;
+-1.07938e-05;
+5.03904e-06;
+-1.14551e-05;
+5.33497e-06;
+-1.1057e-05;
+2.39173e-06;
+-1.01977e-05;
+2.22842e-06;
+-1.08385e-05;
+-3.93909e-06;
+-1.09119e-05;
+-6.59537e-07;
+-1.51572e-05;
+6.80143e-06;
+-1.97651e-05;
+7.92292e-06;
+-1.3582e-05;
+-8.81729e-07;
+-1.73372e-05;
+7.52177e-06;
+-1.63794e-05;
+2.79237e-06;
+-1.33659e-05;
+2.38778e-06;
+-5e-05;
+7.93536e-06;
+-1.92072e-05;
+-9.96979e-07;
+-1.93357e-05;
+3.15138e-06;
+-1.54595e-05;
+-9.11795e-06;
+-1.68142e-05;
+-5.1397e-06;
+-1.4113e-05;
+-4.66534e-06;
+-1.63265e-05;
+-9.49749e-07;
+-5e-05;
+-1.0631e-06;
+-5e-05;
+-1.02021e-05;
+-1.98778e-05;
+-9.99727e-06;
+];
+
+d2 = 0;
+
+sum(sum(abs(f-f2)>1e-14))
+
+
+sum(sum(abs(A-A2)>1e-14))
+
+
+alpha0 = 100;
+V = 1e-5; %1e-20;
+V0 = 5e-5;
+L = 2.25e-5;
+
+alpha = @(t) log(exp(alpha0 - V.*t./L ) - V0 * (exp(-t./L.*V)-1)./V);
+
+x = [0:0.001:0.1];
+alpha(x)
+plot(x, alpha(x));
+
+%abs(d-d2)
+
+% x = [...
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% -1.2337e-10;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% -1.2337e-10;
+% 0;
+% -1.2337e-10;
+% 0;
+% 0;
+% 0;
+% ];
+% 
+% ignore = [...
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0; 
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 1;
+% 1; 
+% 0;
+% 0; 
+% 1; 
+% 1; 
+% 0;
+% 0; 
+% 1;
+% 1; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 1;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0; 
+% 0; 
+% 0;
+% 0; 
+% 0; 
+% 0; 
+% 0; 
+% 0; 
+% 1;
+% 0; 
+% 1;
+% 0; 
+% 0;
+% 0; 
+% ];
+% 
+% newA = A;
+% newf = f;
+% for i=1:length(ignore)
+%     if (ignore(i)==1)
+%         for j=1:length(ignore)
+%             newA(i,j) = (i==j);
+%         end
+%         newf(i) = x(i);
+%     end
+% end
+% newA
+% newf
+
+% ignore = zeros(1, length(f));
+% for i=1:length(ignore)
+%     if (A(i,i)==0)
+%         ignore(1, i) = 1;
+%     end
+% end
+% 
+% indices = [];
+% for i=1:length(ignore)
+%     if (ignore(i)==0)
+%         indices = [indices i];
+%     end
+% end
+% 
+% len = length(indices);
+% newA = zeros(len, len);
+% newf = zeros(1, len);
+% for i=1:len
+%     newf(i) = f(indices(i));
+% 
+%     for j=1:len
+%         newA(i,j) = A(indices(i), indices(j));
+%     end
+% end
+% newf;
+% newA
+
+%inv(A);
+
+% sol = newA\newf
+% 
+% tnnmgSol = [...
+%     -1.03349;
+% -0.0502635;
+% 0.0161336;
+% -0.00852914;
+% -0.0130659;
+% 0.00660115;
+% -0.434915;
+% -0.0128241;
+% -0.00893317;
+% 0.00884256;
+% -0.0110232;
+% -0.10936;
+% 1.14978;
+% 0.0581132;
+% 0.019598;
+% 0.190758;
+% 0.478393;
+% -0.00433049;
+% 0.0513446;
+% 0.353218;
+% -0.00482644;
+% 0.11758;
+% 0;
+% 0;
+% -0.013677;
+% 0.0512563;
+% 0;
+% 0;
+% -0.00514197;
+% -0.117664;
+% 0;
+% 0;
+% 0.130034;
+% 0.574041;
+% 0.0136965;
+% 0.214953;
+% 0.0335986;
+% -0.00703477;
+% 0.0438202;
+% 0.346506;
+% -0.000896805;
+% 0.127837;
+% 0.041946;
+% 0.217409;
+% -1.2337e-10;
+% 0.112528;
+% -0.00828618;
+% 0.00453345;
+% -0.00898709;
+% 0.0781683;
+% 0.092355;
+% -0.556143;
+% -0.00443686;
+% -0.140953;
+% 0.0424371;
+% -0.250082;
+% -0.00336878;
+% 0.00105471;
+% -1.2337e-10;
+% 0.0188452;
+% -1.2337e-10;
+% -0.14111;
+% 0.00298499;
+% -0.190904;
+% ];
+% 
+% sol-tnnmgSol
\ No newline at end of file

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

+add_subdirectory("data-structures")
+add_subdirectory("factories")
+add_subdirectory("io")
+add_subdirectory("multi-threading")
+add_subdirectory("problem-data")
+add_subdirectory("spatial-solving")
+add_subdirectory("tests")
+add_subdirectory("time-stepping")
+add_subdirectory("utils")
+
+add_custom_target(tectonic_dune SOURCES
+  assemblers.hh
+  assemblers.cc 
+  explicitgrid.hh
+  explicitvectors.hh
+  gridselector.hh
+)
+
+#install headers
 install(FILES
-  body.hh
-  frictiondata.hh
-  frictionpotential.hh
-  globalfrictiondata.hh
-  globalfriction.hh
-  globalratestatefriction.hh
-  gravity.hh
-  localfriction.hh
-  minimisation.hh
-  myblockproblem.hh
-  mydirectionalconvexfunction.hh
-  quadraticenergy.hh
-  tectonic.hh
-  DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
+  assemblers.hh
+  explicitgrid.hh
+  explicitvectors.hh
+  gridselector.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

src/assemblers.cc → dune/tectonic/assemblers.cc

@@ -4,7 +4,8 @@
 
 #include <dune/istl/scaledidmatrix.hh>
 
-#include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
+#include <dune/fufem/assemblers/boundaryoperatorassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/normalstressboundaryassembler.hh>
@@ -12,41 +13,49 @@
 #include <dune/fufem/assemblers/localassemblers/variablecoefficientviscosityassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/weightedmassassembler.hh>
-#include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
 #include <dune/fufem/functiontools/p0p1interpolation.hh>
 #include <dune/fufem/quadraturerules/quadraturerulecache.hh>
 
-#include <dune/tectonic/frictionpotential.hh>
-#include <dune/tectonic/globalratestatefriction.hh>
+#include "data-structures/friction/frictionpotential.hh"
+#include "data-structures/friction/globalratestatefriction.hh"
 
 #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),
       vertexAssembler(vertexBasis, vertexBasis) {}
 
+template <class GridView, int dimension>
+template <class LocalBoundaryFunctionalAssemblerType, class GlobalVectorType>
+void MyAssembler<GridView, dimension>::assembleBoundaryFunctional(LocalBoundaryFunctionalAssemblerType& localAssembler,
+                                GlobalVectorType& b,
+                                const BoundaryPatch<GridView>& boundaryPatch,
+                                bool initializeVector) const {
+
+    vertexAssembler.assembleBoundaryFunctional(localAssembler, b, boundaryPatch, initializeVector);
+}
+
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
-    BoundaryPatch<GridView> const &frictionalBoundary,
-    ScalarMatrix &frictionalBoundaryMass) const {
-  BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
-                        LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
-      frictionalBoundaryMassAssembler(frictionalBoundary);
-  vertexAssembler.assembleOperator(frictionalBoundaryMassAssembler,
-                                   frictionalBoundaryMass);
+        BoundaryPatch<GridView> const &frictionalBoundary,
+        ScalarMatrix &frictionalBoundaryMass) const {
+
+  MassAssembler<Grid, LocalVertexBasis, LocalVertexBasis> localMassAssembler;
+  const BoundaryOperatorAssembler<VertexBasis, VertexBasis> boundaryAssembler(vertexBasis, vertexBasis, frictionalBoundary);
+  boundaryAssembler.assemble(localMassAssembler, frictionalBoundaryMass);
 }
 
 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 +67,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 +79,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 +96,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,11 +106,33 @@ 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 {
+        const BoundaryPatch<GridView>& neumannBoundary,
+        const Dune::BitSetVector<dimension>& neumannNodes,
+        Vector& f,
+        const Dune::VirtualFunction<double, double>& neumann,
+        double relativeTime) const {
+
+  typename LocalVector::field_type val = 0;
+  neumann.evaluate(relativeTime, val);
+
+  size_t idx = 0;
+  bool neumannIdx = false;
+  for (; idx<neumannNodes.size() && !neumannIdx; idx++) {
+    for (size_t d=0; d<dimension; d++) {
+        if (neumannNodes[idx][d]) {
+            neumannIdx = true;
+            break;
+        }
+    }
+  }
+  idx--;
+
   LocalVector localNeumann(0);
-  neumann.evaluate(relativeTime, localNeumann[0]);
+  for (size_t i=0; i<localNeumann.size(); i++) {
+      if (neumannNodes[idx][i])
+          localNeumann[i] = val;
+  }
+
   NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
       std::make_shared<ConstantFunction<LocalVector, LocalVector>>(
           localNeumann));
@@ -106,9 +142,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);
@@ -136,41 +174,13 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
         std::fmin(normals[i] * nodalTractionAverage[i], 0) * weights[i];
 }
 
-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>> {
-  // Lumping of the nonlinearity
-  ScalarVector weights;
-  {
-    NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
-        frictionalBoundaryAssembler(std::make_shared<
-            ConstantFunction<LocalVector, typename ScalarVector::block_type>>(
-            1));
-    vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
-                                               weights, frictionalBoundary);
-  }
-  switch (frictionModel) {
-    case Config::Truncated:
-      return std::make_shared<GlobalRateStateFriction<
-          Matrix, Vector, TruncatedRateState, GridView>>(
-          frictionalBoundary, frictionInfo, weights, weightedNormalStress);
-    case Config::Regularised:
-      return std::make_shared<GlobalRateStateFriction<
-          Matrix, Vector, RegularisedRateState, GridView>>(
-          frictionalBoundary, frictionInfo, weights, weightedNormalStress);
-    default:
-      assert(false);
-  }
-}
-
 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);

dune/tectonic/assemblers.hh

+#ifndef SRC_ASSEMBLERS_HH
+#define SRC_ASSEMBLERS_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.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/fufem/boundarypatch.hh>
+
+#include "data-structures/friction/globalfriction.hh"
+#include "data-structures/friction/globalfrictiondata.hh"
+
+#include "data-structures/enums.hh"
+
+template <class GridView, int dimension>
+class MyAssembler {
+public:
+    using GV = GridView;
+    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 CellBasis = P0Basis<GridView, double>;
+    using VertexBasis = P1NodalBasis<GridView, double>;
+
+    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 LocalCellBasis = typename CellBasis::LocalFiniteElement;
+    using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
+
+    Assembler<CellBasis, CellBasis> cellAssembler;
+    Assembler<VertexBasis, VertexBasis> vertexAssembler;
+
+public:
+    MyAssembler(GridView const &gridView);
+
+    template <class LocalBoundaryFunctionalAssemblerType, class GlobalVectorType>
+    void assembleBoundaryFunctional(LocalBoundaryFunctionalAssemblerType& localAssembler,
+                                    GlobalVectorType& b,
+                                    const BoundaryPatch<GridView>& boundaryPatch,
+                                    bool initializeVector=true) const;
+
+    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(
+            const BoundaryPatch<GridView>& neumannBoundary,
+            const Dune::BitSetVector<dimension>& neumannNodes,
+            Vector& f,
+            const Dune::VirtualFunction<double, double>& 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

dune/tectonic/assemblers_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "explicitgrid.hh"
+#include "explicitvectors.hh"
+
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/boundarypatch.hh>
+
+#include "assemblers.hh"
+
+template class MyAssembler<DefLeafGridView, MY_DIM>;
+
+
+using MyNeumannBoundaryAssembler = NeumannBoundaryAssembler<DeformedGrid, typename ScalarVector::block_type>;
+
+template void MyAssembler<DefLeafGridView, MY_DIM>::assembleBoundaryFunctional<MyNeumannBoundaryAssembler, ScalarVector>(
+        MyNeumannBoundaryAssembler& localAssembler,
+        ScalarVector& b,
+        const BoundaryPatch<DefLeafGridView>& boundaryPatch,
+        bool initializeVector) const;

dune/tectonic/data-structures/.cc

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

dune/tectonic/data-structures/CMakeLists.txt

+add_subdirectory("body")
+add_subdirectory("friction")
+add_subdirectory("network")
+
+add_custom_target(tectonic_dune_data-structures SOURCES
+  enumparser.hh
+  enumparser.cc 
+  enums.hh
+  matrices.hh
+  program_state.hh
+)
+
+#install headers
+install(FILES
+  enumparser.hh
+  enums.hh
+  matrices.hh
+  program_state.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/data-structures/body/CMakeLists.txt

+add_custom_target(tectonic_dune_data-structures_body SOURCES
+  body.hh
+  body.cc 
+  bodydata.hh
+  boundarycondition.hh
+)
+
+#install headers
+install(FILES
+  body.hh
+  bodydata.hh
+  boundarycondition.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/data-structures/body/body.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/hybridutilities.hh>
+#include "body.hh"
+
+// -------------------------------
+// --- LeafBody Implementation ---
+// -------------------------------
+template <class GridTEMPLATE, class VectorType>
+LeafBody<GridTEMPLATE, VectorType>::LeafBody(
+        const std::shared_ptr<BodyData<dim>>& bodyData,
+        const std::shared_ptr<GridTEMPLATE>& hostGrid) :
+    bodyData_(bodyData),
+    hostGrid_(hostGrid),
+    deformation_(std::make_shared<DeformationFunction>(hostGrid_->leafGridView())),
+    grid_(std::make_shared<GridType>(*hostGrid_, *deformation_)),
+    gridView_(grid_->leafGridView()),
+    assembler_(std::make_shared<Assembler>(gridView_)),
+    matrices_()
+{
+
+    boundaryConditions_.resize(0);
+
+    externalForce_ = [&](double relativeTime, VectorType &ell) {
+        // Problem formulation: right-hand side
+        std::vector<std::shared_ptr<BoundaryCondition>> leafNeumannConditions;
+        boundaryConditions("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(),
+                                        *leafNeumannCondition->boundaryNodes(),
+                                        b,
+                                        *leafNeumannCondition->boundaryFunction(),
+                                        relativeTime);
+            ell += b;
+        }
+    };
+}
+
+
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::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 forces
+    assembler_->assembleBodyForce(bodyData_->getGravityField(), gravityFunctional_);
+}
+
+// setter and getter
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::data() const
+-> const std::shared_ptr<BodyData<dim>>& {
+
+    return bodyData_;
+}
+
+
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::setDeformation(const VectorType& def) {
+    deformation_->setDeformation(def);
+}
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::deformation() const
+-> DeformationFunction& {
+
+    return *deformation_;
+}
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::grid() const
+-> std::shared_ptr<GridType> {
+
+    return grid_;
+}
+
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::gridView() const
+-> const GridView& {
+
+    return gridView_;
+}
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::nVertices() const
+-> size_t {
+    return gridView_.size(dim);
+}
+
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::assembler() const
+-> const std::shared_ptr<Assembler>& {
+
+    return assembler_;
+}
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::matrices() const
+-> const Matrices<Matrix,1>& {
+
+    return matrices_;
+}
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::externalForce() const
+-> const ExternalForce& {
+
+    return externalForce_;
+}
+
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition) {
+    boundaryConditions_.emplace_back(boundaryCondition);
+}
+
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryConditions(
+        const std::string& tag,
+        std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const {
+
+    selectedConditions.resize(0);
+    for (size_t i=0; i<boundaryConditions_.size(); i++) {
+        if (boundaryConditions_[i]->tag() == tag)
+            selectedConditions.emplace_back(boundaryConditions_[i]);
+    }
+}
+
+template <class GridTEMPLATE, class VectorType>
+auto LeafBody<GridTEMPLATE, VectorType>::boundaryConditions() const
+-> const std::vector<std::shared_ptr<BoundaryCondition>>& {
+
+    return boundaryConditions_;
+}
+
+template <class GridTEMPLATE, class VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryPatchNodes(
+        const std::string& tag,
+        BoundaryPatchNodes& nodes) const {
+
+    std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(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 VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryNodes(
+        const std::string& tag,
+        BoundaryNodes& nodes) const {
+
+    std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(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 VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryFunctions(
+        const std::string& tag,
+        BoundaryFunctions& functions) const {
+
+    std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(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 VectorType>
+void LeafBody<GridTEMPLATE, VectorType>::boundaryPatches(
+        const std::string& tag,
+        BoundaryPatches& patches) const {
+
+    std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(tag, _boundaryConditions);
+
+    patches.resize(_boundaryConditions.size());
+
+    for (size_t i=0; i<patches.size(); i++) {
+        patches[i] = _boundaryConditions[i]->boundaryPatch().get();
+    }
+}
+
+
+// ---------------------------
+// --- Body Implementation ---
+// ---------------------------
+
+// setter and getter
+
+template <class GridView>
+auto Body<GridView>::data() const
+-> const std::shared_ptr<BodyData<dim>>& {
+
+    return bodyData_;
+}
+
+
+template <class GridView>
+auto Body<GridView>::grid() const
+-> std::shared_ptr<GridType> {
+
+    return grid_;
+}
+
+template <class GridView>
+auto Body<GridView>::gridView() const
+-> const GridView& {
+
+    return gridView_;
+}
+
+template <class GridView>
+auto Body<GridView>::nVertices() const
+-> size_t {
+    return gridView_.size(dim);
+}
+
+
+template <class GridView>
+void Body<GridView>::addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition) {
+
+    boundaryConditions_.push_back(boundaryCondition);
+}
+
+
+template <class GridView>
+void Body<GridView>::boundaryConditions(
+        const std::string& tag,
+        std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const {
+
+    selectedConditions.resize(0);
+    for (size_t i=0; i<boundaryConditions_.size(); i++) {
+        if (boundaryConditions_[i]->tag() == tag)
+            selectedConditions.push_back(boundaryConditions_[i]);
+    }
+}
+
+
+template <class GridView>
+auto Body<GridView>::boundaryConditions() const
+-> const std::vector<std::shared_ptr<BoundaryCondition>>& {
+
+    return boundaryConditions_;
+}
+
+
+template <class GridView>
+void Body<GridView>::boundaryPatchNodes(
+        const std::string& tag,
+        BoundaryPatchNodes& nodes) const {
+
+    std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(tag, _boundaryConditions);
+
+    nodes.resize(_boundaryConditions.size());
+
+    for (size_t i=0; i<nodes.size(); i++) {
+        nodes[i] = _boundaryConditions[i]->boundaryPatch()->getVertices();
+    }
+}
+
+template <class GridView>
+void Body<GridView>::boundaryNodes(
+        const std::string& tag,
+        BoundaryNodes& nodes) const {
+
+    std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(tag, _boundaryConditions);
+
+    nodes.resize(_boundaryConditions.size());
+
+    for (size_t i=0; i<nodes.size(); i++) {
+        nodes[i] = _boundaryConditions[i]->boundaryNodes().get();
+    }
+}
+
+template <class GridView>
+void Body<GridView>::boundaryFunctions(
+        const std::string& tag,
+        BoundaryFunctions& functions) const {
+
+    std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(tag, _boundaryConditions);
+
+    functions.resize(_boundaryConditions.size());
+
+    for (size_t i=0; i<functions.size(); i++) {
+        functions[i] = _boundaryConditions[i]->boundaryFunction().get();
+    }
+}
+
+template <class GridView>
+void Body<GridView>::boundaryPatches(
+        const std::string& tag,
+        BoundaryPatches& patches) const {
+
+    std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
+    this->boundaryConditions(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"

dune/tectonic/data-structures/body/body.hh

+#ifndef SRC_DATA_STRUCTURES_BODY_HH
+#define SRC_DATA_STRUCTURES_BODY_HH
+
+#include <functional>
+#include <type_traits>
+
+#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 "../../assemblers.hh"
+#include "../enums.hh"
+#include "../matrices.hh"
+
+#include "boundarycondition.hh"
+#include "bodydata.hh"
+
+template <class HostGridType, class VectorType>
+class LeafBody {
+public:
+    static const int dim = HostGridType::dimensionworld;
+
+    using DeformationFunction = DeformationFunction<typename HostGridType::LeafGridView, VectorType>;
+    using GridType = Dune::GeometryGrid<HostGridType, DeformationFunction>;
+    using GridView = typename GridType::LeafGridView;
+
+    using Function = Dune::VirtualFunction<double, double>;
+    using ExternalForce = std::function<void(double, VectorType &)>;
+
+    using Assembler = MyAssembler<GridView, dim>;
+    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 BoundaryCondition = BoundaryCondition<GridView, dim>;
+
+    using BoundaryFunctions = std::vector<const Function* >;
+    using BoundaryNodes = std::vector<const Dune::BitSetVector<dim>* >;
+    using BoundaryPatchNodes = std::vector<const Dune::BitSetVector<1>* >;
+    using BoundaryPatches = std::vector<const typename BoundaryCondition::BoundaryPatch* >;
+
+    using StateEnergyNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
+
+private:
+    std::shared_ptr<BodyData<dim>> bodyData_;
+
+    std::shared_ptr<HostGridType> hostGrid_;
+    std::shared_ptr<DeformationFunction> deformation_;
+    std::shared_ptr<GridType> grid_;
+    GridView gridView_;
+
+    std::shared_ptr<Assembler> assembler_;
+    Matrices<Matrix,1> matrices_;
+
+    ExternalForce externalForce_;
+
+    VectorType gravityFunctional_;
+
+    // boundary conditions
+    std::vector<std::shared_ptr<BoundaryCondition>> boundaryConditions_;
+
+public:
+    LeafBody(
+            const std::shared_ptr<BodyData<dim>>& bodyData,
+            const std::shared_ptr<HostGridType>& hostGrid);
+
+    void assemble();
+
+    // setter and getter
+    auto data() const -> const std::shared_ptr<BodyData<dim>>&;
+
+    void setDeformation(const VectorType& def);
+    auto deformation() const -> DeformationFunction&;
+    auto grid() const -> std::shared_ptr<GridType>;
+
+    auto gridView() const -> const GridView&;
+    auto nVertices() const -> size_t;
+
+    auto assembler() const -> const std::shared_ptr<Assembler>&;
+    auto matrices() const -> const Matrices<Matrix,1>&;
+
+    auto externalForce() const -> const ExternalForce&;
+
+    void addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition);
+    void boundaryConditions(
+            const std::string& tag,
+            std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const;
+    auto boundaryConditions() const -> const std::vector<std::shared_ptr<BoundaryCondition>>&;
+
+    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;
+};
+
+template <class GridView>
+class Body {
+public:
+    enum {dim = GridView::dimensionworld};
+
+    using GridType = typename GridView::Grid;
+
+    using Function = Dune::VirtualFunction<double, double>;
+
+    using BoundaryCondition = BoundaryCondition<GridView, dim>;
+
+    using BoundaryFunctions = std::vector<const Function* >;
+    using BoundaryNodes = std::vector<const Dune::BitSetVector<dim>* >;
+    using BoundaryPatchNodes = std::vector<const Dune::BitSetVector<1>* >;
+    using BoundaryPatches = std::vector<const typename BoundaryCondition::BoundaryPatch* >;
+
+private:
+    std::shared_ptr<BodyData<dim>> bodyData_;
+
+    std::shared_ptr<GridType> grid_;
+    const size_t level_;
+    GridView gridView_;
+
+    // boundary conditions
+    std::vector<std::shared_ptr<BoundaryCondition>> boundaryConditions_;
+
+public:
+    template <class GridView_ = GridView, std::enable_if_t<std::is_same<GridView_, typename GridType::LeafGridView>::value, int> = 0>
+    Body(
+            const std::shared_ptr<BodyData<dim>>& bodyData,
+            const std::shared_ptr<GridType>& grid) :
+        bodyData_(bodyData),
+        grid_(grid),
+        level_(grid_->maxLevel()),
+        gridView_(grid_->leafGridView()) {
+
+        boundaryConditions_.resize(0);
+    }
+
+    template <class GridView_ = GridView, std::enable_if_t<std::is_same<GridView_, typename GridType::LevelGridView>::value, int> = 0>
+    Body(
+            const std::shared_ptr<BodyData<dim>>& bodyData,
+            const std::shared_ptr<GridType>& grid,
+            const size_t level) :
+        bodyData_(bodyData),
+        grid_(grid),
+        level_(level),
+        gridView_(grid_->levelGridView(level_)) {
+
+        boundaryConditions_.resize(0);
+    }
+
+
+    // setter and getter
+    auto data() const -> const std::shared_ptr<BodyData<dim>>&;
+    auto grid() const -> std::shared_ptr<GridType>;
+
+    //template <class GridView_ = GridView, std::enable_if_t<std::is_same<GridView_, typename GridType::LevelGridView>::value, int> = 0>
+    auto level() const
+    -> size_t {
+
+        return level_;
+    }
+
+    auto gridView() const -> const GridView&;
+    auto nVertices() const -> size_t;
+
+    void addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition);
+    void boundaryConditions(
+            const std::string& tag,
+            std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const;
+    auto boundaryConditions() const -> const std::vector<std::shared_ptr<BoundaryCondition>>&;
+
+    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

dune/tectonic/data-structures/body/body_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+#include "../../explicitvectors.hh"
+
+template class LeafBody<Grid, Vector>;
+
+template class Body<DefLeafGridView>;
+template class Body<DefLevelGridView>;

dune/tectonic/body.hh → dune/tectonic/data-structures/body/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/data-structures/body/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>
+
+#include <dune/tectonic/utils/tobool.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),
+      boundaryFunction_(function)
+  {
+      setBoundaryPatch(patch);
+  }
+
+  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;
+
+      auto nodes = patch->getVertices();
+      boundaryNodes_ = std::make_shared<Dune::BitSetVector<dims>>(nodes->size(), false);
+      for (size_t i=0; i<nodes->size(); i++) {
+          if (toBool((*nodes)[i])) {
+             for (size_t d=0; d<dims; d++) {
+                 (*boundaryNodes_)[i][d] = true;
+             }
+          }
+      }
+  }
+
+  void setBoundaryPatch(const BoundaryPatch& patch) {
+      auto patch_ptr = std::make_shared<BoundaryPatch>(patch);
+      setBoundaryPatch(patch_ptr);
+  }
+
+  void setBoundaryFunction(std::shared_ptr<Function> function) {
+      boundaryFunction_ = function;
+  }
+
+  void set(std::shared_ptr<BoundaryPatch> patch, std::shared_ptr<Function> function) {
+      setBoundaryPatch(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_;
+  }
+};
+
+#endif

src/enumparser.hh → dune/tectonic/data-structures/enumparser.hh

-#ifndef SRC_ENUMPARSER_HH
-#define SRC_ENUMPARSER_HH
+#ifndef SRC_DATA_STRUCTURES_ENUMPARSER_HH
+#define SRC_DATA_STRUCTURES_ENUMPARSER_HH
 
 // Copyright Carsten Graeser 2012
 

dune/tectonic/data-structures/friction/CMakeLists.txt

+add_custom_target(tectonic_dune_data-structures_friction SOURCES
+  frictioncouplingpair.hh
+  frictiondata.hh 
+  frictionpotential.hh
+  globalfriction.hh
+  globalfrictiondata.hh
+  globalratestatefriction.hh
+  localfriction.hh
+)
+
+#install headers
+install(FILES
+  frictioncouplingpair.hh
+  frictiondata.hh 
+  frictionpotential.hh
+  globalfriction.hh
+  globalfrictiondata.hh
+  globalratestatefriction.hh
+  localfriction.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)