Assembler.hpp

#pragma once

#include <list>
#include <memory>

#include <amdis/ContextGeometry.hpp>
#include <amdis/BoundaryOperator.hpp>
#include <amdis/Operator.hpp>
#include <amdis/algorithm/Map.hpp>

namespace AMDiS
{
  namespace tag
  {
    template <class E> struct element_operator { using type = E; };
    template <class I> struct intersection_operator { using type = I; };
  }

  template <class GridView, class Container, class... Nodes>
  class LocalAssembler
  {
    using Element = typename GridView::template Codim<0>::Entity;
    using Intersection = typename GridView::Intersection;

    using ElementTraits = Impl::OperatorTraits<GridView,Element,Container>;
    using IntersectionTraits = Impl::OperatorTraits<GridView,Intersection,Container>;

    using ElementLocalOperator = LocalOperator<ElementTraits, Nodes...>;

    using IntersectionLocalOperator = LocalOperator<IntersectionTraits, Nodes...>;

  public:
    LocalAssembler(std::vector<ElementLocalOperator> element,
                   std::vector<IntersectionLocalOperator> intersection,
                   std::vector<BoundaryLocalOperator<IntersectionLocalOperator,Intersection>> boundary)
      : element_(std::move(element))
      , intersection_(std::move(intersection))
      , boundary_(std::move(boundary))
    {}

    bool empty() const
    {
      return element_.empty() && boundary_.empty() && intersection_.empty();
    }

    void bind(Element const& elem)
    {
      for (auto& lop : element_) lop.bind(elem);
      for (auto& lop : intersection_) lop.bind(elem);
      for (auto& lop : boundary_) lop.bind(elem);
    }

    void unbind()
    {
      for (auto& lop : element_) lop.unbind();
      for (auto& lop : intersection_) lop.unbind();
      for (auto& lop : boundary_) lop.unbind();
      assembled_ = true;
    }

    void assemble(GlobalContext<GridView> const& context, Nodes const&... nodes,
                  Container& matVec) const
    {
      // do not assemble in case nothing to do
      if ((assembled_ && !changing_) || empty())
        return;

      // create a context for the element
      ContextGeometry elementContext{context.element(), context.element(), context.geometry()};

      // assemble element operators
      for (auto const& lop : element_)
        lop.assemble(elementContext, nodes..., matVec);

      if (intersection_.empty() && (boundary_.empty() || !context.element().hasBoundaryIntersections()))
        return;

      // assemble intersection operators
      for (auto const& is : intersections(context.gridView(), context.element()))
      {
        // create a context for the intersection
        ContextGeometry intersectionContext{is, context.element(), context.geometry()};

        if (is.boundary()) {
          // assemble boundary operators
          for (auto& lop : boundary_) {
            lop.assemble(intersectionContext, nodes..., matVec);
          }
        } else {
          // assemble intersection operators
          for (auto& lop : intersection_)
            lop.assemble(intersectionContext, nodes..., matVec);
        }
      }
    }

  private:
    std::vector<ElementLocalOperator> element_;

    std::vector<IntersectionLocalOperator> intersection_;

    std::vector<BoundaryLocalOperator<IntersectionLocalOperator,Intersection>> boundary_;

    bool assembled_ = false;

    bool changing_ = true;
  };


  template <class GridView, class Container, class... Nodes>
  class Assembler
  {
    using Element = typename GridView::template Codim<0>::Entity;
    using Intersection = typename GridView::Intersection;

    using ElementTraits = Impl::OperatorTraits<GridView,Element,Container>;
    using IntersectionTraits = Impl::OperatorTraits<GridView,Intersection,Container>;

    using ElementOperator = Operator<ElementTraits, Nodes...>;

    using IntersectionOperator = Operator<IntersectionTraits, Nodes...>;

  public:
    Assembler () = default;

    void update(GridView const& gridView)
    {
      for (auto& e : element_)
        e.update(gridView);
      for (auto& is : intersection_)
        is.update(gridView);
      for (auto& b : boundary_)
        b.update(gridView);
    }

    template <class Op>
    void push(tag::element_operator<Element>, Op&& op)
    {
      element_.emplace_back(FWD(op));
    }

    template <class Op>
    void push(tag::intersection_operator<Intersection>, Op&& op)
    {
      intersection_.emplace_back(FWD(op));
    }

    template <class Op>
    void push(BoundarySubset<Intersection> b, Op&& op)
    {
      boundary_.push_back({FWD(op), std::move(b)});
    }

    friend LocalAssembler<GridView,Container,Nodes...> localAssembler(Assembler const& a)
    {
      auto toLocalOperator = [](auto const& op) { return localOperator(op); };
      return {Recursive::map(toLocalOperator, a.element_),
              Recursive::map(toLocalOperator, a.intersection_),
              Recursive::map(toLocalOperator, a.boundary_)};
    }

  private:
    std::vector<ElementOperator> element_{};

    std::vector<IntersectionOperator> intersection_{};

    std::vector<BoundaryOperator<IntersectionOperator,Intersection>> boundary_{};
  };


  template <class RowBasis, class ColBasis, class ElementMatrix>
  using MatrixAssembler
    = Assembler<typename RowBasis::GridView,ElementMatrix,
        typename RowBasis::LocalView::TreeCache,
        typename ColBasis::LocalView::TreeCache>;

  template <class Basis, class ElementVector>
  using VectorAssembler
    = Assembler<typename Basis::GridView,ElementVector,
        typename Basis::LocalView::TreeCache>;

} // end namespace AMDiS