#ifndef DUNE_TECTONIC_MYDIRECTIONALCONVEXFUNCTION_HH
#define DUNE_TECTONIC_MYDIRECTIONALCONVEXFUNCTION_HH

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

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

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, Dune::Solvers::Interval<double> &D) const {
    Vector uxv = u;
    Arithmetic::addProduct(uxv, x, v);
    phi.directionalSubDiff(uxv, v, D);
    auto const Axmb = A * x - b;
    D[0] += Axmb;
    D[1] += Axmb;
  }

  void domain(Dune::Solvers::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;

  Dune::Solvers::Interval<double> dom;
};

/*
  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, class Nonlinearity>
MyDirectionalConvexFunction<Nonlinearity> restrict(Matrix const &A,
                                                   Vector const &b,
                                                   Vector const &u,
                                                   Vector const &v,
                                                   Nonlinearity const &phi) {
  return MyDirectionalConvexFunction<Nonlinearity>(
      Arithmetic::Axy(A, v, v), Arithmetic::bmAxy(A, b, u, v), phi, u, v);
}

#endif