#include "ageinglawstateupdater.hh"
#include "../tobool.hh"

template <class ScalarVector, class Vector>
AgeingLawStateUpdater<ScalarVector, Vector>::AgeingLawStateUpdater(
    ScalarVector const &_alpha_initial, Dune::BitSetVector<1> const &_nodes,
    double _L, double _V0)
    : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}

template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::nextTimeStep() {
  alpha_o = alpha;
}

template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::setup(double _tau) {
  tau = _tau;
}

/*
  Compute [ 1-\exp(c*x) ] / x under the assumption that x is
  non-negative
*/
double liftSingularity(double c, double x) {
  if (x <= 0)
    return -c;
  else
    return -std::expm1(c * x) / x;
}

template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::solve(
    Vector const &velocity_field) {
  for (size_t i = 0; i < nodes.size(); ++i) {
    if (not toBool(nodes[i]))
      continue;

    double const V = velocity_field[i].two_norm();
    double const mtoL = -tau / L;
    alpha[i] = std::log(std::exp(alpha_o[i] + V * mtoL) +
                        V0 * liftSingularity(mtoL, V));
  }
}

template <class ScalarVector, class Vector>
void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha(
    ScalarVector &_alpha) {
  _alpha = alpha;
}

template <class ScalarVector, class Vector>
std::shared_ptr<StateUpdater<ScalarVector, Vector>>
AgeingLawStateUpdater<ScalarVector, Vector>::clone() const {
  return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(*this);
}