[Extend] Allow normal stress to vary in space

dune/tectonic/frictiondata.hh

@@ -3,19 +3,17 @@
 #include <dune/common/parametertree.hh>
 
 struct FrictionData {
-  FrictionData(Dune::ParameterTree const &parset, double normalStress)
+  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")),
-        normalStress(normalStress) {}
+        mu0(parset.get<double>("mu0")) {}
 
   double const L;
   double const V0;
   double const a;
   double const b;
   double const mu0;
-  double const normalStress;
 };
 #endif

dune/tectonic/frictionpotential.hh

@@ -28,15 +28,16 @@ class FrictionPotentialWrapper {
 
 class FrictionPotential : public FrictionPotentialWrapper {
 public:
-  FrictionPotential(double coefficient, FrictionData const &_fd)
-      : fd(_fd), weightTimesNormalStress(coefficient * (-fd.normalStress)) {}
+  FrictionPotential(double coefficient, double _normalStress,
+                    FrictionData const &_fd)
+      : fd(_fd), weight(coefficient), normalStress(_normalStress) {}
 
   double differential(double V) const {
     assert(V >= 0.0);
     if (V <= V_cutoff)
       return 0.0;
 
-    return weightTimesNormalStress * fd.a * (std::log(V) - logV_m);
+    return weight * (-normalStress) * fd.a * (std::log(V) - logV_m);
   }
 
   double second_deriv(double V) const {
@@ -44,7 +45,7 @@ class FrictionPotential : public FrictionPotentialWrapper {
     if (V <= V_cutoff)
       return 0;
 
-    return weightTimesNormalStress * fd.a / V;
+    return weight * (-normalStress) * (fd.a / V);
   }
 
   double regularity(double V) const {
@@ -64,7 +65,8 @@ class FrictionPotential : public FrictionPotentialWrapper {
 
 private:
   FrictionData const &fd;
-  double const weightTimesNormalStress;
+  double const weight;
+  double const normalStress;
   double logV_m;
   double V_cutoff; // velocity at which mu falls short of mumin
 };

dune/tectonic/globalruinanonlinearity.hh

@@ -23,7 +23,8 @@ class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
 
 public:
   GlobalRuinaNonlinearity(Dune::BitSetVector<1> const &frictionalNodes,
-                          ScalarVector const &weights, FrictionData const &fd)
+                          ScalarVector const &weights, FrictionData const &fd,
+                          ScalarVector const &normalStress)
       : restrictions(frictionalNodes.size()) {
     auto trivialNonlinearity =
         std::make_shared<Friction>(std::make_shared<TrivialFunction>());
@@ -32,7 +33,8 @@ class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
         restrictions[i] = trivialNonlinearity;
         continue;
       }
-      auto const fp = std::make_shared<FrictionPotential>(weights[i], fd);
+      auto const fp =
+          std::make_shared<FrictionPotential>(weights[i], normalStress[i], fd);
       restrictions[i] = std::make_shared<Friction>(fp);
     }
   }

src/assemblers.cc

@@ -92,7 +92,8 @@ 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, FrictionData const &fd,
+    ScalarVector const &normalStress)
     -> std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector>> {
   // Lump negative normal stress (kludge)
   ScalarVector weights;
@@ -110,7 +111,7 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
   Dune::BitSetVector<1> frictionalNodes;
   frictionalBoundary.getVertices(frictionalNodes);
   return std::make_shared<GlobalRuinaNonlinearity<Matrix, Vector>>(
-      frictionalNodes, weights, fd);
+      frictionalNodes, weights, fd, normalStress);
 }
 
 template <class GridView, int dimension>

src/assemblers.hh

@@ -68,8 +68,8 @@ template <class GridView, int dimension> class MyAssembler {
 
   std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector>>
   assembleFrictionNonlinearity(
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      FrictionData const &fd);
+      BoundaryPatch<GridView> const &frictionalBoundary, FrictionData const &fd,
+      ScalarVector const &normalStress);
 
   void assembleVonMisesStress(double youngModulus, double poissonRatio,
                               Vector const &u, ScalarVector &stress);

src/one-body-sample.cc

@@ -188,18 +188,19 @@ int main(int argc, char *argv[]) {
     };
     Vector ell(fineVertexCount);
     computeExternalForces(0.0, ell);
-    double normalStress = (myGeometry.A[1] - myGeometry.C[1]) *
-                          parset.get<double>("gravity") *
-                          parset.get<double>("body.density");
-    FrictionData const frictionData(parset.sub("boundary.friction"),
-                                    normalStress);
+
+    ScalarVector normalStress(fineVertexCount);
+    normalStress = (myGeometry.A[1] - myGeometry.C[1]) *
+                   parset.get<double>("gravity") *
+                   parset.get<double>("body.density");
+    FrictionData const frictionData(parset.sub("boundary.friction"));
 
     // {{{ Initial conditions
     ScalarVector alpha_initial(fineVertexCount);
     alpha_initial =
         std::log(parset.get<double>("boundary.friction.initialState"));
     auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
-        frictionalBoundary, frictionData);
+        frictionalBoundary, frictionData, normalStress);
     myGlobalNonlinearity->updateLogState(alpha_initial);
 
     using LinearFactory = SolverFactory<