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"
+
+#include <dune/tectonic/utils/debugutils.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);
+    //print(*matrices_.elasticity, "elasticity");
+
+    assembler_->assembleViscosity(bodyData_->getShearViscosityField(), bodyData_->getBulkViscosityField(), *matrices_.damping);
+    //print(*matrices_.damping, "viscosity");
+
+    assembler_->assembleMass(bodyData_->getDensityField(), *matrices_.mass);
+    //print(*matrices_.mass, "mass");
+
+    // assemble forces
+    assembler_->assembleBodyForce(bodyData_->getGravityField(), gravityFunctional_);
+    //print(gravityFunctional_, "gravity");
+}
+
+// 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>::hostGrid() const
+-> std::shared_ptr<GridTEMPLATE> {
+
+    return hostGrid_;
+}
+
+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 hostGrid() const -> std::shared_ptr<HostGridType>;
+
+    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.cc → dune/tectonic/data-structures/enumparser.cc

@@ -32,6 +32,12 @@ Config::FrictionModel StringToEnum<Config::FrictionModel>::convert(
   if (s == "Regularised")
     return Config::Regularised;
 
+  if (s == "Tresca")
+    return Config::Tresca;
+
+  if (s == "None")
+    return Config::None;
+
   DUNE_THROW(Dune::Exception, "failed to parse enum");
 }
 

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
 

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

@@ -2,7 +2,7 @@
 #define SRC_ENUMS_HH
 
 struct Config {
-  enum FrictionModel { Truncated, Regularised };
+  enum FrictionModel { Truncated, Regularised, Tresca, None };
   enum stateModel { AgeingLaw, SlipLaw };
   enum scheme { Newmark, BackwardEuler };
   enum PatchType { Rectangular, Trapezoidal };

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)

dune/tectonic/data-structures/friction/frictioncouplingpair.hh

+#ifndef SRC_FRICTIONCOUPLINGPAIR_HH
+#define SRC_FRICTIONCOUPLINGPAIR_HH
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/contact/common/couplingpair.hh>
+
+#include "globalfrictiondata.hh"
+
+template <class GridType, class LocalVectorType, class field_type = double>
+class FrictionCouplingPair : public Dune::Contact::CouplingPair<GridType, GridType, field_type>{
+private:
+    static const int dim = GridType::dimensionworld;
+
+    using Base = Dune::Contact::CouplingPair<GridType,GridType,field_type>;
+    using LocalVector = LocalVectorType;
+
+    // friction data
+    std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch_;
+    std::shared_ptr<GlobalFrictionData<dim>> frictionData_;
+
+public:
+  void setWeakeningPatch(std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch) {
+      weakeningPatch_ = weakeningPatch;
+  }
+
+  void setFrictionData(std::shared_ptr<GlobalFrictionData<dim>> frictionData) {
+      frictionData_ = frictionData;
+  }
+
+  const auto& weakeningPatch() const {
+      return *weakeningPatch_;
+  }
+
+  const auto& frictionData() const {
+      return *frictionData_;
+  }
+};
+#endif

dune/tectonic/frictionpotential.hh → dune/tectonic/data-structures/friction/frictionpotential.hh

@@ -9,7 +9,7 @@
 #include <dune/common/exceptions.hh>
 #include <dune/common/function.hh>
 
-#include <dune/tectonic/frictiondata.hh>
+#include "frictiondata.hh"
 
 class FrictionPotential {
 public:
@@ -37,30 +37,53 @@ class TruncatedRateState : public FrictionPotential {
     if (V <= Vmin)
       return 0.0;
 
-    return fd.a * std::log(V / Vmin);
+    //std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;
+
+    auto res = fd.a * std::log(V / Vmin);
+    if (std::isinf(res)) {
+        //std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;
+    }
+
+    return res;
   }
 
   double differential(double V) const override {
+    //std::cout << "differential: " <<  weight * fd.C - weightedNormalStress * coefficientOfFriction(V) << std::endl;
+    //if (V <= Vmin or regularity(V)>10e8)
+    //  return 0.0;
+
     return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
   }
 
   double second_deriv(double V) const override {
+    //std::cout << "second deriv: " << ((V<=Vmin) ? 0.0 : -weightedNormalStress * (fd.a / V)) << std::endl;
+
     if (V <= Vmin)
-      return 0;
+      return 0.0;
 
-    return -weightedNormalStress * (fd.a / V);
+    return - weightedNormalStress * (fd.a / V);
   }
 
   double regularity(double V) const override {
+    //std::cout << "regularity: " << ((std::abs(V - Vmin) < 1e-14) ? std::numeric_limits<double>::infinity() : std::abs(second_deriv(V))) << std::endl;
+
     if (std::abs(V - Vmin) < 1e-14) // TODO
       return std::numeric_limits<double>::infinity();
 
     return std::abs(second_deriv(V));
   }
 
+  double evaluate(double V) const override {
+    if (V <= Vmin)
+        return 0.0;
+
+      return weight * fd.C * V - weightedNormalStress * fd.a * V * (std::log(V / Vmin) - 1);
+  }
+
   void updateAlpha(double alpha) override {
     double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
     Vmin = fd.V0 / std::exp(logrest);
+    //Vmin = 0.0;
   }
 
 private:
@@ -92,6 +115,10 @@ class RegularisedRateState : public FrictionPotential {
     return std::abs(second_deriv(V));
   }
 
+  double evaluate(double V) const override {
+      return weight * fd.C * V - weightedNormalStress * 4 * fd.a * Vmin * std::pow(std::asinh(0.25 * V / Vmin), 2.0);
+  }
+
   void updateAlpha(double alpha) override {
     double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
     Vmin = fd.V0 / std::exp(logrest);
@@ -104,8 +131,48 @@ class RegularisedRateState : public FrictionPotential {
   double Vmin;
 };
 
+class Tresca : public FrictionPotential {
+public:
+  Tresca(double _weight, double _weightedNormalStress, FrictionData _fd)
+      : fd(_fd), weightedNormalStress(_weightedNormalStress) {}
+
+  double coefficientOfFriction(double V) const override {
+    return fd.mu0;
+  }
+
+  double differential(double V) const override {
+    //if (std::isinf(regularity(V)))
+    //    return 0.0;
+    return - weightedNormalStress * coefficientOfFriction(V);
+  }
+
+  double second_deriv(double V) const override {
+    return 0;
+  }
+
+  double regularity(double V) const override {
+    if (std::abs(V) < 1e-14) // TODO
+        return std::numeric_limits<double>::infinity();
+
+    return std::abs(second_deriv(V));
+  }
+
+  double evaluate(double V) const override {
+      return - weightedNormalStress * fd.mu0 * V;
+  }
+
+  void updateAlpha(double alpha) override {}
+
+private:
+  FrictionData const fd;
+  double const weightedNormalStress;
+};
+
 class ZeroFunction : public FrictionPotential {
 public:
+  template <typename... Args>
+  ZeroFunction(Args... args) {}
+
   double evaluate(double) const override { return 0; }
 
   double coefficientOfFriction(double s) const override { return 0; }

dune/tectonic/globalfriction.hh → dune/tectonic/data-structures/friction/globalfriction.hh

@@ -9,9 +9,10 @@
 
 #include <dune/solvers/common/interval.hh>
 
-#include <dune/tectonic/localfriction.hh>
+#include "localfriction.hh"
 
-template <class Matrix, class Vector> class GlobalFriction {
+template <class Matrix, class Vector>
+class GlobalFriction {
 protected:
   using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
 
@@ -38,8 +39,9 @@ template <class Matrix, class Vector> class GlobalFriction {
   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)
+    for (size_t i = 0; i < v.size(); ++i) {
       restriction(i).addHessian(v[i], hessian[i][i]);
+    }
   }
 
   void directionalDomain(Vector const &, Vector const &,
@@ -76,11 +78,13 @@ template <class Matrix, class Vector> class GlobalFriction {
 
   ScalarVector coefficientOfFriction(Vector const &x) const {
     ScalarVector ret(x.size());
-    for (size_t i = 0; i < x.size(); ++i)
+    for (size_t i = 0; i < x.size(); ++i) {
       ret[i] = restriction(i).coefficientOfFriction(x[i]);
+    }
     return ret;
   }
 
-  void virtual updateAlpha(ScalarVector const &alpha) = 0;
+  void virtual updateAlpha(const std::vector<ScalarVector>& alpha) = 0;
+
 };
 #endif

dune/tectonic/globalfrictiondata.hh → dune/tectonic/data-structures/friction/globalfrictiondata.hh

@@ -4,7 +4,7 @@
 #include <dune/common/function.hh>
 #include <dune/common/fvector.hh>
 
-#include <dune/tectonic/frictiondata.hh>
+#include "frictiondata.hh"
 
 template <class DomainType>
 double evaluateScalarFunction(
@@ -27,6 +27,7 @@ template <int dimension> class GlobalFrictionData {
       Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
                             Dune::FieldVector<double, 1>>;
 
+public:
   double virtual const &C() const = 0;
   double virtual const &L() const = 0;
   double virtual const &V0() const = 0;

dune/tectonic/data-structures/friction/globalratestatefriction.hh

+#ifndef DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
+#define DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
+
+#include <vector>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include "../../spatial-solving/contact/dualmortarcoupling.hh"
+
+#include "globalfrictiondata.hh"
+#include "globalfriction.hh"
+#include "frictioncouplingpair.hh"
+
+#include "../../utils/geocoordinate.hh"
+#include "../../utils/index-in-sorted-range.hh"
+
+template <class Matrix, class Vector, class ScalarFriction, class GridType>
+class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
+public:
+  using GlobalFriction<Matrix, Vector>::block_size;
+  using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;
+
+private:
+  using Base = GlobalFriction<Matrix, Vector>;
+
+  using field_type = typename Vector::field_type;
+  using typename Base::ScalarVector;
+  using typename Base::LocalVectorType;
+
+  using FrictionCoupling = FrictionCouplingPair<GridType, LocalVectorType, field_type>;
+  using ContactCoupling =  DualMortarCoupling<field_type, GridType>;
+
+  size_t bodyIndex(const size_t globalIdx) {
+     size_t i = offSet_.size()-1;
+
+     for (; i>0; ) {
+         if (globalIdx >= offSet_[i])
+             break;
+         i--;
+     }
+     return i;
+  }
+
+public:
+  GlobalRateStateFriction(const std::vector<std::shared_ptr<ContactCoupling>>& contactCouplings, // contains nonmortarBoundary
+                          const std::vector<std::shared_ptr<FrictionCoupling>>& couplings, // contains frictionInfo
+                          const std::vector<ScalarVector>& weights,
+                          const std::vector<ScalarVector>& weightedNormalStress)
+      : restrictions_(), localToGlobal_(), zeroFriction_() {
+
+    assert(contactCouplings.size() == couplings.size());
+    assert(weights.size() == weightedNormalStress.size());
+
+    const auto nBodies = weights.size();
+    offSet_.resize(nBodies, 0);
+    for (size_t i=1; i<nBodies; i++) {
+        offSet_[i] = offSet_[i-1] + weights[i-1].size();
+    }
+
+    std::vector<std::vector<int>> nonmortarBodies(nBodies); // first index body, second index coupling
+    for (size_t i=0; i<contactCouplings.size(); i++) {
+        const auto nonmortarIdx = couplings[i]->gridIdx_[0];
+        nonmortarBodies[nonmortarIdx].emplace_back(i);
+    }
+
+    for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
+        const auto& couplingIndices = nonmortarBodies[bodyIdx];
+
+        if (couplingIndices.size()==0)
+            continue;
+
+        const auto gridView = contactCouplings[couplingIndices[0]]->nonmortarBoundary().gridView();
+
+        Dune::MultipleCodimMultipleGeomTypeMapper<
+            decltype(gridView), Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
+
+        for (auto it = gridView.template begin<block_size>(); it != gridView.template end<block_size>(); ++it) {
+            const auto i = vertexMapper.index(*it);
+
+            for (size_t j=0; j<couplingIndices.size(); j++) {
+                const auto couplingIdx = couplingIndices[j];
+
+                if (not contactCouplings[couplingIdx]->nonmortarBoundary().containsVertex(i))
+                  continue;
+
+                localToGlobal_.emplace_back(offSet_[bodyIdx] + i);
+                restrictions_.emplace_back(weights[bodyIdx][i], weightedNormalStress[bodyIdx][i],
+                                          couplings[j]->frictionData()(geoToPoint(it->geometry())));
+                break;
+            }
+        }
+    }
+  }
+
+  void updateAlpha(const std::vector<ScalarVector>& alpha) override {
+    //print(alpha, "alpha");
+    for (size_t j = 0; j < restrictions_.size(); ++j) {
+      const auto globalDof = localToGlobal_[j];
+      const auto bodyIdx = bodyIndex(globalDof);
+      size_t bodyDof = globalDof - offSet_[bodyIdx];
+
+      restrictions_[j].updateAlpha(alpha[bodyIdx][bodyDof]);
+    }
+  }
+
+  /*
+    Return a restriction of the outer function to the i'th node.
+  */
+  LocalNonlinearity const &restriction(size_t i) const override {
+    auto const index = indexInSortedRange(localToGlobal_, i);
+    if (index == localToGlobal_.size())
+      return zeroFriction_;
+    return restrictions_[index];
+  }
+
+private:
+  std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions_;
+  std::vector<size_t> offSet_; // index off-set by body id
+  std::vector<size_t> localToGlobal_;
+  WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction_;
+};
+#endif

dune/tectonic/localfriction.hh → dune/tectonic/data-structures/friction/localfriction.hh

@@ -10,7 +10,7 @@
 #include <dune/fufem/arithmetic.hh>
 #include <dune/solvers/common/interval.hh>
 
-#include <dune/tectonic/frictionpotential.hh>
+#include "frictionpotential.hh"
 
 template <size_t dimension> class LocalFriction {
 public:
@@ -19,6 +19,7 @@ template <size_t dimension> class LocalFriction {
   using VectorType = Dune::FieldVector<double, dimension>;
   using MatrixType = Dune::FieldMatrix<double, dimension, dimension>;
 
+  double virtual operator()(VectorType const &x) const = 0;
   void virtual updateAlpha(double alpha) = 0;
   double virtual regularity(VectorType const &x) const = 0;
   double virtual coefficientOfFriction(VectorType const &x) const = 0;
@@ -39,34 +40,51 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
   using VectorType = typename LocalFriction<dimension>::VectorType;
   using MatrixType = typename LocalFriction<dimension>::MatrixType;
 
+private:
+    VectorType removeNormal(const VectorType& x) {
+        VectorType res = x;
+        res[0] = 0.0;
+        return res;
+    }
+
 public:
   template <typename... Args>
   WrappedScalarFriction(Args... args)
       : func_(args...) {}
 
+  double operator()(VectorType const &x) const override {
+    auto tangential_x = removeNormal(x);
+    return func_.evaluate(tangential_x.two_norm());
+  }
+
   void updateAlpha(double alpha) override { func_.updateAlpha(alpha); }
 
   double regularity(VectorType const &x) const override {
-    double const xnorm = x.two_norm();
-    if (xnorm <= 0.0)
-      return std::numeric_limits<double>::infinity();
+    auto tangential_x = removeNormal(x);
+    double xnorm = tangential_x.two_norm();
+    if (xnorm < 0.0) {
+        return std::numeric_limits<double>::infinity();
+    }
 
     return func_.regularity(xnorm);
   }
 
   double coefficientOfFriction(VectorType const &x) const override {
-    return func_.coefficientOfFriction(x.two_norm());
+    auto tangential_x = removeNormal(x);
+    return func_.coefficientOfFriction(tangential_x.two_norm());
   }
 
   // directional subdifferential: at u on the line u + t*v
   // u and v are assumed to be non-zero
   void directionalSubDiff(VectorType const &x, VectorType const &v,
                           Dune::Solvers::Interval<double> &D) const override {
-    double const xnorm = x.two_norm();
+
+    auto tangential_x = removeNormal(x);
+    double const xnorm = tangential_x.two_norm();
     if (xnorm <= 0.0)
       D[0] = D[1] = func_.differential(0.0) * v.two_norm();
     else
-      D[0] = D[1] = func_.differential(xnorm) * (x * v) / xnorm;
+      D[0] = D[1] = func_.differential(xnorm) * (tangential_x * v) / xnorm;
   }
 
   /** Formula for the derivative:
@@ -86,17 +104,26 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
       \f}
   */
   void addHessian(VectorType const &x, MatrixType &A) const override {
-    double const xnorm2 = x.two_norm2();
+    //std::cout << A << std::endl;
+    //std::cout << x << std::endl;
+
+ /*   double const xnorm2 = x.two_norm2();
     double const xnorm = std::sqrt(xnorm2);
     if (xnorm2 <= 0.0)
       return;
 
+    //std::cout << xnorm << " " << xnorm2 << std::endl;
+
     double const H1 = func_.differential(xnorm);
     double const H2 = func_.second_deriv(xnorm);
 
+    //std::cout << H1 << " " << H2 << std::endl;
+
     double const tensorweight = (H2 - H1 / xnorm) / xnorm2;
     double const idweight = H1 / xnorm;
 
+    //std::cout << tensorweight << " " << idweight << std::endl;
+
     for (size_t i = 0; i < dimension; ++i)
       for (size_t j = 0; j < i; ++j) {
         double const entry = tensorweight * x[i] * x[j];
@@ -107,16 +134,47 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
     for (size_t k = 0; k < dimension; ++k) {
       double const entry = tensorweight * x[k] * x[k];
       A[k][k] += entry + idweight;
+    }*/
+    auto tangential_x = removeNormal(x);
+    const double xnorm = tangential_x.two_norm();
+
+    if (std::isinf(func_.regularity(xnorm)) or xnorm<=0.0)
+      return;
+
+    VectorType y = tangential_x;
+    y /= xnorm;
+
+    double const H1 = func_.differential(xnorm);
+    double const H2 = func_.second_deriv(xnorm);
+
+    double const tensorweight = (H2 - H1 / xnorm);
+    double const idweight = H1 / xnorm;
+
+    //std::cout << tensorweight << " " << idweight << std::endl;
+
+    for (size_t i = 0; i < dimension; ++i)
+      for (size_t j = 0; j < i; ++j) {
+        double const entry = tensorweight * y[i] * y[j];
+        A[i][j] += entry;
+        A[j][i] += entry;
+      }
+
+    for (size_t k = 0; k < dimension; ++k) {
+      double const entry = tensorweight * y[k] * y[k];
+      A[k][k] += entry + idweight;
     }
+
+    //std::cout << A << std::endl;
   }
 
   void addGradient(VectorType const &x, VectorType &y) const override {
-    double const xnorm = x.two_norm();
+    auto tangential_x = removeNormal(x);
+    double const xnorm = tangential_x.two_norm();
     if (std::isinf(func_.regularity(xnorm)))
       return;
 
     if (xnorm > 0.0)
-      Arithmetic::addProduct(y, func_.differential(xnorm) / xnorm, x);
+      Arithmetic::addProduct(y, func_.differential(xnorm) / xnorm, tangential_x);
   }
 
   void directionalDomain(VectorType const &, VectorType const &,

dune/tectonic/data-structures/globalfrictioncontainer.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+/*
+#include "globalfrictioncontainer.hh"
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::operator[](std::size_t i)
+-> IndexObject& {
+    return globalFriction_[i];
+}
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::operator[](std::size_t i) const
+-> const IndexObject& {
+    return globalFriction_[i];
+}
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::size() const
+-> size_t {
+    return globalFriction_.size();
+}
+
+template <class BaseGlobalFriction, size_t depth>
+void GlobalFrictionContainer<BaseGlobalFriction, depth>::resize(std::list<size_t> list) {
+    assert(list.size() <= depth);
+
+    if (list.size() == 0) {
+        globalFriction_.resize(0);
+    } else {
+        globalFriction_.resize(list.front());
+        list.pop_front();
+
+        for (size_t i=0; i<size(); i++) {
+            globalFriction_[i].resize(list);
+        }
+    }
+}
+
+template <class BaseGlobalFriction, size_t depth>
+template <class VectorContainer>
+void GlobalFrictionContainer<BaseGlobalFriction, depth>::updateAlpha(const VectorContainer& newAlpha) {
+    assert(newAlpha.size() == size());
+
+    for (size_t i=0; i<size(); i++) {
+        globalFriction_[i].updateAlpha(newAlpha[i]);
+    }
+}
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::globalFriction()
+-> GlobalFriction& {
+    return globalFriction_;
+}
+
+template <class BaseGlobalFriction, size_t depth>
+auto GlobalFrictionContainer<BaseGlobalFriction, depth>::globalFriction() const
+-> const GlobalFriction& {
+    return globalFriction_;
+}
+
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::operator[](std::size_t i)
+-> IndexObject& {
+    return globalFriction_[i];
+}
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::operator[](std::size_t i) const
+-> const IndexObject& {
+    return globalFriction_[i];
+}
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::size() const
+-> size_t {
+    return globalFriction_.size();
+}
+
+template <class BaseGlobalFriction>
+void GlobalFrictionContainer<BaseGlobalFriction, 1>::resize(std::list<size_t> newSize) {
+    if (newSize.size() > 0) {
+        globalFriction_.resize(newSize.front(), nullptr);
+    } else {
+        globalFriction_.resize(0);
+    }
+}
+
+template <class BaseGlobalFriction>
+template <class Vector>
+void GlobalFrictionContainer<BaseGlobalFriction, 1>::updateAlpha(const Vector& newAlpha) {
+    for (size_t i=0; i<size(); i++) {
+        globalFriction_[i]->updateAlpha(newAlpha);
+    }
+}
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::globalFriction()
+-> GlobalFriction& {
+    return globalFriction_;
+}
+
+template <class BaseGlobalFriction>
+auto GlobalFrictionContainer<BaseGlobalFriction, 1>::globalFriction() const
+-> const GlobalFriction& {
+    return globalFriction_;
+}
+
+
+#include "globalfrictioncontainer_tmpl.cc"
+*/