// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_MULTIMESH_LOCALGEOMETRY_HH
#define DUNE_MULTIMESH_LOCALGEOMETRY_HH

#include <dune/common/fmatrix.hh>
#include <dune/common/typetraits.hh>
#include <dune/grid/common/geometry.hh>

namespace Dune
{
  template <int mydim, int coorddim, class HostGrid>
  class MultiMeshLocalGeometry
      : public GeometryDefaultImplementation<mydim, coorddim, HostGrid, MultiMeshLocalGeometry>
  {
  private:
    using ctype = typename HostGrid::ctype;

  public:
    // The codimension of this entitypointer wrt the host grid
    enum { codimension = HostGrid::dimension - mydim };
    enum { dimensionworld = HostGrid::dimensionworld };

    // select appropriate hostgrid geometry via typeswitch
    using HostLocalGeometry = typename HostGrid::Traits::template Codim<codimension>::LocalGeometry;

    /// entity in the host grid
    using HostGridEntity = typename HostGrid::Traits::template Codim<codimension>::Entity;

    //! type of jacobian transposed
    using JacobianInverseTransposed = typename HostLocalGeometry::JacobianInverseTransposed;
    using JacobianTransposed = typename HostLocalGeometry::JacobianTransposed;


    /// Constructor from host entities
    MultiMeshLocalGeometry (const HostGridEntity& sourceEntity, const HostGridEntity& targetEntity)
    {
      hostLocalGeometries_.reserve(sourceEntity.level() - targetEntity.level());
      auto entity = sourceEntity;
      for (int l = sourceEntity.level(); l > targetEntity.level(); --l) {
        assert( entity.hasFather() );
        hostLocalGeometries_.push_back(entity.geometryInFather());
        if (l > targetEntity.level()+1)
          entity = entity.father();
      }
    }


    /// Return the element type identifier
    GeometryType type () const
    {
      return hostLocalGeometries_[0].type();
    }

    /// Return wether we have an affine mapping
    bool affine () const
    {
      return std::all_of(hostLocalGeometries_.begin(), hostLocalGeometries_.end(),
        [](auto const& localGeo) { return localGeo.affine(); });
    }

    /// Returns true if the point is in the current element
    bool checkInside (const FieldVector<ctype, mydim>& local) const
    {
      return hostLocalGeometries_[0].checkInside(local);
    }

    /// Return the number of corners of this element. Corners are numbered 0...n-1
    int corners () const
    {
      return hostLocalGeometries_[0].corners();
    }

    /// Access to coordinates of corners. Index is the number of the corner
    const FieldVector<ctype, coorddim> corner (int i) const
    {
      FieldVector<ctype, coorddim> c = hostLocalGeometries_[0].corner(i);
      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j)
        c = hostLocalGeometries_[j].global(c);
      return c;
    }


    /// Maps a local coordinate within reference element to global coordinate in element
    FieldVector<ctype, coorddim> global (const FieldVector<ctype, mydim>& local) const
    {
      auto c = hostLocalGeometries_[0].global(local);
      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j)
        c = hostLocalGeometries_[j].global(c);
      return c;
    }

    /// Maps a global coordinate within the element to a local coordinate in its reference element
    FieldVector<ctype, mydim> local (const FieldVector<ctype, coorddim>& global) const
    {
      auto c = hostLocalGeometries_.back().local(global);
      std::size_t n = hostLocalGeometries_.size()-1;
      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j)
        c = hostLocalGeometries_[n-j].local(c);
      return c;
    }

    /// Return the factor appearing in the integral transformation formula
    ctype integrationElement (const FieldVector<ctype, mydim>& local) const
    {
      ctype dx = 1;
      auto c = local;
      for (std::size_t j = 0; j < hostLocalGeometries_.size(); ++j) {
        dx *= hostLocalGeometries_[j].integrationElement(c);
        if (j < hostLocalGeometries_.size()-1)
          c = hostLocalGeometries_[j].global(local);
      }
      return dx;
    }

    /// The Jacobian matrix of the mapping from the reference element to this element
    JacobianInverseTransposed jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const
    {
      auto J = hostLocalGeometries_[0].jacobianInverseTransposed(local);
      auto c = hostLocalGeometries_[0].global(local);
      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j) {
        J.leftmultiply(hostLocalGeometries_[j].jacobianInverseTransposed(c));
        if (j < hostLocalGeometries_.size()-1)
          c = hostLocalGeometries_[j].global(local);
      }
      return J;
    }

    /// Return the transposed of the Jacobian
    JacobianTransposed jacobianTransposed (const FieldVector<ctype, mydim>& local) const
    {
      auto J = hostLocalGeometries_[0].jacobianTransposed(local);
      auto c = hostLocalGeometries_[0].global(local);
      for (std::size_t j = 1; j < hostLocalGeometries_.size(); ++j) {
        J.rightmultiply(hostLocalGeometries_[j].jacobianTransposed(c));
        if (j < hostLocalGeometries_.size()-1)
          c = hostLocalGeometries_[j].global(local);
      }
      return J;
    }

  private:
    std::vector<HostLocalGeometry> hostLocalGeometries_;
  };

}  // end namespace Dune

#endif // DUNE_MULTIMESH_LOCALGEOMETRY_HH