-
Elias Pipping authoredElias Pipping authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
nicefunction.hh 5.21 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:
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; }
void virtual evaluate(double const &x, double &y) const = 0;
};
class RuinaFunction : public NiceFunction {
public:
RuinaFunction(double coefficient, double a, double mu, double eta,
double normalStress, double b, double state, double L, double h)
: coefficient(coefficient),
a(a),
mu(mu),
eta(eta),
normalStress(normalStress),
compound_state(b * (state - std::log(eta * L))),
h(h),
rho(exp(-(mu + compound_state) / a)) {}
// TODO: untested
void virtual evaluate(double const &x, double &y) const {
assert(false);
double const arg = eta * x / h;
double const expstar = (arg == 0) ? 0 : arg * std::log(arg) - arg;
double const gamma =
normalStress * (mu + a * std::log(eta) + compound_state);
y = a * normalStress * expstar + gamma * arg;
y *= coefficient * h * normalStress / eta;
}
double virtual leftDifferential(double s) const {
double const arg = eta * s / h;
if (arg == 0)
return 0;
double const gamma =
normalStress * (mu + a * std::log(eta) + compound_state);
return coefficient * (a * normalStress * std::log(arg) + gamma);
}
/* see above */
double virtual rightDifferential(double s) const {
return leftDifferential(s);
}
double virtual second_deriv(double s) const {
return coefficient * a * normalStress / s;
}
double virtual regularity(double s) const {
if (s == 0)
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 compound_state;
double const h;
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 += std::log(
++i); // factorials grow to fast so we compute this incrementally
y += std::log(fl + 2) * (x - fl);
}
double virtual leftDifferential(double x) const {
double const fl = floor(x);
if (x - fl < 1e-14)
return std::log(fl + 1);
else
return std::log(fl + 2);
}
double virtual rightDifferential(double x) const {
double const c = ceil(x);
if (c - x < 1e-14)
return std::log(c + 2);
else
return std::log(c + 1);
}
};
}
#endif