From 27efa025f784c88363680ab08531b264f430b01b Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Sun, 13 Nov 2011 23:00:55 +0100
Subject: [PATCH] Add Ruina's approach
---
dune/tectonic/globalruinanonlinearity.hh | 44 ++++++++++++++++++
dune/tectonic/nicefunction.hh | 58 ++++++++++++++++++++++++
src/one-body-sample.cc | 22 +++++++--
src/one-body-sample.parset | 2 +
4 files changed, 122 insertions(+), 4 deletions(-)
create mode 100644 dune/tectonic/globalruinanonlinearity.hh
diff --git a/dune/tectonic/globalruinanonlinearity.hh b/dune/tectonic/globalruinanonlinearity.hh
new file mode 100644
index 00000000..aed66653
--- /dev/null
+++ b/dune/tectonic/globalruinanonlinearity.hh
@@ -0,0 +1,44 @@
+/* -*- mode:c++; mode:semantic -*- */
+
+#ifndef DUNE_TECTONIC_GLOBAL_RUINA_NONLINEARITY_HH
+#define DUNE_TECTONIC_GLOBAL_RUINA_NONLINEARITY_HH
+
+#include <dune/common/fvector.hh>
+
+#include "localnonlinearity.hh"
+#include "globalnonlinearity.hh"
+
+namespace Dune {
+template <int dim>
+class GlobalRuinaNonlinearity
+ : public Dune::GlobalNonlinearity<dim, Dune::RuinaFunction> {
+public:
+ GlobalRuinaNonlinearity(
+ std::vector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
+ std::vector<double> const &a,
+ std::vector<double> const &coefficientOfFriction,
+ std::vector<double> const &eta, std::vector<double> const &normalStress)
+ : nodalIntegrals(nodalIntegrals),
+ a(a),
+ coefficientOfFriction(coefficientOfFriction),
+ eta(eta),
+ normalStress(normalStress) {}
+
+ /*
+ Return a restriction of the outer function to the i'th node.
+ */
+ RuinaFunction virtual restriction(int i) const {
+ return Dune::RuinaFunction(nodalIntegrals[i][0], a[i],
+ coefficientOfFriction[i], eta[i],
+ normalStress[i]);
+ }
+
+private:
+ std::vector<Dune::FieldVector<double, 1>> nodalIntegrals;
+ std::vector<double> a;
+ std::vector<double> coefficientOfFriction;
+ std::vector<double> eta;
+ std::vector<double> normalStress;
+};
+}
+#endif
diff --git a/dune/tectonic/nicefunction.hh b/dune/tectonic/nicefunction.hh
index a3614c03..cc582278 100644
--- a/dune/tectonic/nicefunction.hh
+++ b/dune/tectonic/nicefunction.hh
@@ -16,6 +16,64 @@ class NiceFunction : public VirtualFunction<double, double> {
virtual double rightDifferential(double s) const = 0;
};
+class RuinaFunction : public NiceFunction {
+public:
+ RuinaFunction() {}
+
+ RuinaFunction(double coefficient, double a, double mu, double eta,
+ double normalStress)
+ : coefficient(coefficient),
+ a(a),
+ mu(mu),
+ rho(exp(-mu / a)),
+ eta(eta),
+ normalStress(normalStress) {}
+
+ /*
+ If mu and sigma_n denote the coefficient of friction and the
+ normal stress, respectively, this function is given by
+
+ 1/eta sigma_n h(max(eta id, rho)) + a rho + mu
+
+ with the constants a and b from Ruina's model and
+
+ h(beta) = beta (a (log beta - 1) + mu)
+
+ as well as
+
+ rho = exp(-mu/a)
+ */
+ void virtual evaluate(double const &x, double &y) const {
+ double const arg = std::max(eta * x, rho);
+ double const h = arg * (a * (std::log(arg) - 1) + mu);
+ y = 1 / eta * normalStress * h + a * rho + mu;
+ y *= coefficient;
+ }
+
+ double virtual leftDifferential(double s) const {
+ if (eta * s <= rho)
+ return 0;
+
+ return normalStress * (a * std::log(s) + mu);
+ }
+
+ double virtual rightDifferential(double s) const {
+ if (eta * s <= rho)
+ return 0;
+
+ return normalStress * (a * std::log(s) + mu);
+ }
+
+private:
+ double coefficient;
+ double a;
+ double mu;
+ double eta;
+ double normalStress;
+
+ double rho;
+};
+
class LinearFunction : public NiceFunction {
public:
LinearFunction() {}
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 3b0b6146..720f80be 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -42,6 +42,7 @@
#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
#include <dune/tectonic/globallaursennonlinearity.hh>
+#include <dune/tectonic/globalruinanonlinearity.hh>
#include <dune/tectonic/myconvexproblem.hh>
#include <dune/tectonic/myblockproblem.hh>
@@ -186,6 +187,8 @@ int main(int argc, char *argv[]) {
typedef MyConvexProblem<OperatorType, VectorType, Dune::LinearFunction>
MyConvexProblemType;
+ // typedef MyConvexProblem<OperatorType, VectorType, Dune::RuinaFunction>
+ // MyConvexProblemType;
typedef MyBlockProblem<MyConvexProblemType> MyBlockProblemType;
GenericNonlinearGS<MyBlockProblemType> nonlinearGSStep;
nonlinearGSStep.ignoreNodes_ = &ignoreNodes;
@@ -216,18 +219,29 @@ int main(int argc, char *argv[]) {
stiffnessMatrix.umv(u3, b3);
// {{{ Assemble terms for the nonlinearity
- std::vector<double> normalStress;
- normalStress.resize(grid.size(grid.maxLevel(), dim));
- std::fill(normalStress.begin(), normalStress.end(),
- parset.get<double>("boundary.normalstress"));
+ std::vector<double> a;
+ a.resize(grid.size(grid.maxLevel(), dim));
+ std::fill(a.begin(), a.end(), parset.get<double>("boundary.a"));
std::vector<double> coefficientOfFriction;
coefficientOfFriction.resize(grid.size(grid.maxLevel(), dim));
std::fill(coefficientOfFriction.begin(), coefficientOfFriction.end(),
parset.get<double>("boundary.mu"));
+ std::vector<double> eta;
+ eta.resize(grid.size(grid.maxLevel(), dim));
+ std::fill(eta.begin(), eta.end(), parset.get<double>("boundary.eta"));
+
+ std::vector<double> normalStress;
+ normalStress.resize(grid.size(grid.maxLevel(), dim));
+ std::fill(normalStress.begin(), normalStress.end(),
+ parset.get<double>("boundary.normalstress"));
+
Dune::GlobalLaursenNonlinearity<dim, Dune::LinearFunction>
myGlobalNonlinearity(coefficientOfFriction, normalStress, nodalIntegrals);
+ // Dune::GlobalRuinaNonlinearity<dim>
+ // myGlobalNonlinearity(nodalIntegrals, a, coefficientOfFriction, eta,
+ // normalStress);
// }}}
{
diff --git a/src/one-body-sample.parset b/src/one-body-sample.parset
index ef20f2b1..cb061d7c 100644
--- a/src/one-body-sample.parset
+++ b/src/one-body-sample.parset
@@ -12,5 +12,7 @@ maxiterations = 100000
tolerance = 1e-6
[boundary]
+a = 0.0015
normalstress = 0.1
mu = 0.75
+eta = 1
\ No newline at end of file
--
GitLab