OperatorList.hpp

#pragma once

#include <list>
#include <memory>

#include <amdis/BoundarySubset.hpp>
#include <amdis/ContextGeometry.hpp>
#include <amdis/Operator.hpp>
#include <amdis/typetree/TreeContainer.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 I>
    struct boundary_operator
        : public BoundarySubset<I>
    {
      using type = I;

      boundary_operator(BoundarySubset<I> const& bs)
        : BoundarySubset<I>{bs}
      {}

      using BoundarySubset<I>::BoundarySubset;
    };
  }

  template <class GridView, class Container>
  class OperatorLists
  {
    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>;

    template <class Op>
    struct DataElement
    {
      Op op;
      BoundarySubset<Intersection> bs;

      friend auto localOperator(DataElement const& de)
      {
        return DataElement<decltype(localOperator(de.op))>{localOperator(de.op), de.bs};
      }
    };

    template <class... Nodes>
    struct LocalData
    {
      bool empty() const
      {
        return element_.empty() && boundary_.empty() && intersection_.empty();
      }

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

      void unbind()
      {
        for (auto& op : element_) op.unbind();
        for (auto& op : intersection_) op.unbind();
        for (auto& data : boundary_) data.op.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& op : element_)
          op.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& data : boundary_) {
              if (data.bs(is))
                data.op.assemble(intersectionContext, nodes..., matVec);
            }
          } else {
            // assemble intersection operators
            for (auto& op : intersection_)
              op.assemble(intersectionContext, nodes..., matVec);
          }
        }
      }

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

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

      std::vector<ElementOperator> element_;

      std::vector<IntersectionOperator> intersection_;

      std::vector<DataElement<IntersectionOperator>> boundary_;

      bool assembled_ = false;

      bool changing_ = true;
    };


    template <class... Nodes>
    struct Data
    {
      void update(GridView const& gridView)
      {
        for (auto& e : element_)
          e.update(gridView);
        for (auto& is : intersection_)
          is.update(gridView);
        for (auto& b : boundary_)
          b.op.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(tag::boundary_operator<Intersection> b, Op&& op)
      {
        boundary_.push_back({FWD(op), b});
      }

      friend LocalData<Nodes...> localOperator(Data const& d)
      {
        LocalData<Nodes...> data;
        for (auto const& e : d.element_)
          data.element_.push_back(localOperator(e));
        for (auto const& is : d.intersection_)
          data.intersection_.push_back(localOperator(is));
        for (auto const& b : d.boundary_)
          data.boundary_.push_back(localOperator(b));
        return data;
      }

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

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

      std::vector<ElementOperator> element_;

      std::vector<IntersectionOperator> intersection_;

      std::vector<DataElement<IntersectionOperator>> boundary_;
    };

  public:

    template <class RowNode, class ColNode>
    using MatData = Data<RowNode, ColNode>;

    template <class RowNode, class ColNode>
    using LocalMatData = LocalData<RowNode, ColNode>;

    template <class Node>
    using VecData = Data<Node>;

    template <class Node>
    using LocalVecData = LocalData<Node>;
  };


  template <class RowBasis, class ColBasis, class ElementMatrix>
  using MatrixOperators
    = TypeTree::TreeMatrix<
        OperatorLists<typename RowBasis::GridView,ElementMatrix>::template MatData,
        typename RowBasis::LocalView::TreeCache,
        typename ColBasis::LocalView::TreeCache>;

  template <class Basis, class ElementVector>
  using VectorOperators
    = TypeTree::TreeContainer<
        OperatorLists<typename Basis::GridView,ElementVector>::template VecData,
        typename Basis::LocalView::TreeCache>;

} // end namespace AMDiS