// 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 MatrixType, class VectorType>
double computeDirectionalA(MatrixType const &A, VectorType const &v) {
  return Arithmetic::Axy(A, v, v);
}

template <class MatrixType, class VectorType>
double computeDirectionalb(MatrixType const &A, VectorType const &b,
                           VectorType const &u, VectorType const &v) {
  VectorType tmp = b;
  Arithmetic::addProduct(tmp, -1.0, A, u);
  return tmp * v;
}

template <class NonlinearityType> class MyDirectionalConvexFunction {
public:
  using VectorType = typename NonlinearityType::VectorType;
  using MatrixType = typename NonlinearityType::MatrixType;

  MyDirectionalConvexFunction(double A, double b, NonlinearityType const &phi,
                              VectorType const &u, VectorType 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 {
    VectorType 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:
  NonlinearityType const &phi;
  VectorType const &u;
  VectorType const &v;

  Interval<double> dom;
};

#endif