[Bug fix] Fix normal stress computation

* Use PointTractionBoundaryAssembler: The normal stress is properly calculated here
* Incorporate weights into normal stress

dune/tectonic/pointtractionboundaryassembler.hh

+#ifndef DUNE_TECTONIC_POINTTRACTIONBOUNDARYASSEMBLER_HH
+#define DUNE_TECTONIC_POINTTRACTIONBOUNDARYASSEMBLER_HH
+
+// Based on
+// dune/fufem/assemblers/localassemblers/normalstressboundaryassembler.hh
+
+#include <dune/common/fvector.hh>
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/fufem/assemblers/localboundaryassembler.hh>
+#include <dune/fufem/symmetrictensor.hh>
+#include <dune/fufem/functions/virtualgridfunction.hh>
+#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
+
+/** \brief Computes the surface traction of a linear elastic material,
+    pointwise. This is done in such a manner that the sum over all
+    vertices equals the (approximate) integral */
+template <class GridType>
+class PointTractionBoundaryAssembler
+    : public LocalBoundaryAssembler<
+          GridType,
+          Dune::FieldVector<typename GridType::ctype, GridType::dimension>> {
+  static const int dim = GridType::dimension;
+  typedef typename GridType::ctype ctype;
+  typedef typename Dune::FieldVector<ctype, dim> T;
+
+public:
+  typedef typename LocalBoundaryAssembler<GridType, T>::LocalVector LocalVector;
+  typedef VirtualGridFunction<
+      GridType, Dune::FieldVector<ctype, GridType::dimension>> GridFunction;
+
+  //! Create assembler for grid
+  PointTractionBoundaryAssembler(double E, double nu,
+                                 const GridFunction* displacement, int order)
+      : E_(E), nu_(nu), displacement_(displacement), order_(order) {}
+
+  /** \brief Assemble the normal stress at a boundary face.*/
+  template <class TrialLocalFE, class BoundaryIterator>
+  void assemble(const BoundaryIterator& it, LocalVector& localVector,
+                const TrialLocalFE& tFE) {
+
+    localVector = 0.0;
+
+    // geometry of the boundary face
+    const typename BoundaryIterator::Intersection::Geometry segmentGeometry =
+        it->geometry();
+
+    // get quadrature rule
+    const int order = (it->inside()->type().isSimplex())
+                          ? 2 * (order_ - 1)
+                          : 2 * (order_ * dim - 1);
+
+    // get quadrature rule
+    const Dune::QuadratureRule<ctype, dim - 1>& quad =
+        QuadratureRuleCache<ctype, dim - 1>::rule(segmentGeometry.type(), order,
+                                                  false);
+
+    // loop over quadrature points
+    for (size_t pt = 0; pt < quad.size(); ++pt) {
+      // get quadrature point
+      const Dune::FieldVector<ctype, dim - 1>& quadPos = quad[pt].position();
+
+      // get integration factor
+      const ctype integrationElement =
+          segmentGeometry.integrationElement(quadPos);
+
+      // get outer unit normal at quad pos
+      Dune::FieldVector<ctype, dim> normalAtQuadPos =
+          it->unitOuterNormal(quadPos);
+
+      // position of the quadrature point within the element
+      const Dune::FieldVector<ctype, dim> elementQuadPos =
+          it->geometryInInside().global(quadPos);
+
+      // evaluate the displacement gradient at quad point of the element
+      typename GridFunction::DerivativeType localDispGradient;
+      if (displacement_->isDefinedOn(*it->inside()))
+        displacement_->evaluateDerivativeLocal(*it->inside(), elementQuadPos,
+                                               localDispGradient);
+      else
+        displacement_->evaluateDerivative(segmentGeometry.global(quadPos),
+                                          localDispGradient);
+
+      // Compute strain tensor from the displacement gradient
+      SymmetricTensor<dim> strain;
+      computeStrain(localDispGradient, strain);
+
+      // and the stress
+      SymmetricTensor<dim> stress = hookeTimesStrain(strain);
+
+      normalAtQuadPos *= quad[pt].weight() * integrationElement;
+      typename LocalVector::block_type tractionAtQuadPos;
+      stress.umv(normalAtQuadPos, tractionAtQuadPos);
+
+      // Divide the traction
+      tractionAtQuadPos /= segmentGeometry.corners();
+
+      for (size_t k = 0; k < localVector.size(); k++)
+        localVector[k] += tractionAtQuadPos;
+    }
+  }
+
+  ~PointTractionBoundaryAssembler() {}
+
+private:
+  /** \brief Young's modulus */
+  double E_;
+
+  /** \brief The Poisson ratio */
+  double nu_;
+
+  /** \brief The displacement.*/
+  const GridFunction* displacement_;
+
+  /** \brief The gridfunction is usually corresponding to some finite element
+   * solution. This is its order.*/
+  int order_;
+
+  /** \brief Compute linear strain tensor from the deformation gradient. */
+  void computeStrain(const Dune::FieldMatrix<double, dim, dim>& grad,
+                     SymmetricTensor<dim>& strain) const {
+    for (int i = 0; i < dim; ++i) {
+      strain(i, i) = grad[i][i];
+      for (int j = i + 1; j < dim; ++j)
+        strain(i, j) = 0.5 * (grad[i][j] + grad[j][i]);
+    }
+  }
+
+  /** \brief Compute the stress tensor from strain tensor. */
+  SymmetricTensor<dim> hookeTimesStrain(
+      const SymmetricTensor<dim>& strain) const {
+
+    SymmetricTensor<dim> stress = strain;
+    stress *= E_ / (1 + nu_);
+
+    double f = E_ * nu_ / ((1 + nu_) * (1 - 2 * nu_)) * strain.trace();
+
+    for (int i = 0; i < dim; i++)
+      stress[i] += f;
+
+    return stress;
+  }
+};
+
+#endif

src/assemblers.cc

@@ -12,13 +12,13 @@
 #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/weightedmassassembler.hh>
 #include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/computestress.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
 #include <dune/fufem/quadraturerules/quadraturerulecache.hh>
 
 #include <dune/tectonic/frictionpotential.hh>
 #include <dune/tectonic/globalratestatefriction.hh>
+#include <dune/tectonic/pointtractionboundaryassembler.hh>
 
 #include "assemblers.hh"
 
@@ -104,21 +104,27 @@ void MyAssembler<GridView, dimension>::assembleNeumann(
 }
 
 template <class GridView, int dimension>
-void MyAssembler<GridView, dimension>::assembleNormalStress(
+void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
     BoundaryPatch<GridView> const &frictionalBoundary,
-    ScalarVector &normalStress, double youngModulus, double poissonRatio,
-    Vector const &displacement) const {
+    ScalarVector &weightedNormalStress, double youngModulus,
+    double poissonRatio, Vector const &displacement) const {
 
+  BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
+                                                              displacement);
   Vector traction;
-  Stress<GridView>::getElasticSurfaceNormalStress // misnomer(!)
-      (frictionalBoundary, displacement, traction, youngModulus, poissonRatio);
+  // Note: We assume v = 0 here so that there is no viscous
+  // contribution to the stress.
+  PointTractionBoundaryAssembler<Grid> weightedNormalStressBoundaryAssembler(
+      youngModulus, poissonRatio, &displacementFunction, 1);
+  vertexAssembler.assembleBoundaryFunctional(
+      weightedNormalStressBoundaryAssembler, traction, frictionalBoundary);
 
   std::vector<typename Vector::block_type> normals;
   frictionalBoundary.getNormals(normals);
   for (size_t i = 0; i < traction.size(); ++i) {
-    normalStress[i] = normals[i] * traction[i];
-    if (normalStress[i] > 0.0) {
-      normalStress[i] = 0.0;
+    weightedNormalStress[i] = normals[i] * traction[i];
+    if (weightedNormalStress[i] > 0.0) {
+      weightedNormalStress[i] = 0.0;
       std::cout << "Warning: Manually reducing positive normal stress to zero."
                 << std::endl;
     }
@@ -130,9 +136,9 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
     Config::FrictionModel frictionModel,
     BoundaryPatch<GridView> const &frictionalBoundary,
     GlobalFrictionData<dimension> const &frictionInfo,
-    ScalarVector const &normalStress) const
+    ScalarVector const &weightedNormalStress) const
     -> std::shared_ptr<GlobalFriction<Matrix, Vector>> {
-  // Lump negative normal stress (kludge)
+  // Lumping of the nonlinearity
   ScalarVector weights;
   {
     NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
@@ -142,10 +148,6 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
     vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
                                                weights, frictionalBoundary);
   }
-  ScalarVector weightedNormalStress(weights.size());
-  for (size_t i = 0; i < weights.size(); ++i)
-    weightedNormalStress[i] = weights[i] * normalStress[i];
-
   switch (frictionModel) {
     case Config::Truncated:
       return std::make_shared<GlobalRateStateFriction<

src/assemblers.hh

@@ -73,16 +73,16 @@ template <class GridView, int dimension> class MyAssembler {
                        Dune::VirtualFunction<double, double> const &neumann,
                        double relativeTime) const;
 
-  void assembleNormalStress(BoundaryPatch<GridView> const &frictionalBoundary,
-                            ScalarVector &normalStress, double youngModulus,
-                            double poissonRatio,
-                            Vector const &displacement) const;
+  void assembleWeightedNormalStress(
+      BoundaryPatch<GridView> const &frictionalBoundary,
+      ScalarVector &weightedNormalStress, double youngModulus,
+      double poissonRatio, Vector const &displacement) const;
 
   std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity(
       Config::FrictionModel frictionModel,
       BoundaryPatch<GridView> const &frictionalBoundary,
       GlobalFrictionData<dimension> const &frictionInfo,
-      ScalarVector const &normalStress) const;
+      ScalarVector const &weightedNormalStress) const;
 
   void assembleVonMisesStress(double youngModulus, double poissonRatio,
                               Vector const &u, ScalarVector &stress) const;

src/program_state.hh

@@ -19,7 +19,7 @@ template <class Vector, class ScalarVector> class ProgramState {
         v(leafVertexCount),
         a(leafVertexCount),
         alpha(leafVertexCount),
-        normalStress(leafVertexCount) {}
+        weightedNormalStress(leafVertexCount) {}
 
   // Set up initial conditions
   template <class Matrix, class GridView>
@@ -96,9 +96,9 @@ template <class Vector, class ScalarVector> class ProgramState {
     alpha = parset.get<double>("boundary.friction.initialAlpha");
 
     // Initial normal stress
-    myAssembler.assembleNormalStress(frictionalBoundary, normalStress,
-                                     body.getYoungModulus(),
-                                     body.getPoissonRatio(), u);
+    myAssembler.assembleWeightedNormalStress(
+        frictionalBoundary, weightedNormalStress, body.getYoungModulus(),
+        body.getPoissonRatio(), u);
   }
 
 public:
@@ -106,7 +106,7 @@ template <class Vector, class ScalarVector> class ProgramState {
   Vector v;
   Vector a;
   ScalarVector alpha;
-  ScalarVector normalStress;
+  ScalarVector weightedNormalStress;
   double relativeTime;
   double relativeTau;
   size_t timeStep;
@@ -121,7 +121,7 @@ namespace serialization {
     ar &ps.v;
     ar &ps.a;
     ar &ps.alpha;
-    ar &ps.normalStress;
+    ar &ps.weightedNormalStress;
     ar &ps.relativeTime;
     ar &ps.relativeTau;
     ar &ps.timeStep;

src/sand-wedge.cc

@@ -225,7 +225,7 @@ int main(int argc, char *argv[]) {
         parset.sub("boundary.friction"), weakPatch);
     auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
         parset.get<Config::FrictionModel>("boundary.friction.frictionModel"),
-        frictionalBoundary, frictionInfo, programState.normalStress);
+        frictionalBoundary, frictionInfo, programState.weightedNormalStress);
     myGlobalFriction->updateAlpha(programState.alpha);
 
     Vector vertexCoordinates(leafVertexCount);