// Copied from dune/tnnmg/problem-classes/directionalconvexfunction.hh
// Allows phi to be const

#ifndef MY_DIRECTIONAL_CONVEX_FUNCTION_HH
#define MY_DIRECTIONAL_CONVEX_FUNCTION_HH

#include <dune/fufem/arithmetic.hh>
#include <dune/fufem/interval.hh>

/*
  1/2 <A(u + hv),u + hv> - <b, u + hv>
  = 1/2 <Av,v> h^2 - <b - Au, v> h + const.

  localA = <Av,v>
  localb = <b - Au, v>
*/

template <class Matrix, class Vector>
double computeDirectionalA(Matrix const &A, Vector const &v) {
  return Arithmetic::Axy(A, v, v);
}

template <class Matrix, class Vector>
double computeDirectionalb(Matrix const &A, Vector const &b, Vector const &u,
                           Vector const &v) {
  return Arithmetic::bmAxy(A, b, u, v);
}

template <class Nonlinearity> class MyDirectionalConvexFunction {
public:
  using Vector = typename Nonlinearity::VectorType;
  using Matrix = typename Nonlinearity::MatrixType;

  MyDirectionalConvexFunction(double A, double b, Nonlinearity const &phi,
                              Vector const &u, Vector const &v)
      : A(A), b(b), phi(phi), u(u), v(v) {
    phi.directionalDomain(u, v, dom);
  }

  double quadraticPart() const { return A; }

  double linearPart() const { return b; }

  void subDiff(double x, Interval<double> &D) const {
    Vector uxv = u;
    Arithmetic::addProduct(uxv, x, v);
    phi.directionalSubDiff(uxv, v, D);
    D[0] += A * x - b;
    D[1] += A * x - b;
  }

  void domain(Interval<double> &domain) const {
    domain[0] = this->dom[0];
    domain[1] = this->dom[1];
  }

  double A;
  double b;

private:
  Nonlinearity const &phi;
  Vector const &u;
  Vector const &v;

  Interval<double> dom;
};

#endif