Assembler.hpp

#pragma once

// std c++ headers
#include <memory>

// AMDiS includes
#include <dune/amdis/QuadratureRules.hpp>
#include <dune/amdis/common/TypeDefs.hpp>

namespace AMDiS
{
  // the LocalContext is either an Codim=0-EntityType or an IntersectionType
  template <class GridView, class LocalContext>
  class Assembler
  {
    static constexpr int dim = GridView::dimension;
    using ctype = typename GridView::ctype;

    using LocalQuadratureRules = Dune::QuadratureRules<ctype, LocalContext::mydimension>;
    using QuadratureRule = QuadratureRuleFactory_t<LocalContext, ctype, dim>;
    using Geometry = typename Impl::Get<LocalContext>::Geometry;

  public:
    /// Constructor, initializes the geometry of the element and a
    /// quadrature-rule wrapper
    template <class Operator>
    Assembler(Operator& op, LocalContext const& element, int degree, FirstOrderType type = GRD_PHI)
      : geometry(get_geometry(element))
    {
      int order = op.getQuadratureDegree(geometry.type(), geometry, degree, type);
      auto const& quad = LocalQuadratureRules::rule(geometry.type(), order);

      quadrature.reset(new QuadratureRule(element, quad));
    }

    /// \brief Returns reference to the transformed quadrature rule
    /**
     * For intersections this corresponds to the points in the intersection
     * geometry, but coordinates extended to the whole inside-entity. For
     * elements this is a wrapper around the classical element quadrature rule.
     * Access the underlying dune-quadrature rule, with `quadrature->getRule()`.
     **/
    QuadratureRule const& getQuadrature() const
    {
      return *quadrature;
    }

    /// Return the geometry of the Object
    /**
     * The geometry is either the element geometry, or the geometry of the
     * inside-element for intersections.
     **/
    Geometry const& getGeometry() const
    {
      return geometry;
    }

  protected:

    // transform coords from intersection coords to element coords
    template <class Coordinate>
    decltype(auto) map(Coordinate const& p) const
    {
      return get_position<LocalContext>(geometry, p);
    }

  private:
    /// transformed quadrature rule
    std::shared_ptr<QuadratureRule const> quadrature;

    /// geometry() for entities or geometryInInside() for intersections
    Geometry geometry;

    // TODO: create geometry just once for each element, e.g. in ProblemStat when traversing the grid
  };

} // end namespace AMDiS