Forked from
agnumpde / dune-tectonic
966 commits behind the upstream repository.
-
Elias Pipping authored
We now re-assemble the global nonlinearity for every step
Elias Pipping authoredWe now re-assemble the global nonlinearity for every step
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
nicefunction.hh 5.58 KiB
#ifndef NICE_FUNCTION_HH
#define NICE_FUNCTION_HH
#include <algorithm>
#include <cmath>
#include <limits>
#include <dune/common/exceptions.hh>
#include <dune/common/function.hh>
namespace Dune {
class NiceFunction : public VirtualFunction<double, double> {
public:
virtual ~NiceFunction() {}
double virtual leftDifferential(double s) const = 0;
double virtual rightDifferential(double s) const = 0;
double virtual second_deriv(double x) const {
DUNE_THROW(NotImplemented, "second derivative not implemented");
}
double virtual regularity(double s) const {
DUNE_THROW(NotImplemented, "regularity not implemented");
}
// Whether H(|.|) is smooth at zero
bool virtual smoothesNorm() const { return false; }
};
class RuinaFunction : public NiceFunction {
public:
RuinaFunction(double coefficient, double a, double mu, double eta,
double normalStress, double state)
: coefficient(coefficient),
a(a),
mu(mu),
rho(exp(-mu / a)),
eta(eta),
normalStress(normalStress),
state(state) {}
/*
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 + state);
y = 1 / eta * normalStress * h + a * rho + mu + state;
y *= coefficient;
}
/*
(leaving some terms aside): with s > rho
1/eta d/dx [ a * (s log s - s) + mu s ] where s = eta x
= 1/eta [ a * (log (eta x) * eta) + eta mu ]
= a * log(eta x) + mu
*/
double virtual leftDifferential(double s) const {
if (eta * s <= rho)
return 0;
return coefficient * normalStress * (a * std::log(eta * s) + mu + state);
}
/* see above */
double virtual rightDifferential(double s) const {
if (eta * s <= rho)
return 0;
return coefficient * normalStress * (a * std::log(eta * s) + mu + state);
}
/*
d/dx a * log(eta x) + mu
= a * 1/(eta x) * eta
= a/x
*/
double virtual second_deriv(double s) const {
// includes the case eta * s = rho for which there is no second derivative
if (eta * s <= rho)
return 0;
else
return coefficient * normalStress * (a / s);
}
double virtual regularity(double s) const {
if (eta * s == rho)
return std::numeric_limits<double>::infinity();
return std::abs(second_deriv(s));
}
private:
double const coefficient;
double const a;
double const mu;
double const eta;
double const normalStress;
double const state;
double const rho;
};
class LinearFunction : public NiceFunction {
public:
LinearFunction(double a) : coefficient(a) {}
void virtual evaluate(double const &x, double &y) const {
y = coefficient * x;
}
double virtual leftDifferential(double s) const { return coefficient; }
double virtual rightDifferential(double s) const { return coefficient; }
double virtual second_deriv(double) const { return 0; }
double virtual regularity(double s) const { return 0; }
private:
double const coefficient;
};
template <int slope> class SampleFunction : public NiceFunction {
public:
void virtual evaluate(double const &x, double &y) const {
y = (x < 1) ? x : (slope * (x - 1) + 1);
}
double virtual leftDifferential(double s) const {
return (s <= 1) ? 1 : slope;
}
double virtual rightDifferential(double s) const {
return (s < 1) ? 1 : slope;
}
};
class SteepFunction : public NiceFunction {
public:
void virtual evaluate(double const &x, double &y) const { y = 100 * x; }
double virtual leftDifferential(double s) const { return 100; }
double virtual rightDifferential(double s) const { return 100; }
};
class TrivialFunction : public NiceFunction {
public:
void virtual evaluate(double const &x, double &y) const { y = 0; }
double virtual leftDifferential(double) const { return 0; }
double virtual rightDifferential(double) const { return 0; }
double virtual second_deriv(double) const { return 0; }
double virtual regularity(double) const { return 0; }
bool virtual smoothesNorm() const { return true; }
};
// slope in [n-1,n] is n
class HorribleFunction : public NiceFunction {
public:
void virtual evaluate(double const &x, double &y) const {
double const fl = floor(x);
double const sum = fl * (fl + 1) / 2;
y = sum + (fl + 1) * (x - fl);
}
double virtual leftDifferential(double x) const {
double const fl = floor(x);
if (x - fl < 1e-14)
return fl;
else
return fl + 1;
}
double virtual rightDifferential(double x) const {
double const c = ceil(x);
if (c - x < 1e-14)
return c + 1;
else
return c;
}
};
// slope in [n-1,n] is log(n+1)
class HorribleFunctionLogarithmic : public NiceFunction {
public:
void virtual evaluate(double const &x, double &y) const {
y = 0;
size_t const fl = floor(x);
for (size_t i = 1; i <= fl;)
y += log(++i); // factorials grow to fast so we compute this incrementally
y += log(fl + 2) * (x - fl);
}
double virtual leftDifferential(double x) const {
double const fl = floor(x);
if (x - fl < 1e-14)
return log(fl + 1);
else
return log(fl + 2);
}
double virtual rightDifferential(double x) const {
double const c = ceil(x);
if (c - x < 1e-14)
return log(c + 2);
else
return log(c + 1);
}
};
}
#endif