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/grid/common/grid.hh>

#include <amdis/AdaptInfo.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/CreatorInterface.hpp>
#include <amdis/DirichletBC.hpp>
//#include <amdis/Estimator.hpp>
#include <amdis/Flag.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/LinearAlgebra.hpp>
#include <amdis/Marker.cpp>
#include <amdis/Mesh.hpp>
#include <amdis/ProblemStatBase.hpp>
#include <amdis/ProblemStatTraits.hpp>
#include <amdis/StandardProblemIteration.hpp>

#include <amdis/common/TupleUtility.hpp>
#include <amdis/common/TypeDefs.hpp>
#include <amdis/common/Utility.hpp>

#include <amdis/GridFunctions.hpp>
#include <amdis/gridfunctions/DOFVectorView.hpp>

#include <amdis/io/FileWriterInterface.hpp>

#include <amdis/utility/TreePath.hpp>

namespace AMDiS
{
  template <class Traits>
  class ProblemStat
      : public ProblemStatBase
      , public StandardProblemIteration
  {
    using Self = ProblemStat;

  public: // typedefs and static constants

    using GlobalBasis = typename Traits::GlobalBasis;
    using GridView = typename GlobalBasis::GridView;
    using Grid     = typename GridView::Grid;
    using Element  = typename GridView::template Codim<0>::Entity;
    using WorldVector = typename Element::Geometry::GlobalCoordinate;

    /// Dimension of the mesh
    static constexpr int dim = Grid::dimension;

    /// Dimension of the world
    static constexpr int dow = Grid::dimensionworld;

    using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
    using SystemVector = DOFVector<Traits, double>;

    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;

  public:
    /**
     * \brief Constructor. Takes the name of the problem that is used to
     * access values corresponding to this problem in the parameter file.
     **/
    explicit ProblemStat(std::string name)
      : StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
      , name_(std::move(name))
    {}

    /// Constructor taking additionally a reference to a mesh that is used
    /// instead of the default created mesh, \ref ProblemStat
    ProblemStat(std::string name, Grid& grid)
      : ProblemStat(std::move(name))
    {
      grid_ = Dune::stackobject_to_shared_ptr(grid);
      Parameters::get(name_ + "->mesh", gridName_);
    }

    /// \brief Constructor taking a grid reference and a basis reference.
    /// Stores pointers to both.
    ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
      : ProblemStat(std::move(name), grid)
    {
      globalBasis_ = Dune::stackobject_to_shared_ptr(globalBasis);
      initGlobalBasis(*globalBasis_);
    }


    /**
     * \brief Initialisation of the problem.
     *
     * Parameters read in initialize() for problem with name 'PROB'
     *   MESH[0]->global refinements:  nr of initial global refinements
     **/
    void initialize(Flag initFlag,
                    Self* adoptProblem = nullptr,
                    Flag adoptFlag = INIT_NOTHING);


    /// Adds an operator to \ref A.
    /** @{ */
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(Operator const& op, RowTreePath, ColTreePath,
                           double* factor = nullptr, double* estFactor = nullptr);

    // operator evaluated on the boundary
    template <class Operator, class RowTreePath, class ColTreePath>
    void addMatrixOperator(BoundaryType b, Operator const& op, RowTreePath, ColTreePath,
                           double* factor = nullptr, double* estFactor = nullptr);
    /** @} */


    /// Adds an operator to \ref rhs.
    /** @{ */
    template <class Operator, class TreePath>
    void addVectorOperator(Operator const& op, TreePath,
                           double* factor = nullptr, double* estFactor = nullptr);

    // operator evaluated on the boundary
    template <class Operator, class TreePath>
    void addVectorOperator(BoundaryType b, Operator const& op, TreePath,
                           double* factor = nullptr, double* estFactor = nullptr);
    /** @} */


    /// Adds a Dirichlet boundary condition
    template <class Predicate, class RowTreePath, class ColTreePath, class Values>
    void addDirichletBC(Predicate const& predicate,
                        RowTreePath row, ColTreePath col,
                        Values const& values);

  public:

    /// Implementation of \ref ProblemStatBase::solve
    virtual void solve(AdaptInfo& adaptInfo,
                       bool createMatrixData = true,
                       bool storeMatrixData = false) override;


    /// Implementation of \ref ProblemStatBase::estimate.
    virtual void estimate(AdaptInfo& adaptInfo) override;


    /// Implementation of \ref ProblemStatBase::markElements.
    virtual Flag markElements(AdaptInfo& adaptInfo) override;


    /// Implementation of \ref ProblemStatBase::refineMesh.
    virtual Flag refineMesh(AdaptInfo& adaptInfo) override;


    /// Implementation of \ref ProblemStatBase::buildAfterCoarse
    virtual void buildAfterCoarsen(AdaptInfo& adaptInfo,
                                   Flag flag,
                                   bool asmMatrix = true,
                                   bool asmVector = true) override;


    /// Writes output files.
    void writeFiles(AdaptInfo& adaptInfo, bool force = false);


  public: // get-methods

    /// Returns a pointer to system-matrix, \ref systemMatrix_
    std::shared_ptr<SystemMatrix> getSystemMatrix()       { return systemMatrix_; }
    std::shared_ptr<SystemMatrix> getSystemMatrix() const { return systemMatrix_; }

    /// Returns a pointer to the solution vector, \ref solution_
    std::shared_ptr<SystemVector> getSolutionVector() const
    {
      return solution_;
    }

    /// Return a mutable view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {})
    {
      auto&& tp = makeTreePath(path);
      return makeDOFVectorView(*solution_, tp);
    }

    /// Return a const view to a solution component
    template <class TreePath = RootTreePath>
    auto getSolution(TreePath const& path = {}) const
    {
      auto&& tp = makeTreePath(path);
      return makeDOFVectorView(*solution_, tp);
    }


    /// Return a point to the rhs system-vector, \ref rhs
    std::shared_ptr<SystemVector> getRhsVector()       { return rhs_; }
    std::shared_ptr<SystemVector> getRhsVector() const { return rhs_; }


    /// Return a pointer to the linear solver, \ref linearSolver
    std::shared_ptr<LinearSolverType> getSolver() { return linearSolver_; }

    void setSolver(std::shared_ptr<LinearSolverType> const& solver)
    {
      linearSolver_ = solver;
    }

    /// Return a pointer to the grid, \ref grid
    std::shared_ptr<Grid> getGrid() { return grid_; }

    /// Set the mesh. Stores pointer to passed reference and initializes feSpaces
    /// matrices and vectors, as well as the file-writer.
    void setGrid(Grid& grid)
    {
      grid_ = Dune::stackobject_to_shared_ptr(grid);

      createGlobalBasis();
      createMatricesAndVectors();
      createFileWriter();
    }


    /// Return the gridView of the leaf-level
    auto leafGridView() { return grid_->leafGridView(); }

    /// Return the gridView of levle `level`
    auto levelGridView(int level) { return grid_->levelGridView(level); }

    /// Return the \ref feSpaces
    std::shared_ptr<GlobalBasis> const& getGlobalBasis() { return globalBasis_; }


    /// Implementation of \ref ProblemStatBase::getName
    virtual std::string getName() const override
    {
      return name_;
    }

  protected: // initialization methods

    void createGrid()
    {
      gridName_ = "";
      Parameters::get(name_ + "->mesh", gridName_);
      test_exit(!gridName_.empty(), "No mesh name specified for '", name_, "->mesh'!");

      grid_ = MeshCreator<Grid>::create(gridName_);

      msg("Create grid:");
      msg("#elements = "   , grid_->size(0));
      msg("#faces/edges = ", grid_->size(1));
      msg("#vertices = "   , grid_->size(dim));
      msg("");
    }

    void createGlobalBasis()
    {
      globalBasis_ = std::make_shared<GlobalBasis>(makeGlobalBasis<GlobalBasis>(grid_->leafGridView()));
      initGlobalBasis(*globalBasis_);
    }

    void initGlobalBasis(GlobalBasis const& globalBasis)
    {
      localView_ = std::make_shared<typename GlobalBasis::LocalView>(globalBasis.localView());
      matrixOperators_.init(localView_->tree(), tag::store{});
      rhsOperators_.init(localView_->tree(), tag::store{});
      constraints_.init(localView_->tree(), tag::store{});
    }

    void createMatricesAndVectors()
    {
      systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_, "mat");
      solution_ = std::make_shared<SystemVector>(*globalBasis_, "solution");
      rhs_ = std::make_shared<SystemVector>(*globalBasis_, "rhs");

      estimates.resize(leafGridView().indexSet().size(0));
      for (std::size_t i = 0; i < estimates.size(); i++) {
        estimates[i].resize(nComponents);
        for (std::size_t j = 0; j < estimates[i].size(); j++) {
          estimates[i][j] = 0.0; // TODO: Remove when estimate() is implemented
        }
      }
    }

    void createSolver()
    {
      std::string solverName = "cg";
      Parameters::get(name_ + "->solver->name", solverName);

      auto solverCreator
        = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));

      linearSolver_ = solverCreator->create(name_ + "->solver");
    }

    void createEstimator()
    {/*
      for (std::size_t i = 0, i < nComponents, i++) {
        std::string estName = "";
        Parameters::get(name + "->estimator->name[" + std::to_string(i) + "]", estName);
        auto estimatorCreator
          = named(CreatorMap<Estimator>::getCreator(estName, name + "->estimator->name[" + std::to_string(i) + "]"));
        if (estimatorCreator) {
	        estimatorCreator->setName(estName);
	        estimatorCreator->setRow(i);
	        estimatorCreator->setSolution(solution->getDOFVector(i));
	        estimator[i] = estimatorCreator->create();
        }

        if (estimator[i]) {
          for (std::size_t j = 0; j < nComponents; j++)
	          estimator[i]->addSystem((*systemMatrix)[i][j],
              solution->getDOFVector(j),
				      rhs->getDOFVector(i));
        }
      }*/
    }


    void createMarker()
    {
      int nMarkersCreated = 0;

      marker.resize(nComponents);

      for (std::size_t i = 0; i < nComponents; i++) {
        marker[i] = Marker<Grid>::
          createMarker(name + "->marker[" + std::to_string(i) + "]", i,
                       estimates, componentGrids[i]);

        if (marker[i]) {
	        nMarkersCreated++;

	        // If there is more than one marker, and all components are defined
	        // on the same mesh, the maximum marking has to be enabled.
 	        if (nMarkersCreated > 1 && nGrids == 1)
 	          marker[i]->setMaximumMarking(true);
        }
      }
    }


    void createFileWriter();


  public: // implementation of iteration interface methods

    /**
      * \brief Determines the execution order of the single adaption steps.
      *
      * If adapt is true, mesh adaption will be performed. This allows to avoid
      * mesh adaption, e.g. in timestep adaption loops of timestep adaptive
      * strategies.
      *
      * Implementation of \ref StandardProblemIteration::oneIteration.
      **/
    virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
    {
      return StandardProblemIteration::oneIteration(adaptInfo, toDo);
    }

    /// Implementation of \ref ProblemStatBase::buildBeforeRefine.
    virtual void buildBeforeRefine(AdaptInfo&, Flag) override { /* Does nothing. */ }

    /// Implementation of \ref ProblemStatBase::buildBeforeCoarsen.
    virtual void buildBeforeCoarsen(AdaptInfo&, Flag) override { /* Does nothing. */ }

    /// Implementation of \ref ProblemStatBase::estimate.
    virtual void estimate(AdaptInfo& adaptInfo) override { /* Does nothing. */ }

    /// Implementation of \ref ProblemStatBase::refineMesh.
    virtual Flag refineMesh(AdaptInfo& adaptInfo) override { return 0; }

    /// Implementation of \ref ProblemStatBase::coarsenMesh.
    virtual Flag coarsenMesh(AdaptInfo& adaptInfo) override { return 0; }

    /// Implementation of \ref ProblemStatBase::markElements.
    virtual Flag markElements(AdaptInfo& adaptInfo) override { return 0; }


  private:
    /// Name of this problem.
    std::string name_;

    /// Grid of this problem.
    std::shared_ptr<Grid> grid_; // TODO: generalize to multi-mesh problems

    /// Number of grids
    int nGrids = 1;

    /// Name of the mesh
    std::string gridName_ = "none";

    /// FE spaces of this problem.
    std::shared_ptr<GlobalBasis> globalBasis_;
    std::shared_ptr<typename GlobalBasis::LocalView> localView_;

    /// A FileWriter object
    std::list<std::shared_ptr<FileWriterInterface>> filewriter_;

    /// Pointer to the adaptation markers
    std::vector<Marker<Grid>* > marker;

    /// Pointer to the estimators for this problem
//    std::vector<Estimator*> estimator;

    /// An object of the linearSolver Interface
    std::shared_ptr<LinearSolverType> linearSolver_;

    /// A block-matrix that is filled during assembling
    std::shared_ptr<SystemMatrix> systemMatrix_;

    /// A block-vector with the solution components
    std::shared_ptr<SystemVector> solution_;

    /// A vector with the local element error estimates
    std::vector<std::vector<double> > estimates;

    /// A block-vector (load-vector) corresponding to the right.hand side
    /// of the equation, filled during assembling
    std::shared_ptr<SystemVector> rhs_;


  private: // some internal data-structures

    MatrixOperators<GlobalBasis> matrixOperators_;
    VectorOperators<GlobalBasis> rhsOperators_;
    Constraints<GlobalBasis> constraints_;
  };

#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
  // Deduction rule
  template <class Grid, class GlobalBasis>
  ProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
    -> ProblemStat<DefaultProblemTraits<GlobalBasis>>;
#endif

  // Generator for ProblemStat with given Grid and GlobalBasis
  template <class Grid, class GlobalBasis>
  ProblemStat<DefaultProblemTraits<GlobalBasis>>
  makeProblemStat(std::string name, Grid& grid, GlobalBasis& globalBasis)
  {
    return {std::move(name), grid, globalBasis};
  }


#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
  extern template class ProblemStat<YaspGridBasis<2,1,2>>; // 2 components with different polynomial degree
#endif

} // end namespace AMDiS

#include "ProblemStat.inc.hpp"