BiLinearForm.hpp

#pragma once

#include <cmath>

#include <amdis/LinearAlgebra.hpp>
#include <amdis/Observer.hpp>
#include <amdis/OperatorList.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/FlatMatrix.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
#include <amdis/typetree/TreePath.hpp>

namespace AMDiS
{
  template <class RB, class CB, class T = double, class Traits = BackendTraits<RB>>
  class BiLinearForm
      : public MatrixFacade<T, Traits::template MatrixImpl>
      , private ObserverSequence<event::adapt,2>
  {
    using Self = BiLinearForm;
    using Super = MatrixFacade<T, Traits::template MatrixImpl>;

  public:
    using RowBasis = RB;

    using ColBasis = CB;

    using CoefficientType = T;

    using ElementMatrix = FlatMatrix<CoefficientType>;

    using SparsityPattern = typename Traits::SparsityPattern;

  public:
    BiLinearForm(std::shared_ptr<RB> const& rowBasis, std::shared_ptr<CB> const& colBasis)
      : Super(*rowBasis, *colBasis)
      , ObserverSequence<event::adapt,2>(*rowBasis, *colBasis)
      , symmetry_(SymmetryStructure::unknown)
      , rowBasis_(rowBasis)
      , colBasis_(colBasis)
      , updatePattern_(true)
    {
      auto const rowSize = rowBasis_->localView().maxSize();
      auto const colSize = colBasis_->localView().maxSize();
      elementMatrix_.resize(rowSize, colSize);
    }

    template <class RB_, class CB_,
      REQUIRES(Concepts::Similar<RB_,RB>),
      REQUIRES(Concepts::Similar<CB_,CB>)>
    BiLinearForm(RB_&& rowBasis, CB_&& colBasis)
      : BiLinearForm(Dune::wrap_or_move(FWD(rowBasis)), Dune::wrap_or_move(FWD(colBasis)))
    {}

    template <class RB_ = RB, class CB_ = CB,
      REQUIRES(std::is_same_v<RB_,CB_>)>
    explicit BiLinearForm(std::shared_ptr<RB> const& rowBasis)
      : BiLinearForm(rowBasis, rowBasis)
    {}

    template <class RB_,
      REQUIRES(Concepts::Similar<RB_,RB>)>
    explicit BiLinearForm(RB_&& rowBasis)
      : BiLinearForm(Dune::wrap_or_move(FWD(rowBasis)))
    {}

    RowBasis const& rowBasis() const { return *rowBasis_; }
    ColBasis const& colBasis() const { return *colBasis_; }


    template <class ContextTag, class Operator, class Row, class Col>
    void addOperator(ContextTag contextTag, Operator const& op, Row row, Col col);

    // Add an operator to be assembled on the elements of the grid
    template <class Operator, class Row = RootTreePath, class Col = RootTreePath>
    void addOperator(Operator const& op, Row row = {}, Col col = {})
    {
      using E = typename RowBasis::LocalView::Element;
      addOperator(tag::element_operator<E>{}, op, row, col);
    }

    // Add an operator to be assembled on the intersections of the grid
    template <class Operator, class Row = RootTreePath, class Col = RootTreePath>
    void addIntersectionOperator(Operator const& op, Row row = {}, Col col = {})
    {
      using I = typename RowBasis::LocalView::GridView::Intersection;
      addOperator(tag::intersection_operator<I>{}, op, row, col);
    }


    void setSymmetryStructure(SymmetryStructure symm)
    {
      updatePattern_ = updatePattern_ || (symmetry_ != symm);
      symmetry_ = symm;
    }

    void setSymmetryStructure(std::string str)
    {
      setSymmetryStructure(symmetryStructure(str));
    }

    template <class RowLocalView, class ColLocalView, class LocalOperators,
      REQUIRES(Concepts::LocalView<RowLocalView>),
      REQUIRES(Concepts::LocalView<ColLocalView>)>
    void assemble(RowLocalView const& rowLocalView,
                  ColLocalView const& colLocalView,
                  LocalOperators& localOperators);

    void assemble();


    void init(bool forcePatternRebuild = false)
    {
      if (updatePattern_ || forcePatternRebuild)
      {
        Super::init(SparsityPattern{*rowBasis_, *colBasis_, symmetry_});
        updatePattern_ = false;
      }
      else
      {
        Super::init();
      }
    }

    void updateImpl(event::adapt e, index_t<0> i) override { updateImpl2(*rowBasis_); }
    void updateImpl(event::adapt e, index_t<1> i) override { updateImpl2(*colBasis_); }

    auto& operators() { return operators_; }

  private:
    template <class Basis>
    void updateImpl2(Basis const& basis)
    {
      if (!updatePattern_)
        Recursive::forEach(operators_, [&](auto& op) { op.update(basis.gridView()); });
      updatePattern_ = true;
    }

  protected:
    SymmetryStructure symmetry_;

    ElementMatrix elementMatrix_;

    MatrixOperators<RowBasis,ColBasis,ElementMatrix> operators_;

    std::shared_ptr<RowBasis const> rowBasis_;
    std::shared_ptr<ColBasis const> colBasis_;

  private:
    // The pattern is rebuilt on calling init if this flag is set
    bool updatePattern_;
  };


  // deduction guides
  template <class RB, class CB>
  BiLinearForm(RB&&, CB&&)
    -> BiLinearForm<Underlying_t<RB>, Underlying_t<CB>>;

  template <class RB>
  BiLinearForm(RB&&)
    -> BiLinearForm<Underlying_t<RB>, Underlying_t<RB>>;


  template <class T = double, class RB, class CB>
  auto makeBiLinearForm(RB&& rowBasis, CB&& colBasis)
  {
    using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<CB>, T>;
    return BLF{FWD(rowBasis), FWD(colBasis)};
  }

  template <class T = double, class RB>
  auto makeBiLinearForm(RB&& rowBasis)
  {
    using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<RB>, T>;
    return BLF{FWD(rowBasis)};
  }

} // end namespace AMDiS

#include <amdis/BiLinearForm.inc.hpp>