transformedGlobalRateStateFriction v0

dune/tectonic/CMakeLists.txt

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

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

src/tests/gridgluefrictiontest.cc

@@ -32,6 +32,8 @@
 
 #include "../utils/debugutils.hh"
 
+#include <dune/tectonic/transformedglobalratestatefriction.hh>
+
 #include "common.hh"
 
 const int dim = 2;