#ifndef DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH
#define DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH

#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>

#include <dune/geometry/quadraturerules.hh>

#include <dune/elasticity/assemblers/localenergy.hh>
#include <dune/elasticity/materials/localdensity.hh>

namespace Dune::Elasticity {

template<class GridView, class LocalFiniteElement, class field_type=double>
class LocalIntegralEnergy
  : public Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,GridView::dimension>>>
{
  // grid types
  using DT = typename GridView::Grid::ctype;
  using Entity = typename GridView::template Codim<0>::Entity;

  // some other sizes
  enum {gridDim=GridView::dimension};
  enum {dim=GridView::dimension};

public:

  /** \brief Constructor with a local energy density
    */
  LocalIntegralEnergy(const std::shared_ptr<LocalDensity<dim,field_type,DT>>& ld)
  : localDensity_(ld)
  {}

  /** \brief Virtual destructor */
  virtual ~LocalIntegralEnergy()
  {}

  /** \brief Assemble the energy for a single element */
  field_type energy (const Entity& e,
                     const LocalFiniteElement& localFiniteElement,
                     const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;

protected:
  const std::shared_ptr<LocalDensity<dim,field_type,DT>> localDensity_ = nullptr;

};



template <class GridView, class LocalFiniteElement, class field_type>
field_type
LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>::
energy(const Entity& element,
       const LocalFiniteElement& localFiniteElement,
       const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const
{
  assert(element.type() == localFiniteElement.type());

  field_type energy = 0;

  // store gradients of shape functions and base functions
  std::vector<FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
  std::vector<FieldVector<DT,gridDim> > gradients(localFiniteElement.size());

  int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
                                               : localFiniteElement.localBasis().order() * gridDim;

  const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);

  for (const auto& qp : quad)
  {
    const DT integrationElement = element.geometry().integrationElement(qp.position());

    const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(qp.position());

    // Global position
    auto x = element.geometry().global(qp.position());

    // Get gradients of shape functions
    localFiniteElement.localBasis().evaluateJacobian(qp.position(), referenceGradients);

    // compute gradients of base functions
    for (size_t i=0; i<gradients.size(); ++i)
      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);

    // Deformation gradient
    FieldMatrix<field_type,gridDim,gridDim> deformationGradient(0);
    for (size_t i=0; i<gradients.size(); i++)
      for (int j=0; j<gridDim; j++)
        deformationGradient[j].axpy(localConfiguration[i][j], gradients[i]);

    // Integrate the energy density
    energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
  }

  return energy;
}

}  // namespace Dune::Elasticity


#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH