From 26e4db8c209ff6fdd2301bbe72d007883ba30f9b Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 9 Mar 2015 17:37:34 +0100
Subject: [PATCH] [Bug fix] Fix normal stress computation
* Use PointTractionBoundaryAssembler: The normal stress is properly calculated here
* Incorporate weights into normal stress
---
.../pointtractionboundaryassembler.hh | 146 ++++++++++++++++++
src/assemblers.cc | 32 ++--
src/assemblers.hh | 10 +-
src/program_state.hh | 12 +-
src/sand-wedge.cc | 2 +-
5 files changed, 175 insertions(+), 27 deletions(-)
create mode 100644 dune/tectonic/pointtractionboundaryassembler.hh
diff --git a/dune/tectonic/pointtractionboundaryassembler.hh b/dune/tectonic/pointtractionboundaryassembler.hh
new file mode 100644
index 00000000..648ec7e2
--- /dev/null
+++ b/dune/tectonic/pointtractionboundaryassembler.hh
@@ -0,0 +1,146 @@
+#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
diff --git a/src/assemblers.cc b/src/assemblers.cc
index 1dc59593..7eb63cf8 100644
--- a/src/assemblers.cc
+++ b/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<
diff --git a/src/assemblers.hh b/src/assemblers.hh
index 9fd600fb..b3227a9a 100644
--- a/src/assemblers.hh
+++ b/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;
diff --git a/src/program_state.hh b/src/program_state.hh
index 6ce0eba0..0c544041 100644
--- a/src/program_state.hh
+++ b/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;
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index e4baf5b5..c8795737 100644
--- a/src/sand-wedge.cc
+++ b/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);
--
GitLab