ProblemStat.hpp

#pragma once

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

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

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

#include <dune/amdis/AdaptInfo.hpp>
#include <dune/amdis/CreatorInterface.hpp>
#include <dune/amdis/DirichletBC.hpp>
#include <dune/amdis/FileWriter.hpp>
#include <dune/amdis/Flag.hpp>
#include <dune/amdis/Initfile.hpp>
#include <dune/amdis/LinearAlgebra.hpp>
#include <dune/amdis/Mesh.hpp>
#include <dune/amdis/Operator.hpp>
#include <dune/amdis/ProblemStatBase.hpp>
#include <dune/amdis/ProblemStatTraits.hpp>
#include <dune/amdis/StandardProblemIteration.hpp>

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

namespace AMDiS
{
  struct BoundaryType { int b; };

  template <class Traits>
  class ProblemStatSeq
      : public ProblemStatBase
  {
    using Self = ProblemStatSeq;

  public: // typedefs and static constants

    using FeSpaces = typename Traits::FeSpaces;
    using Mesh     = typename Traits::Mesh;
    using MeshView = typename Mesh::LeafGridView;

    using Codim0   = typename MeshView::template Codim<0>;
    using Element  = typename Codim0::Entity;

    using ElementVector = Impl::ElementVector;
    using ElementMatrix = Impl::ElementMatrix;


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

    /// Dimension of the world
    static constexpr int dow = Traits::dimworld;

    /// Number of problem components
    static constexpr int nComponents = Traits::nComponents;


    template <std::size_t I>
    using FeSpace = std::tuple_element_t<I, FeSpaces>;

    using WorldVector = typename Codim0::Geometry::GlobalCoordinate;

    using SystemVectorType = SystemVector<FeSpaces>;
    using SystemMatrixType = SystemMatrix<FeSpaces>;

    using ElementOperator = Operator<MeshView, Element>;
    using IntersectionOperator = Operator<MeshView, typename MeshView::Intersection>;

    using LinearSolverType = LinearSolverInterface<typename SystemMatrixType::MultiMatrix,
                                                   typename SystemVectorType::MultiVector>;

    template <class OperatorType>
    struct Scaled
    {
      Scaled(std::shared_ptr<OperatorType> const& op, double* factor = nullptr, double* estFactor = nullptr, BoundaryType b = {0})
        : op(op)
        , factor(factor)
        , estFactor(estFactor)
      {}

      std::shared_ptr<OperatorType> op;
      double* factor;
      double* estFactor;
      BoundaryType b;
    };

  public:
    /**
     * \brief Constructor. Takes the name of the problem that is used to
     * access values correpsonding to this püroblem in the parameter file.
     *
     * Parameters read by ProblemStatSeq, with name 'PROB'
     *   PROB->names: \ref componentNames
     **/
    ProblemStatSeq(std::string name)
      : name(name)
    {
      Parameters::get(name + "->names", componentNames);
      for (std::size_t i = componentNames.size(); i < nComponents; ++i)
        componentNames.push_back("solution[" + std::to_string(i) + "]");
    }

    ProblemStatSeq(std::string name, Mesh& mesh)
      : ProblemStatSeq(name)
    {
      this->mesh = std::shared_ptr<Mesh>(&mesh, optional_delete(false));
      componentMeshes.resize(nComponents, this->mesh.get());

      meshName = "";
      Parameters::get(name + "->mesh", meshName);
    }


    /**
     * \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 OperatorType>
    void addMatrixOperator(OperatorType const& op, int i, int j,
                           double* factor = nullptr, double* estFactor = nullptr);

    // operator evaluated on the boundary
    void addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j,
                           double* factor = nullptr, double* estFactor = nullptr);

    void addMatrixOperator(std::map< std::pair<int,int>,
                                     std::shared_ptr<ElementOperator> > ops);
    void addMatrixOperator(std::map< std::pair<int,int>,
                                     std::pair<std::shared_ptr<ElementOperator>, bool> > ops);
    /** @} */


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

    // operator evaluated on the boundary
    void addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i,
                           double* factor = nullptr, double* estFactor = nullptr);

    void addVectorOperator(std::map< int, std::shared_ptr<ElementOperator> > ops);
    void addVectorOperator(std::map< int, std::pair<std::shared_ptr<ElementOperator>, bool> > ops);
    /** @} */

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


    /// Implementation of \ref ProblemStatBase::solve
    virtual void solve(AdaptInfo& adaptInfo,
                       bool createMatrixData = true,
                       bool storeMatrixData = false) 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 nr of components \ref nComponents
    static constexpr int getNumComponents()
    {
      return nComponents;
    }


    /// Return the \ref systemMatrix
    auto getSystemMatrix()       { return systemMatrix; }
    auto getSystemMatrix() const { return systemMatrix; }


    /// Return the \ref solution
    auto getSolution()       { return solution; }
    auto getSolution() const { return solution; }

    /// Return the i'th \ref solution component
    template <std::size_t I = 0>
    decltype(auto) getSolution(const index_t<I> _i = {})
    {
      return (*solution)[_i];
    }


    /// Return the \ref rhs
    auto getRhs()       { return rhs; }
    auto getRhs() const { return rhs; }


    /// Return the \ref linearSolver
    auto getSolver() { return linearSolver; }

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

    /// Return the \ref mesh
    auto const& getMesh() { return mesh; }

    void setMesh(Mesh& mesh_)
    {
      mesh = std::shared_ptr<Mesh>(&mesh_, optional_delete(false));
      std::fill(componentMeshes.begin(), componentMeshes.end(), this->mesh.get());

      createFeSpaces();
      createMatricesAndVectors();
      createFileWriter();
    }


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

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

    /// Return the \ref feSpaces
    auto getFeSpaces() { return feSpaces; }


    /// Return the I'th \ref feSpaces component
    template <std::size_t I = 0>
    FeSpace<I> const& getFeSpace(const index_t<I> = {}) const
    {
      return std::get<I>(*feSpaces);
    }

    template <std::size_t I = 0>
    FeSpace<I>& getFeSpace(const index_t<I> = {})
    {
      return std::get<I>(*feSpaces);
    }

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

    /// Return a vector of names of the problem components
    std::vector<std::string> getComponentNames() const
    {
      return componentNames;
    }

  protected: // initialization methods

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

      mesh = MeshCreator<Mesh>::create(meshName);

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

    void createFeSpaces()
    {
      feSpaces = std::make_shared<FeSpaces>(constructTuple<FeSpaces>(leafGridView()));
    }

    void createMatricesAndVectors()
    {
      systemMatrix = std::make_shared<SystemMatrixType>(*feSpaces);
      solution = std::make_shared<SystemVectorType>(*feSpaces, componentNames);

      auto rhsNames = std::vector<std::string>(nComponents, "rhs");
      rhs = std::make_shared<SystemVectorType>(*feSpaces, rhsNames);
    }

    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 createFileWriter()
    {
      filewriter = std::make_shared<FileWriter<Traits>>(name + "->output", leafGridView(), componentNames);
    }


  protected: // sub-methods to assemble DOFMatrix

    template <class LhsData, class RhsData, class GV>
    void assemble(LhsData lhs, RhsData rhs, GV const& gridView);

    template <class RowView, class ColView>
    bool getElementMatrix(RowView const& rowView,
                          ColView const& colView,
                          ElementMatrix& elementMatrix,
                          std::list<Scaled<ElementOperator>>& operators,
                          std::list<Scaled<IntersectionOperator>>& boundary_operators,
                          std::list<Scaled<IntersectionOperator>>& intersection_operators);

    template <class RowView>
    bool getElementVector(RowView const& rowView,
                          ElementVector& elementVector,
                          std::list<Scaled<ElementOperator>>& operators,
                          std::list<Scaled<IntersectionOperator>>& boundary_operators,
                          std::list<Scaled<IntersectionOperator>>& intersection_operators);


    template <class Matrix, class RowIndexSet, class ColIndexSet>
    void addElementMatrix(Matrix& matrix,
                          RowIndexSet const& rowIndexSet,
                          ColIndexSet const& colIndexSet,
                          ElementMatrix const& elementMatrix);

    template <class Vector, class IndexSet>
    void addElementVector(Vector& vector,
                          IndexSet const& indexSet,
                          ElementVector const& elementvector);


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

    /// Stores the names for all components. Is used for naming the solution
    /// vectors, \ref solution.
    std::vector<std::string> componentNames;

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

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

    /// Pointer to the meshes for the different problem components
    std::vector<Mesh*> componentMeshes;

    /// FE spaces of this problem.
    std::shared_ptr<FeSpaces> feSpaces; // eventuell const

    /// A FileWriter object
    std::shared_ptr<FileWriter<Traits>> filewriter;

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

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

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

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


    template <class T>
    using MatrixEntries = Dune::FieldMatrix<T, nComponents, nComponents>;

    template <class T>
    using VectorEntries = Dune::FieldVector<T, nComponents>;

    using BoundaryCondition = std::list<std::shared_ptr<DirichletBC<WorldVector>>>;

    /// A map (i,j) -> list<DirichleBC> string a boundary conditions for
    /// each matrix block
    MatrixEntries<BoundaryCondition> dirichletBc;

    /// A map (i,j) -> list<ElementOperator> string the differential operators for
    /// each matrix block
    MatrixEntries<std::list<Scaled<ElementOperator>>> matrixOperators;
    MatrixEntries<std::list<Scaled<IntersectionOperator>>> matrixBoundaryOperators;
    MatrixEntries<std::list<Scaled<IntersectionOperator>>> matrixIntersectionOperators;
    MatrixEntries<bool> matrixAssembled; // if false, do reassemble
    MatrixEntries<bool> matrixChanging;  // if true, or vectorAssembled false, do reassemble

    /// A map (i) -> list<ElementOperator> string the differential operators for
    /// each rhs block
    VectorEntries<std::list<Scaled<ElementOperator>>> vectorOperators;
    VectorEntries<std::list<Scaled<IntersectionOperator>>> vectorBoundaryOperators;
    VectorEntries<std::list<Scaled<IntersectionOperator>>> vectorIntersectionOperators;
    VectorEntries<bool> vectorAssembled; // if false, do reassemble
    VectorEntries<bool> vectorChanging;  // if true, or vectorAssembled false, do reassemble

    template <int R, int C>
    struct MatrixData
    {
      using DOFMatrixType =
        std::tuple_element_t<C, std::tuple_element_t<R, typename SystemMatrixType::DOFMatrices>>;

      DOFMatrixType& matrix;
      std::list<Scaled<ElementOperator>>& operators;
      std::list<Scaled<IntersectionOperator>>& boundary_operators;
      std::list<Scaled<IntersectionOperator>>& intersection_operators;
      bool assemble = true;
    };

    template <int I>
    struct VectorData
    {
      using DOFVectorType = std::tuple_element_t<I, typename SystemVectorType::DOFVectors>;

      DOFVectorType& vector;
      std::list<Scaled<ElementOperator>>& operators;
      std::list<Scaled<IntersectionOperator>>& boundary_operators;
      std::list<Scaled<IntersectionOperator>>& intersection_operators;
      bool assemble = true;
    };
  };


#ifndef AMDIS_NO_EXTERN_PROBLEMSTAT
//   extern template class ProblemStatSeq<ProblemStatTraits<1>>;
  extern template class ProblemStatSeq<ProblemStatTraits<2,2,1>>;
  extern template class ProblemStatSeq<ProblemStatTraits<2,2,2>>;
//   extern template class ProblemStatSeq<ProblemStatTraits<3>>;
#endif

  namespace Impl
  {
    template <class ProblemStatType>
    struct ProblemStat
        : public ProblemStatType
        , public StandardProblemIteration
    {
      using ProblemStatType::getName;

      /// Constructor
      explicit ProblemStat(std::string nameStr)
        : ProblemStatType(nameStr)
        , StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
      {}

      template <class Grid>
      ProblemStat(std::string nameStr, Grid& grid)
        : ProblemStatType(nameStr, grid)
        , StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
      {}

      /**
       * \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
      {
        for (std::size_t i = 0; i < ProblemStatType::getNumComponents(); ++i)
          if (adaptInfo.spaceToleranceReached(i))
            adaptInfo.allowRefinement(false, i);
          else
            adaptInfo.allowRefinement(true, i);

        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}; }
    };

  } // end namespace Impl

  template <class Traits>
  using ProblemStat = Impl::ProblemStat<ProblemStatSeq<Traits>>;

} // end namespace AMDiS

#include "ProblemStat.inc.hpp"