From 14a40c3a0a6fa158277f9f363825b76bda78342c Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Sat, 1 Jun 2013 22:12:55 +0200
Subject: [PATCH] Update the state
As a consequence, the nonlinearity is only assembled once
---
dune/tectonic/frictionpotential.hh | 28 ++++++++++++------------
dune/tectonic/globalnonlinearity.hh | 5 +++++
dune/tectonic/globalruinanonlinearity.hh | 12 ++++++----
dune/tectonic/localfriction.hh | 5 +++++
src/assemblers.cc | 5 ++---
src/assemblers.hh | 3 +--
src/assemblers_tmpl.cc | 3 +--
src/one-body-sample.cc | 6 ++---
8 files changed, 39 insertions(+), 28 deletions(-)
diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh
index e2a8c87f..c88afb32 100644
--- a/dune/tectonic/frictionpotential.hh
+++ b/dune/tectonic/frictionpotential.hh
@@ -27,26 +27,15 @@ class FrictionPotentialWrapper {
// Between 0 and this point, the function is constantly zero (and
// thus so are all derivatives)
double virtual smallestPositivePoint() const = 0;
+ void virtual updateState(double) = 0;
};
// phi(V) = V log(V/V_m) - V + V_m if V >= V_m
// = 0 otherwise
-// with V_m = V_0 exp(-K/a),
-// i.e. K = -a log(V_m / V_0) = mu_0 + b log(V_0 / L) + b alpha
class FrictionPotential : public FrictionPotentialWrapper {
public:
- FrictionPotential(double coefficient, FrictionData const &fd, double state)
- : FrictionPotentialWrapper(),
- coefficientProduct(coefficient * fd.a * fd.normalStress),
- // state is assumed to be logarithmic
- V_m(fd.V0 *
- std::exp(-(fd.mu0 + fd.b * (state + std::log(fd.V0 / fd.L))) /
- fd.a))
- // We could also compute V_m as
- // V_0 * std::exp(-(mu_0 + b * state)/a)
- // * std::pow(V_0 / L, -b/a)
- // which would avoid the std::exp(std::log())
- {}
+ FrictionPotential(double coefficient, FrictionData const &fd)
+ : fd(fd), coefficientProduct(coefficient * fd.a * fd.normalStress) {}
double evaluate(double V) const {
assert(V >= 0);
@@ -89,7 +78,16 @@ class FrictionPotential : public FrictionPotentialWrapper {
double smallestPositivePoint() const { return V_m; }
+ // V_m = V_0 exp(-K/a),
+ // with K = -a log(V_m / V_0) = mu_0 + b [ alpha + log(V_0 / L) ]
+ void updateState(double state) {
+ // state is assumed to be logarithmic
+ V_m = fd.V0 *
+ std::exp(-(fd.mu0 + fd.b * (state + std::log(fd.V0 / fd.L))) / fd.a);
+ }
+
private:
+ FrictionData const &fd;
double const coefficientProduct;
double V_m;
};
@@ -107,6 +105,8 @@ class TrivialFunction : public FrictionPotentialWrapper {
double smallestPositivePoint() const {
return std::numeric_limits<double>::infinity();
}
+
+ void updateState(double) {}
};
}
#endif
diff --git a/dune/tectonic/globalnonlinearity.hh b/dune/tectonic/globalnonlinearity.hh
index 246cc7f1..0e1b71e0 100644
--- a/dune/tectonic/globalnonlinearity.hh
+++ b/dune/tectonic/globalnonlinearity.hh
@@ -14,6 +14,9 @@
namespace Dune {
template <class MatrixTypeTEMPLATE, class VectorTypeTEMPLATE>
class GlobalNonlinearity {
+private:
+ using dataref = BlockVector<FieldVector<double, 1>> const &;
+
public:
using MatrixType = MatrixTypeTEMPLATE;
using VectorType = VectorTypeTEMPLATE;
@@ -74,6 +77,8 @@ class GlobalNonlinearity {
auto const res = restriction(i);
return res->regularity(x);
}
+
+ virtual void updateState(dataref state) = 0;
};
}
#endif
diff --git a/dune/tectonic/globalruinanonlinearity.hh b/dune/tectonic/globalruinanonlinearity.hh
index 81b8d902..f57e447c 100644
--- a/dune/tectonic/globalruinanonlinearity.hh
+++ b/dune/tectonic/globalruinanonlinearity.hh
@@ -24,20 +24,24 @@ class GlobalRuinaNonlinearity
using GlobalNonlinearity<MatrixType, VectorType>::dim;
GlobalRuinaNonlinearity(Dune::BitSetVector<1> const &frictionalNodes,
- dataref nodalIntegrals, FrictionData const &fd,
- dataref state)
+ dataref nodalIntegrals, FrictionData const &fd)
: restrictions(nodalIntegrals.size()) {
auto trivialNonlinearity =
make_shared<LocalFriction<dim>>(make_shared<TrivialFunction>());
for (size_t i = 0; i < restrictions.size(); ++i) {
restrictions[i] =
frictionalNodes[i][0]
- ? make_shared<LocalFriction<dim>>(make_shared<FrictionPotential>(
- nodalIntegrals[i], fd, state[i]))
+ ? make_shared<LocalFriction<dim>>(
+ make_shared<FrictionPotential>(nodalIntegrals[i], fd))
: trivialNonlinearity;
}
}
+ void updateState(dataref state) override {
+ for (size_t i = 0; i < restrictions.size(); ++i)
+ restrictions[i]->updateState(state[i]);
+ }
+
/*
Return a restriction of the outer function to the i'th node.
*/
diff --git a/dune/tectonic/localfriction.hh b/dune/tectonic/localfriction.hh
index 01138327..bfac47f6 100644
--- a/dune/tectonic/localfriction.hh
+++ b/dune/tectonic/localfriction.hh
@@ -25,6 +25,11 @@ template <int dimension> class LocalFriction {
return func->evaluate(x.two_norm());
}
+ void updateState(double state) {
+ func->updateState(state);
+ smp = func->smallestPositivePoint();
+ }
+
double regularity(VectorType const &x) const {
double const xnorm = x.two_norm();
if (xnorm < 1e-14 && xnorm >= smp)
diff --git a/src/assemblers.cc b/src/assemblers.cc
index 74f3a4e8..1449acc2 100644
--- a/src/assemblers.cc
+++ b/src/assemblers.cc
@@ -52,11 +52,10 @@ Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType>>
assemble_nonlinearity(
Dune::BitSetVector<1> const &frictionalNodes,
Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
- FrictionData const &fd,
- Dune::BlockVector<Dune::FieldVector<double, 1>> const &state) {
+ FrictionData const &fd) {
return Dune::make_shared<
Dune::GlobalRuinaNonlinearity<MatrixType, VectorType>>(
- frictionalNodes, nodalIntegrals, fd, state);
+ frictionalNodes, nodalIntegrals, fd);
}
#include "assemblers_tmpl.cc"
diff --git a/src/assemblers.hh b/src/assemblers.hh
index 22b19e69..053a336b 100644
--- a/src/assemblers.hh
+++ b/src/assemblers.hh
@@ -28,6 +28,5 @@ Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType>>
assemble_nonlinearity(
Dune::BitSetVector<1> const &frictionalNodes,
Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
- FrictionData const &fd,
- Dune::BlockVector<Dune::FieldVector<double, 1>> const &state);
+ FrictionData const &fd);
#endif
diff --git a/src/assemblers_tmpl.cc b/src/assemblers_tmpl.cc
index 87a44928..325eaee8 100644
--- a/src/assemblers_tmpl.cc
+++ b/src/assemblers_tmpl.cc
@@ -35,5 +35,4 @@ template Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType>>
assemble_nonlinearity<MatrixType, VectorType>(
Dune::BitSetVector<1> const &frictionalNodes,
Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
- FrictionData const &fd,
- Dune::BlockVector<Dune::FieldVector<double, 1>> const &state);
+ FrictionData const &fd);
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 50459522..61bcf846 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -321,6 +321,8 @@ int main(int argc, char *argv[]) {
auto const nodalIntegrals = assemble_frictional<GridView, SmallVector>(
leafView, p1Assembler, frictionalNodes);
+ auto myGlobalNonlinearity = assemble_nonlinearity<MatrixType, VectorType>(
+ frictionalNodes, *nodalIntegrals, frictionData);
// {{{ Initial conditions
VectorType u_initial(finestSize);
@@ -407,9 +409,7 @@ int main(int argc, char *argv[]) {
size_t iterationCounter;
auto solveVelocityProblem = [&](VectorType &_problem_iterate,
SingletonVectorType const &_alpha) {
- auto myGlobalNonlinearity =
- assemble_nonlinearity<MatrixType, VectorType>(
- frictionalNodes, *nodalIntegrals, frictionData, _alpha);
+ myGlobalNonlinearity->updateState(_alpha);
using MyConvexProblemType = MyConvexProblem<MatrixType, VectorType>;
MyConvexProblemType const myConvexProblem(
--
GitLab