Assembler.hpp

#pragma once

#include <memory>
#include <tuple>

#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>

#include <amdis/DirichletBC.hpp>
#include <amdis/LinearAlgebra.hpp>
#include <amdis/LocalAssemblerBase.hpp>
#include <amdis/common/Mpl.hpp>
#include <amdis/common/TypeDefs.hpp>

namespace AMDiS
{
  template <class Traits>
  class Assembler
  {
    using GlobalBasis = typename Traits::GlobalBasis;

    /// The grid view the global FE basis lives on
    using GridView = typename GlobalBasis::GridView;

  public:

    /// Constructor, stores a shared-pointer to the feSpaces
    Assembler(GlobalBasis& globalBasis,
              MatrixOperators<GlobalBasis>& matrixOperators,
              VectorOperators<GlobalBasis>& rhsOperators,
              Constraints<GlobalBasis>& constraints)
      : globalBasis_(globalBasis)
      , matrixOperators_(matrixOperators)
      , rhsOperators_(rhsOperators)
      , constraints_(constraints)
    {}

    /// Assemble the linear system
    template <class SystemMatrixType, class SystemVectorType>
    void assemble(
        SystemMatrixType& matrix,
        SystemVectorType& solution,
        SystemVectorType& rhs,
        bool asmMatrix, bool asmVector);


  private:

    /// Sets the system to zero and initializes all operators and boundary conditions
    template <class SystemMatrixType, class SystemVectorType>
    void initMatrixVector(
        SystemMatrixType& matrix,
        SystemVectorType& solution,
        SystemVectorType& rhs,
        bool asmMatrix, bool asmVector) const;

    /// Assemble operators on an element, by passing the element/intersection to
    /// `elementAssembler` functor.
    template <class Element, class Operators, class ElementAssembler>
    void assembleElementOperators(
        Element const& element,
        Operators& operators,
        ElementAssembler const& elementAssembler) const;

    /// Finish insertion into the matrix and assembles boundary conditions
    /// Return the number of nonzeros assembled into the matrix
    template <class SystemMatrixType, class SystemVectorType>
    std::size_t finishMatrixVector(
        SystemMatrixType& matrix,
        SystemVectorType& solution,
        SystemVectorType& rhs,
        bool asmMatrix, bool asmVector) const;


    /// Return the element the LocalViews are bound to
    template <class LocalView0, class... LovalViews>
    auto const& getElement(LocalView0 const& localView, LovalViews const&...) const
    {
      return localView.element();
    }

    /// Return the gridView the localViews are bound to
    template <class LocalView0, class... LovalViews>
    auto const& getGridView(LocalView0 const& localView, LovalViews const&...) const
    {
      return globalBasis_.gridView();
    }


  private:

    GlobalBasis& globalBasis_;
    MatrixOperators<GlobalBasis>& matrixOperators_;
    VectorOperators<GlobalBasis>& rhsOperators_;
    Constraints<GlobalBasis>& constraints_;

    //TODO: add caching of localBases
  };

} // end namespace AMDiS

#include "Assembler.inc.hpp"