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