[Extend] Give friction an overhaul

* Move coefficient computation into RuinaNonlinearity
* Allow some friction data to vary in space

dune/tectonic/frictiondata.hh

@@ -3,12 +3,8 @@
 #include <dune/common/parametertree.hh>
 
 struct FrictionData {
-  FrictionData(Dune::ParameterTree const &parset)
-      : L(parset.get<double>("L")),
-        V0(parset.get<double>("V0")),
-        a(parset.get<double>("a")),
-        b(parset.get<double>("b")),
-        mu0(parset.get<double>("mu0")) {}
+  explicit FrictionData(double _L, double _V0, double _a, double _b, double _mu0)
+      : L(_L), V0(_V0), a(_a), b(_b), mu0(_mu0) {}
 
   double const L;
   double const V0;

dune/tectonic/frictionpotential.hh

@@ -18,6 +18,7 @@ class FrictionPotentialWrapper {
   double virtual differential(double s) const = 0;
   double virtual second_deriv(double x) const = 0;
   double virtual regularity(double s) const = 0;
+  double virtual coefficientOfFriction(double s) const = 0;
 
   double virtual evaluate(double x) const {
     DUNE_THROW(Dune::NotImplemented, "evaluation not implemented");
@@ -28,16 +29,19 @@ class FrictionPotentialWrapper {
 
 class FrictionPotential : public FrictionPotentialWrapper {
 public:
-  FrictionPotential(double coefficient, double _normalStress,
-                    FrictionData const &_fd)
+  FrictionPotential(double coefficient, double _normalStress, FrictionData _fd)
       : fd(_fd), weight(coefficient), normalStress(_normalStress) {}
 
-  double differential(double V) const {
+  double coefficientOfFriction(double V) const {
     assert(V >= 0.0);
     if (V <= V_cutoff)
       return 0.0;
 
-    return weight * (-normalStress) * fd.a * (std::log(V) - logV_m);
+    return fd.a * (std::log(V) - logV_m);
+  }
+
+  double differential(double V) const {
+    return weight * (-normalStress) * coefficientOfFriction(V);
   }
 
   double second_deriv(double V) const {
@@ -64,7 +68,7 @@ class FrictionPotential : public FrictionPotentialWrapper {
   }
 
 private:
-  FrictionData const &fd;
+  FrictionData const fd;
   double const weight;
   double const normalStress;
   double logV_m;
@@ -75,6 +79,8 @@ class TrivialFunction : public FrictionPotentialWrapper {
 public:
   double evaluate(double) const { return 0; }
 
+  double coefficientOfFriction(double s) const { return 0; }
+
   double differential(double) const { return 0; }
 
   double second_deriv(double) const { return 0; }

dune/tectonic/globalfrictiondata.hh

+#ifndef GLOBAL_FRICTION_DATA_HH
+#define GLOBAL_FRICTION_DATA_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/tectonic/frictiondata.hh>
+
+template <class DomainType>
+double evaluateScalarFunction(
+    Dune::VirtualFunction<DomainType, Dune::FieldVector<double, 1>> const &f,
+    DomainType const &x) {
+  Dune::FieldVector<double, 1> ret;
+  f.evaluate(x, ret);
+  return ret;
+};
+
+template <int dimension> class GlobalFrictionData {
+public:
+  FrictionData operator()(Dune::FieldVector<double, dimension> const &x) const {
+    return FrictionData(L(), V0(), evaluateScalarFunction(a(), x),
+                        evaluateScalarFunction(b(), x), mu0());
+  }
+
+protected:
+  using VirtualFunction = Dune::VirtualFunction<
+      Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
+
+  double virtual const &L() const = 0;
+  double virtual const &V0() const = 0;
+  VirtualFunction virtual const &a() const = 0;
+  VirtualFunction virtual const &b() const = 0;
+  double virtual const &mu0() const = 0;
+};
+
+#endif

dune/tectonic/globalruinanonlinearity.hh

@@ -6,13 +6,16 @@
 #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/tectonic/globalfrictiondata.hh>
+
 #include "globalnonlinearity.hh"
 #include "frictionpotential.hh"
 
-template <class Matrix, class Vector>
+template <class Matrix, class Vector, class GridView>
 class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
 public:
   using GlobalNonlinearity<Matrix, Vector>::block_size;
@@ -22,19 +25,28 @@ class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
   using typename GlobalNonlinearity<Matrix, Vector>::ScalarVector;
 
 public:
-  GlobalRuinaNonlinearity(Dune::BitSetVector<1> const &frictionalNodes,
-                          ScalarVector const &weights, FrictionData const &fd,
+  GlobalRuinaNonlinearity(BoundaryPatch<GridView> const &frictionalBoundary,
+                          GridView const &gridView,
+                          GlobalFrictionData<block_size> const &frictionInfo,
+                          // Note: passing the following two makes no sense
+                          ScalarVector const &weights,
                           ScalarVector const &normalStress)
-      : restrictions(frictionalNodes.size()) {
+      : restrictions(normalStress.size()) {
     auto trivialNonlinearity =
         std::make_shared<Friction>(std::make_shared<TrivialFunction>());
-    for (size_t i = 0; i < restrictions.size(); ++i) {
-      if (not frictionalNodes[i][0]) {
+
+    Dune::MultipleCodimMultipleGeomTypeMapper<
+        GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
+    for (auto it = gridView.template begin<block_size>();
+         it != gridView.template end<block_size>(); ++it) {
+      auto const i = vertexMapper.map(*it);
+      auto const coordinate = it->geometry().corner(0);
+      if (not frictionalBoundary.containsVertex(i)) {
         restrictions[i] = trivialNonlinearity;
         continue;
       }
-      auto const fp =
-          std::make_shared<FrictionPotential>(weights[i], normalStress[i], fd);
+      auto const fp = std::make_shared<FrictionPotential>(
+          weights[i], normalStress[i], frictionInfo(coordinate));
       restrictions[i] = std::make_shared<Friction>(fp);
     }
   }
@@ -44,6 +56,13 @@ class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
       restrictions[i]->updateLogState(logState[i]);
   }
 
+  void coefficientOfFriction(Vector const &x, ScalarVector &coefficient) {
+    assert(x.size() == restrictions.size());
+    coefficient.resize(restrictions.size());
+    for (size_t i = 0; i < restrictions.size(); ++i)
+      coefficient[i] = restrictions[i]->coefficientOfFriction(x[i]);
+  }
+
   /*
     Return a restriction of the outer function to the i'th node.
   */

dune/tectonic/localfriction.hh

@@ -51,6 +51,10 @@ template <size_t dimension> class LocalFriction {
     return func->regularity(xnorm);
   }
 
+  double coefficientOfFriction(VectorType const &x) const {
+    return func->coefficientOfFriction(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,

src/assemblers.cc

@@ -92,9 +92,10 @@ void MyAssembler<GridView, dimension>::assembleNeumann(
 
 template <class GridView, int dimension>
 auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
-    BoundaryPatch<GridView> const &frictionalBoundary, FrictionData const &fd,
+    BoundaryPatch<GridView> const &frictionalBoundary,
+    GlobalFrictionData<dimension> const &frictionInfo,
     ScalarVector const &normalStress)
-    -> std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector>> {
+    -> std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector, GridView>> {
   // Lump negative normal stress (kludge)
   ScalarVector weights;
   {
@@ -105,13 +106,8 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
     vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
                                                weights, frictionalBoundary);
   }
-  for (size_t i = 0; i < weights.size(); ++i)
-    assert(weights[i] >= 0.0);
-
-  Dune::BitSetVector<1> frictionalNodes;
-  frictionalBoundary.getVertices(frictionalNodes);
-  return std::make_shared<GlobalRuinaNonlinearity<Matrix, Vector>>(
-      frictionalNodes, weights, fd, normalStress);
+  return std::make_shared<GlobalRuinaNonlinearity<Matrix, Vector, GridView>>(
+      frictionalBoundary, gridView, frictionInfo, weights, normalStress);
 }
 
 template <class GridView, int dimension>

src/assemblers.hh

@@ -13,6 +13,7 @@
 #pragma clang diagnostic pop
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 
+#include <dune/tectonic/globalfrictiondata.hh>
 #include <dune/tectonic/globalruinanonlinearity.hh>
 
 template <class GridView, int dimension> class MyAssembler {
@@ -66,9 +67,10 @@ template <class GridView, int dimension> class MyAssembler {
                        Dune::VirtualFunction<double, double> const &neumann,
                        double relativeTime);
 
-  std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector>>
+  std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector, GridView>>
   assembleFrictionNonlinearity(
-      BoundaryPatch<GridView> const &frictionalBoundary, FrictionData const &fd,
+      BoundaryPatch<GridView> const &frictionalBoundary,
+      GlobalFrictionData<dimension> const &frictionInfo,
       ScalarVector const &normalStress);
 
   void assembleVonMisesStress(double youngModulus, double poissonRatio,

src/friction_writer.cc

@@ -4,20 +4,12 @@
 
 #include "friction_writer.hh"
 
-double computeCOF(FrictionData const &fd, double V, double logState) {
-  double const mu = fd.mu0 + fd.a * std::log(V / fd.V0) +
-                    fd.b * (logState + std::log(fd.V0 / fd.L));
-  return (mu > 0) ? mu : 0;
-}
-
 template <class BitVector>
-FrictionWriter<BitVector>::FrictionWriter(FrictionData const &_fd,
-                                          BitVector const &_boundaryNodes)
+FrictionWriter<BitVector>::FrictionWriter(BitVector const &_boundaryNodes)
     : coefficientWriter("coefficients", std::fstream::out),
       displacementWriter("displacements", std::fstream::out),
       stateWriter("states", std::fstream::out),
       velocityWriter("velocities", std::fstream::out),
-      fd(_fd),
       boundaryNodes(_boundaryNodes) {}
 
 template <class BitVector> FrictionWriter<BitVector>::~FrictionWriter() {
@@ -29,13 +21,14 @@ template <class BitVector> FrictionWriter<BitVector>::~FrictionWriter() {
 
 template <class BitVector>
 template <class ScalarVector, class Vector>
-void FrictionWriter<BitVector>::writeInfo(ScalarVector const &alpha,
+void FrictionWriter<BitVector>::writeInfo(ScalarVector const &coefficient,
+                                          ScalarVector const &alpha,
                                           Vector const &u, Vector const &v) {
   for (size_t i = 0; i < boundaryNodes.size(); ++i) {
     if (!boundaryNodes[i][0])
       continue;
 
-    coefficientWriter << computeCOF(fd, v[i].two_norm(), alpha[i]) << " ";
+    coefficientWriter << coefficient[i] << " ";
     displacementWriter << u[i][0] << " ";
     stateWriter << alpha[i][0] << " ";
     velocityWriter << v[i][0] << " ";

src/friction_writer.hh

@@ -3,23 +3,21 @@
 
 #include <fstream>
 
-#include <dune/tectonic/frictiondata.hh>
-
 template <class BitVector> class FrictionWriter {
 public:
-  FrictionWriter(FrictionData const &_fd, BitVector const &_boundaryNodes);
+  FrictionWriter(BitVector const &_boundaryNodes);
 
   ~FrictionWriter();
 
   template <class ScalarVector, class Vector>
-  void writeInfo(ScalarVector const &alpha, Vector const &u, Vector const &v);
+  void writeInfo(ScalarVector const &coefficient, ScalarVector const &alpha,
+                 Vector const &u, Vector const &v);
 
 private:
   std::fstream coefficientWriter;
   std::fstream displacementWriter;
   std::fstream stateWriter;
   std::fstream velocityWriter;
-  FrictionData const &fd;
   BitVector const &boundaryNodes;
 };
 #endif

src/friction_writer_tmpl.cc

@@ -10,6 +10,6 @@ using BitVector = Dune::BitSetVector<1>;
 
 template class FrictionWriter<BitVector>;
 
-template void FrictionWriter<BitVector>::writeInfo(ScalarVector const &alpha,
-                                                   Vector const &u,
-                                                   Vector const &v);
+template void FrictionWriter<BitVector>::writeInfo(
+    ScalarVector const &coefficient, ScalarVector const &alpha, Vector const &u,
+    Vector const &v);

src/mybody.hh

@@ -13,7 +13,6 @@
 template <int dimension> class MyBody : public Body<dimension> {
   using typename Body<dimension>::ScalarFunction;
   using typename Body<dimension>::VectorField;
-
   using MyConstantFunction = ConstantFunction<
       Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
 

src/myglobalfrictiondata.hh

+#ifndef MY_GLOBAL_FRICTION_DATA_HH
+#define MY_GLOBAL_FRICTION_DATA_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/functions/constantfunction.hh>
+
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "mygeometry.hh"
+
+template <int dimension>
+class MyGlobalFrictionData : public GlobalFrictionData<dimension> {
+private:
+  using typename GlobalFrictionData<dimension>::VirtualFunction;
+
+public:
+  MyGlobalFrictionData(Dune::ParameterTree const &parset, MyGeometry const &tri)
+      : L_(parset.get<double>("L")),
+        V0_(parset.get<double>("V0")),
+        a_(parset.get<double>("a")),
+        b_(parset.get<double>("b")),
+        mu0_(parset.get<double>("mu0")) {}
+
+  double virtual const &L() const override { return L_; }
+  double virtual const &V0() const override { return V0_; }
+  VirtualFunction virtual const &a() const override { return a_; }
+  VirtualFunction virtual const &b() const override { return b_; }
+  double virtual const &mu0() const override { return mu0_; }
+
+private:
+  using MyConstantFunction = ConstantFunction<
+      Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
+
+  double const L_;
+  double const V0_;
+  MyConstantFunction const a_;
+  MyConstantFunction const b_;
+  double const mu0_;
+};
+
+#endif

src/one-body-sample.cc

@@ -76,6 +76,7 @@
 #include "friction_writer.hh"
 #include "mybody.hh"
 #include "mygeometry.hh"
+#include "myglobalfrictiondata.hh"
 #include "mygrid.hh"
 #include "solverfactory.hh"
 #include "state.hh"
@@ -193,14 +194,15 @@ int main(int argc, char *argv[]) {
     normalStress = (myGeometry.A[1] - myGeometry.C[1]) *
                    parset.get<double>("gravity") *
                    parset.get<double>("body.density");
-    FrictionData const frictionData(parset.sub("boundary.friction"));
+    MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"),
+                                            myGeometry);
 
     // {{{ Initial conditions
     ScalarVector alpha_initial(fineVertexCount);
     alpha_initial =
         std::log(parset.get<double>("boundary.friction.initialState"));
     auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
-        frictionalBoundary, frictionData, normalStress);
+        frictionalBoundary, frictionInfo, normalStress);
     myGlobalNonlinearity->updateLogState(alpha_initial);
 
     using LinearFactory = SolverFactory<
@@ -267,8 +269,12 @@ int main(int argc, char *argv[]) {
                          parset.sub("a0.solver"));
     }
     // }}}
-    FrictionWriter<Dune::BitSetVector<1>> writer(frictionData, frictionalNodes);
-    writer.writeInfo(alpha_initial, u_initial, v_initial);
+    FrictionWriter<Dune::BitSetVector<1>> writer(frictionalNodes);
+    {
+      ScalarVector c;
+      myGlobalNonlinearity->coefficientOfFriction(v_initial, c);
+      writer.writeInfo(c, alpha_initial, u_initial, v_initial);
+    }
 
     MyVTKWriter<typename MyAssembler::VertexBasis,
                 typename MyAssembler::CellBasis> const
@@ -310,7 +316,8 @@ int main(int argc, char *argv[]) {
                         u_initial, v_initial, a_initial);
     auto stateUpdater = initStateUpdater<ScalarVector, Vector>(
         parset.get<Config::stateModel>("boundary.friction.stateModel"),
-        alpha_initial, frictionalNodes, frictionData);
+        alpha_initial, frictionalNodes,
+        parset.get<double>("boundary.friction.L"));
 
     Vector v = v_initial;
     ScalarVector alpha(fineVertexCount);
@@ -418,7 +425,11 @@ int main(int argc, char *argv[]) {
       if (printProgress)
         std::cout << std::endl;
 
-      writer.writeInfo(alpha, u, v);
+      {
+        ScalarVector c;
+        myGlobalNonlinearity->coefficientOfFriction(v, c);
+        writer.writeInfo(c, alpha, u, v);
+      }
       iterationWriter << std::endl;
       relaxationWriter << std::endl;
 

src/state.hh

@@ -11,14 +11,14 @@
 template <class ScalarVector, class Vector>
 std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater(
     Config::stateModel model, ScalarVector const &alpha_initial,
-    Dune::BitSetVector<1> const &frictionalNodes, FrictionData const &fd) {
+    Dune::BitSetVector<1> const &frictionalNodes, double L) {
   switch (model) {
     case Config::Dieterich:
       return std::make_shared<DieterichStateUpdater<ScalarVector, Vector>>(
-          alpha_initial, frictionalNodes, fd.L);
+          alpha_initial, frictionalNodes, L);
     case Config::Ruina:
       return std::make_shared<RuinaStateUpdater<ScalarVector, Vector>>(
-          alpha_initial, frictionalNodes, fd.L);
+          alpha_initial, frictionalNodes, L);
     default:
       assert(false);
   }