diff --git a/dune/tectonic/globalruinanonlinearity.hh b/dune/tectonic/globalruinanonlinearity.hh new file mode 100644 index 0000000000000000000000000000000000000000..aed666530af1f2c876e53ad5d7c96991349aa652 --- /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 a3614c03fe242b97015a62813ac0692aeaf0f600..cc582278b260a6f88357cd1ebe47397bee066eb9 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 3b0b6146ff8f763aed1cdfe906e28d28977d92e7..720f80be405ce9e3ecb2da4bb8e45f3b1d9a5156 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 ef20f2b1110bcb955fab37a8347b1b7d78eff89f..cb061d7c87ee268ed4c52572ecacce294dda06bc 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