ProblemStat.hpp

#pragma once

#include <list>
#include <map>
#include <memory>
#include <string>
#include <tuple>
#include <utility>
#include <vector>

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

#include <dune/grid/common/grid.hh>

#include <amdis/AdaptInfo.hpp>
#include <amdis/AdaptiveGrid.hpp>
#include <amdis/BiLinearForm.hpp>
#include <amdis/CreatorInterface.hpp>
#include <amdis/CreatorMap.hpp>
#include <amdis/DirichletBC.hpp>
#include <amdis/DOFVector.hpp>
//#include <amdis/Estimator.hpp>
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
#include <amdis/LinearForm.hpp>
#include <amdis/Marker.hpp>
#include <amdis/MeshCreator.hpp>
#include <amdis/OperatorList.hpp>
#include <amdis/PeriodicBC.hpp>
#include <amdis/ProblemStatBase.hpp>
#include <amdis/ProblemStatTraits.hpp>
#include <amdis/StandardProblemIteration.hpp>
#include <amdis/common/SharedPtr.hpp>
#include <amdis/common/TupleUtility.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/GridFunctions.hpp>
#include <amdis/gridfunctions/DiscreteFunction.hpp>
#include <amdis/io/FileWriterBase.hpp>
#include <amdis/typetree/Concepts.hpp>
#include <amdis/typetree/TreePath.hpp>

namespace AMDiS
{
  // forward declaration
  template <class Traits>
  class ProblemInstat;

  template <class Traits>
  class ProblemStat
      : public ProblemStatBase
      , public StandardProblemIterationAdaptor<ProblemStat<Traits>>
  {
    using Self = ProblemStat;

    friend class ProblemInstat<Traits>;

  public: // typedefs and static constants

    using GlobalBasis = typename Traits::GlobalBasis;
    using GridView    = typename GlobalBasis::GridView;
    using Grid        = AdaptiveGrid_t<typename GridView::Grid>;
    using Element     = typename GridView::template Codim<0>::Entity;
    using WorldVector = typename Element::Geometry::GlobalCoordinate;
    using WorldMatrix = FieldMatrix<typename WorldVector::field_type, WorldVector::dimension, WorldVector::dimension>;

    static constexpr int dim = Grid::dimension;

    static constexpr int dow = Grid::dimensionworld;

  private:
    using LinAlgTraits = typename Traits::LinAlgTraits;
    using Mat = typename LinAlgTraits::template MatrixImpl<typename Traits::CoefficientType>;
    using Vec = typename LinAlgTraits::template VectorImpl<typename Traits::CoefficientType>;
    using PartitionSet = typename LinAlgTraits::PartitionSet;

  public:
    using LinearSolver = LinearSolverInterface<Mat,Vec>;
    using SystemMatrix = BiLinearForm<GlobalBasis, GlobalBasis, typename Traits::CoefficientType, LinAlgTraits>;
    using SystemVector = LinearForm<GlobalBasis, typename Traits::CoefficientType, LinAlgTraits>;
    using SolutionVector = DOFVector<GlobalBasis, typename Traits::CoefficientType, LinAlgTraits>;

  public:
    explicit ProblemStat(std::string const& name)
      : name_(name)
    {}

    template <class Grid_>
    ProblemStat(std::string const& name, Grid_&& grid)
      : ProblemStat(name)
    {
      adoptGrid(wrap_or_share(FWD(grid)));
    }

    template <class Grid_, class Basis_, class B_ = Underlying_t<Basis_>,
      REQUIRES(Concepts::GlobalBasis<B_>)>
    ProblemStat(std::string const& name, Grid_&& grid, Basis_&& globalBasis)
      : ProblemStat(name, FWD(grid))
    {
      adoptGlobalBasis(wrap_or_share(FWD(globalBasis)));
    }

    template <class Grid_, class PBF_, class GV_ = typename Underlying_t<Grid_>::LeafGridView,
      REQUIRES(Concepts::PreBasisFactory<PBF_, GV_, MultiIndex_t<PBF_>>)>
    ProblemStat(std::string const& name, Grid_&& grid, PBF_ const& preBasisFactory)
      : ProblemStat(name, FWD(grid))
    {
      adoptGlobalBasis(makeSharedPtr(GlobalBasis{grid_->leafGridView(), preBasisFactory}));
    }


    void initialize(Flag initFlag, Self* adoptProblem = nullptr, Flag adoptFlag = INIT_NOTHING);


    void restore(Flag initFlag);




    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath row = {}, ColTreePath col = {})
    {
      static constexpr bool isValidTreePath =
        Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, RowTreePath> &&
        Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, ColTreePath>;
      static_assert(isValidTreePath, "Invalid row and/or col treepath passed to addMatrixOperator!");

      if constexpr (isValidTreePath)
        systemMatrix_->addOperator(tag::element_operator<Element>{}, op, row, col);
    }


    template <class Operator, class RowTreePath = RootTreePath, class ColTreePath = RootTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath row = {}, ColTreePath col = {})
    {
      using I = typename GridView::Intersection;
      static constexpr bool isValidTreePath =
        Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, RowTreePath> &&
        Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, ColTreePath>;
      static_assert(isValidTreePath, "Invalid row and/or col treepath passed to addMatrixOperator!");

      if constexpr (isValidTreePath)
        systemMatrix_->addOperator(tag::boundary_operator<I>{*boundaryManager_,b}, op, row, col);
    }


    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(Operator const& op, TreePath path = {})
    {
      static constexpr bool isValidTreePath =
        Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, TreePath>;
      static_assert(isValidTreePath, "Invalid treepath passed to addVectorOperator!");

      if constexpr (isValidTreePath)
        rhs_->addOperator(tag::element_operator<Element>{}, op, path);
    }


    template <class Operator, class TreePath = RootTreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath path = {})
    {
      using I = typename GridView::Intersection;
      static constexpr bool isValidTreePath =
        Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, TreePath>;
      static_assert(isValidTreePath, "Invalid treepath passed to addVectorOperator!");

      if constexpr (isValidTreePath)
        rhs_->addOperator(tag::boundary_operator<I>{*boundaryManager_,b}, op, path);
    }


    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
    void addDirichletBC(Predicate const& predicate,
                        RowTreePath row, ColTreePath col,
                        Values const& values);

    template <class RowTreePath, class ColTreePath, class Values>
    void addDirichletBC(BoundaryType id,
                        RowTreePath row, ColTreePath col,
                        Values const& values);

    template <class Identifier, class Values>
    void addDirichletBC(Identifier&& id, Values&& values)
    {
      addDirichletBC(FWD(id), RootTreePath{}, RootTreePath{}, FWD(values));
    }

    void addPeriodicBC(BoundaryType id, WorldMatrix const& A, WorldVector const& b);
    void addConstraint(BoundaryCondition<SystemMatrix, SolutionVector, SystemVector> constraint)
    {
      constraints_.push_back(constraint);
    }

  public:

    Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
    {
      return StandardProblemIteration::oneIteration(adaptInfo, toDo);
    }

    void buildAfterAdapt(AdaptInfo& adaptInfo,
                         Flag flag,
                         bool asmMatrix = true,
                         bool asmVector = true) override;

    void assemble(AdaptInfo& adaptInfo)
    {
      buildAfterAdapt(adaptInfo, Flag{0}, true, true);
    }

    void solve(AdaptInfo& adaptInfo,
               bool createMatrixData = true,
               bool storeMatrixData = false) override;

    void estimate(AdaptInfo& /*adaptInfo*/) override { /* do nothing. */ }

    Flag adaptGrid(AdaptInfo& adaptInfo) override;

    Flag markElements(AdaptInfo& adaptInfo) override;

    Flag globalCoarsen(int n) override;

    Flag globalRefine(int n) override;

    void writeFiles(AdaptInfo& adaptInfo, bool force = false);


  public: // get-methods

    std::string const& name() const override { return name_; }

    std::shared_ptr<Grid>       grid()       { return grid_; }
    std::shared_ptr<Grid const> grid() const { return grid_; }

    GridView gridView() const { return globalBasis_->gridView(); }

    std::shared_ptr<BoundaryManager<Grid>>       boundaryManager()       { return boundaryManager_; }
    std::shared_ptr<BoundaryManager<Grid> const> boundaryManager() const { return boundaryManager_; }

    std::shared_ptr<GlobalBasis>       globalBasis()       { return globalBasis_; }
    std::shared_ptr<GlobalBasis const> globalBasis() const { return globalBasis_; }

    std::shared_ptr<LinearSolver>       solver()       { return linearSolver_; }
    std::shared_ptr<LinearSolver const> solver() const { return linearSolver_; }

    std::shared_ptr<SystemMatrix>       systemMatrix()       { return systemMatrix_; }
    std::shared_ptr<SystemMatrix const> systemMatrix() const { return systemMatrix_; }

    std::shared_ptr<SolutionVector>       solutionVector()       { return solution_; }
    std::shared_ptr<SolutionVector const> solutionVector() const { return solution_; }

    std::shared_ptr<SystemVector>       rhsVector()       { return rhs_; }
    std::shared_ptr<SystemVector const> rhsVector() const { return rhs_; }



    template <class Range = void, class... Indices>
    auto solution(Indices... ii)
    {
      assert(bool(solution_) && "You have to call initialize() before.");
      return valueOf<Range>(*solution_, ii...);
    }


    template <class Range = void, class... Indices>
    auto solution(Indices... ii) const
    {
      assert(bool(solution_) && "You have to call initialize() before.");
      return valueOf<Range>(*solution_, ii...);
    }


  public: // set-methods

    template <class Solver_>
    void setSolver(Solver_&& solver)
    {
      linearSolver_ = wrap_or_share(FWD(solver));
    }


    template <class Grid_>
    void setGrid(Grid_&& grid)
    {
      adoptGrid(wrap_or_share(FWD(grid)));
      createGlobalBasis();
      createMatricesAndVectors();
      createMarker();
      createFileWriter();
    }



    template <class Marker_>
    void addMarker(Marker_&& m)
    {
      auto marker = wrap_or_share(FWD(m));
      auto it = marker_.emplace(marker->name(), marker);
      if (marker_.size() > 1)
        it.first->second->setMaximumMarking(true);
    }

    void removeMarker(std::string name)
    {
      std::size_t num = marker_.erase(name);
      test_warning(num == 1, "A marker with the given name '{}' does not exist.", name);
    }

    void removeMarker(Marker<Grid> const& marker)
    {
      removeMarker(marker.name());
    }


  protected: // initialization methods

    void createGlobalBasis();
    void createGrid();
    void createMatricesAndVectors();
    void createSolver();
    void createMarker();
    void createFileWriter();

    void adoptGlobalBasis(std::shared_ptr<GlobalBasis> globalBasis)
    {
      globalBasis_ = std::move(globalBasis);
      initGlobalBasis();
    }

    void adoptGrid(std::shared_ptr<Grid> const& grid,
                   std::shared_ptr<BoundaryManager<Grid>> const& boundaryManager)
    {
      grid_ = grid;
      boundaryManager_ = boundaryManager;
      Parameters::get(name_ + "->mesh", gridName_);
    }

    void adoptGrid(std::shared_ptr<Grid> const& grid)
    {
      adoptGrid(grid, std::make_shared<BoundaryManager<Grid>>(grid));
    }

    void adoptGrid(std::shared_ptr<typename Grid::HostGrid> const& hostGrid)
    {
      auto grid = std::make_shared<Grid>(hostGrid);
      adoptGrid(grid, std::make_shared<BoundaryManager<Grid>>(grid));
    }

  private:

    void createGlobalBasisImpl(std::true_type);
    void createGlobalBasisImpl(std::false_type);

    void initGlobalBasis();

  protected:

    std::string name_;

    std::shared_ptr<Grid> grid_;

    std::string gridName_ = "mesh";

    std::shared_ptr<BoundaryManager<Grid>> boundaryManager_;

    std::shared_ptr<GlobalBasis> globalBasis_;

    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;

    std::map<std::string, std::shared_ptr<Marker<Grid>>> marker_;

//    std::vector<Estimator*> estimator;

    std::shared_ptr<LinearSolver> linearSolver_;

    std::shared_ptr<SystemMatrix> systemMatrix_;

    std::shared_ptr<SolutionVector> solution_;

    std::shared_ptr<SystemVector> rhs_;

    std::map<std::string, std::vector<double>> estimates_;

    std::list<BoundaryCondition<SystemMatrix, SolutionVector, SystemVector>> constraints_;
  };


  namespace Impl
  {
    template <class Grid, class B, class = void>
    struct DeducedProblemTraits;

    template <class Grid, class PB>
    struct DeducedProblemTraits<Grid,GlobalBasis<PB>,void>
    {
      using type = DefaultProblemTraits<GlobalBasis<PB>>;
    };

    template <class G, class PBF>
    struct DeducedProblemTraits<G,PBF,
      std::enable_if_t<Concepts::PreBasisFactory<PBF, typename G::LeafGridView, MultiIndex_t<PBF>>>>
    {
      using Grid = AdaptiveGrid_t<G>;
      using GridView = typename Grid::LeafGridView;
      using Basis = decltype(GlobalBasis{std::declval<GridView>(),std::declval<PBF>()});

      using type = DefaultProblemTraits<Basis>;
    };

    template <class Grid, class Basis>
    using DeducedProblemTraits_t = typename DeducedProblemTraits<Grid,Basis>::type;
  }


  // Deduction guide
  template <class Grid, class Basis>
  ProblemStat(std::string name, Grid&& grid, Basis&& globalBasis)
    -> ProblemStat<Impl::DeducedProblemTraits_t<Underlying_t<Grid>,Underlying_t<Basis>>>;


  // mark templates as explicitly instantiated in cpp file
  extern template class ProblemStat<YaspGridBasis<2,1>>;
  extern template class ProblemStat<YaspGridBasis<2,2>>;

} // end namespace AMDiS

#include "ProblemStat.inc.hpp"