#ifndef ELLIPTIC_ENERGY_HH
#define ELLIPTIC_ENERGY_HH

#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/stdstreams.hh>

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

#include "localfriction.hh"

template <size_t dim> class EllipticEnergy {
public:
  using LocalVector = Dune::FieldVector<double, dim>;
  using LocalMatrix = Dune::FieldMatrix<double, dim, dim>;

  using Nonlinearity = LocalFriction<dim>;

  EllipticEnergy(LocalMatrix const &A, LocalVector const &b,
                 std::shared_ptr<Nonlinearity const> phi,
                 typename Dune::BitSetVector<dim>::const_reference ignore)
      : A(A), b(b), phi(phi), ignore(ignore) {}

  double operator()(LocalVector const &v) const {
    return computeEnergy(A, v, b) + (*phi)(v);
  }

  LocalMatrix const &A;
  LocalVector const &b;
  std::shared_ptr<Nonlinearity const> const phi;
  typename Dune::BitSetVector<dim>::const_reference const ignore;

  void descentDirection(LocalVector const &x, LocalVector &y) const {
    if (x.two_norm() >= 0.0) {
      A.mv(x, y);
      y -= b;
      phi->addGradient(x, y);

      for (size_t i = 0; i < dim; ++i)
        if (ignore[i])
          y[i] = 0;
      y *= -1;
    } else {
      A.mv(x, y);
      y -= b;

      for (size_t i = 0; i < dim; ++i)
        if (ignore[i])
          y[i] = 0;
      y *= -1;

      // The contribution of phi to the directional derivative is the
      // same for all directions here! Choose the one that is best of
      // the smooth part and check if we have overall decline
      Dune::Solvers::Interval<double> D;
      phi->directionalSubDiff(x, y, D);
      if (D[1] - y.two_norm2() >= 0.0)
        y = 0.0;
    }
  }
};
#endif