[Cleanup] Use NormalStressBoundaryAssembler again

This works because interpolateP0ToP1 is in dune-fufem now

dune-tectonic.pc.in

@@ -10,6 +10,6 @@ Name: @PACKAGE_NAME@
 Version: @VERSION@
 Description: dune-tectonic module
 URL: http://dune-project.org/
-Requires: dune-common dune-fufem dune-geometry dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Requires: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
 Libs: -L${libdir}
 Cflags: -I${includedir}

dune.module

@@ -5,4 +5,4 @@
 Module: dune-tectonic
 Version: 2.5-1
 Maintainer: elias.pipping@fu-berlin.de
-Depends: dune-common dune-fufem dune-geometry dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Depends: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid

dune/tectonic/CMakeLists.txt

@@ -10,7 +10,6 @@ install(FILES
   minimisation.hh
   myblockproblem.hh
   mydirectionalconvexfunction.hh
-  pointtractionboundaryassembler.hh
   quadraticenergy.hh
   tectonic.hh
   DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

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) {
-    auto const inside = it->inside();
-    localVector = 0.0;
-
-    // geometry of the boundary face
-    const typename BoundaryIterator::Intersection::Geometry segmentGeometry =
-        it->geometry();
-
-    // get quadrature rule
-    const int order =
-        (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(inside))
-        displacement_->evaluateDerivativeLocal(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

@@ -7,6 +7,7 @@
 #include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/normalstressboundaryassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/variablecoefficientviscosityassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
@@ -14,11 +15,11 @@
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
+#include <dune/fufem/functiontools/p0p1interpolation.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"
 
@@ -111,24 +112,28 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
 
   BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
                                                               displacement);
-  Vector traction;
-  // Note: We assume v = 0 here so that there is no viscous
-  // contribution to the stress.
-  PointTractionBoundaryAssembler<Grid> weightedTractionBoundaryAssembler(
+  Vector traction(cellBasis.size());
+  NormalStressBoundaryAssembler<Grid> tractionBoundaryAssembler(
       youngModulus, poissonRatio, &displacementFunction, 1);
-  vertexAssembler.assembleBoundaryFunctional(weightedTractionBoundaryAssembler,
-                                             traction, frictionalBoundary);
-
-  std::vector<typename Vector::block_type> normals;
-  frictionalBoundary.getNormals(normals);
-  for (size_t i = 0; i < traction.size(); ++i) {
-    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;
-    }
+  cellAssembler.assembleBoundaryFunctional(tractionBoundaryAssembler, traction,
+                                           frictionalBoundary);
+
+  auto const nodalTractionAverage =
+      interpolateP0ToP1(frictionalBoundary, traction);
+
+  ScalarVector weights;
+  {
+    NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
+        frictionalBoundaryAssembler(
+            std::make_shared<ConstantFunction<
+                LocalVector, typename ScalarVector::block_type>>(1));
+    vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
+                                               weights, frictionalBoundary);
   }
+  auto const normals = frictionalBoundary.getNormals();
+  for (size_t i = 0; i < vertexBasis.size(); ++i)
+    weightedNormalStress[i] =
+        std::fmin(normals[i] * nodalTractionAverage[i], 0) * weights[i];
 }
 
 template <class GridView, int dimension>