diff --git a/dune-tectonic.pc.in b/dune-tectonic.pc.in
index 160377afb91624a20ad75114a5134c16112ee564..d9c759b84e0fb0de0fa0b8c366b801ea1ac98763 100644
--- a/dune-tectonic.pc.in
+++ b/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}
diff --git a/dune.module b/dune.module
index 45e5a5a5ba81c2db986ebbbcc5c9e9cce1794e80..e5c502803ccabb7cb6bf7be4bac6ac88f2cae0fe 100644
--- a/dune.module
+++ b/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
diff --git a/dune/tectonic/CMakeLists.txt b/dune/tectonic/CMakeLists.txt
index 3c522fe25fe8e01f95d01b18bc1f7a77f15e695b..7e6ac59121e6adcd6e5f28333526e99e862bf3c0 100644
--- a/dune/tectonic/CMakeLists.txt
+++ b/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)
diff --git a/dune/tectonic/pointtractionboundaryassembler.hh b/dune/tectonic/pointtractionboundaryassembler.hh
deleted file mode 100644
index c64b92276c80f158f78af3e68faa168a2129fca1..0000000000000000000000000000000000000000
--- a/dune/tectonic/pointtractionboundaryassembler.hh
+++ /dev/null
@@ -1,145 +0,0 @@
-#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
diff --git a/src/assemblers.cc b/src/assemblers.cc
index a2b6e1ed92eb8a22aa27cc328ac1a9635bde7166..1097625fbd3055c9b7d68720d93b5fbc3952c2c1 100644
--- a/src/assemblers.cc
+++ b/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>