Assembler.hpp

#pragma once

#include <functional>
#include <memory>
#include <type_traits>

#include <dune/common/shared_ptr.hh>
#include <dune/geometry/quadraturerules.hh>

#include <amdis/ContextGeometry.hpp>
#include <amdis/AssemblerInterface.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/typetree/FiniteElementType.hpp>

namespace AMDiS
{
  /// Implementation of interface \ref AssemblerBase
  /**
   * \tparam Traits  See \ref DefaultAssemblerTraits
   **/
  template <class Traits, class Operator, class... Nodes>
  class Assembler
      : public AssemblerInterface<Traits, Nodes...>
  {
    using Super = AssemblerInterface<Traits, Nodes...>;

    using LocalContext = typename Traits::LocalContext;
    using Element = typename Super::Element;
    using Geometry = typename Super::Geometry;
    using ElementContainer = typename Traits::ElementContainer;

  public:

    /// Constructor. Stores a copy of operator `op`.
    explicit Assembler(Operator const& op)
      : op_(Dune::wrap_or_move(op))
    {}

    /// Constructor. Stores a copy of operator `op`.
    explicit Assembler(Operator&& op)
      : op_(Dune::wrap_or_move(std::move(op)))
    {}

    /// Constructor. Stores the reference to the operator.
    explicit Assembler(std::reference_wrapper<Operator> op)
      : op_(Dune::wrap_or_move(op.get()))
    {}

    /// \brief Implementation of \ref AssemblerInterface::bind.
    /**
     * Binds the operator `op_` to the `element` and `geometry` and
     * stores point to the `element` and `geometry`.
     **/
    void bind(Element const& element, Geometry const& geometry) final
    {
      element_ = &element;
      geometry_ = &geometry;
      op_->bind(element, geometry);
    }

    /// \brief Implementation of \ref AssemblerBase::unbind
    /**
     * Unbinds the operator `op_` and sets \ref element_ and \ref geometry_
     * to nullptr.
     **/
    void unbind() final
    {
      op_->unbind();
      geometry_ = nullptr;
      element_ = nullptr;
    }

    /// Implementation of \ref AssemblerBase::assemble
    /**
     * Stores geometry and localGeometry and calls
     * \ref calculateElementVector or \ref calculateElementMatrix on the
     * vector or matrix operator, respectively.
     **/
    void assemble(LocalContext const& localContext, Nodes const&... nodes,
                  ElementContainer& ElementContainer) final
    {
      ContextGeometry<LocalContext> contextGeo{localContext, element(), geometry()};
      assembleImpl(contextGeo, nodes..., ElementContainer);
    }


#ifndef DOXYGEN

  protected: // implementation detail

    // matrix assembling
    template <class CG, class RN, class CN, class Mat>
    void assembleImpl(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
    {
      op_->calculateElementMatrix(contextGeo, rowNode, colNode, elementMatrix);
    }

    // vector assembling
    template <class CG, class Node, class Vec>
    void assembleImpl(CG const& contextGeo, Node const& node, Vec& elementVector)
    {
      op_->calculateElementVector(contextGeo, node, elementVector);
    }

#endif // DOXYGEN

  public:

    /// return the bound entity (of codim 0)
    Element const& element() const
    {
      assert( element_ );
      return *element_;
    }

    /// return the geometry of the bound element
    Geometry const& geometry() const
    {
      assert( geometry_ );
      return *geometry_;
    }


  private:

    /// the stored operator, implementing \ref GridFunctionOperatorBase
    std::shared_ptr<Operator> op_;

    Element const* element_ = nullptr;
    Geometry const* geometry_ = nullptr;
  };


  /// Generate a \ref Assembler on a given LocalContext `LC` (element or intersection)
  template <class Traits, class Operator, class... Nodes>
  auto makeAssembler(Operator&& op, Nodes const&...)
  {
    return Assembler<Traits, Underlying_t<Operator>, Nodes...>{FWD(op)};
  }

} // end namespace AMDiS