ProblemStat.inc.hpp

#pragma once

#include <dune/amdis/AdaptInfo.hpp>
#include <dune/amdis/common/Loops.hpp>
#include <dune/amdis/common/Timer.hpp>

namespace AMDiS
{
  template <class Traits>
  void ProblemStatSeq<Traits>::initialize(Flag initFlag,
                                          Self* adoptProblem,
                                          Flag adoptFlag)
  {
    // === create meshes ===
    if (mesh) {
      warning("mesh already created");
    }
    else {
      if (initFlag.isSet(CREATE_MESH) ||
          (!adoptFlag.isSet(INIT_MESH) &&
          (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) {
        createMesh();
      }

      if (adoptProblem &&
          (adoptFlag.isSet(INIT_MESH) ||
          adoptFlag.isSet(INIT_SYSTEM) ||
          adoptFlag.isSet(INIT_FE_SPACE))) {
        mesh = adoptProblem->getMesh();
      }

      componentMeshes.resize(nComponents, mesh.get());
    }

    if (!mesh)
      warning("no mesh created");


    int globalRefinements = 0;
    Parameters::get(meshName + "->global refinements", globalRefinements);
    if (globalRefinements > 0) {
      mesh->globalRefine(globalRefinements);
      msg("After refinement:");
      msg("#elements = "   , mesh->size(0));
      msg("#faces/edges = ", mesh->size(1));
      msg("#vertices = "   , mesh->size(dim));
      msg("");
    }



    // === create fespace ===
    if (feSpaces) {
      warning("feSpaces already created");
    }
    else {
      if (initFlag.isSet(INIT_FE_SPACE) ||
          (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
        createFeSpaces();
      }

      if (adoptProblem &&
          (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
        feSpaces = adoptProblem->getFeSpaces();
      }
    }

    if (!feSpaces)
      warning("no feSpaces created\n");

    msg("#DOFs = ", this->getFeSpace<0>().size());

    // === create system ===
    if (initFlag.isSet(INIT_SYSTEM))
      createMatricesAndVectors();

    if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
      solution = adoptProblem->getSolution();
      rhs = adoptProblem->getRhs();
      systemMatrix = adoptProblem->getSystemMatrix();
    }

    // === create solver ===
    if (linearSolver) {
      warning("solver already created\n");
    }
    else {
      if (initFlag.isSet(INIT_SOLVER))
        createSolver();

      if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
        test_exit(!linearSolver, "solver already created\n");
        linearSolver = adoptProblem->getSolver();
      }
    }

    if (!linearSolver) {
      warning("no solver created\n");
    }

    // === create file writer ===
    if (initFlag.isSet(INIT_FILEWRITER))
      createFileWriter();

    solution->compress();
  }


  template <class Traits>
    template <class OperatorType>
  void ProblemStatSeq<Traits>::
  addMatrixOperator(OperatorType const& op, int i, int j, double* factor, double* estFactor)
  {
    matrixOperators[i][j].emplace_back(std::make_shared<OperatorType>(op), factor, estFactor);
    matrixChanging[i][j] = true;
  }


  template <class Traits>
  void ProblemStatSeq<Traits>::
  addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j, double* factor, double* estFactor)
  {
    matrixBoundaryOperators[i][j].emplace_back(std::make_shared<IntersectionOperator>(op), factor, estFactor, b);
    matrixChanging[i][j] = true;
  }

  template <class Traits>
  void ProblemStatSeq<Traits>::
  addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops)
  {
    for (auto op : ops){
      int r = op.first.first;
      int c = op.first.second;
      matrixOperators[r][c].emplace_back(std::make_shared<ElementOperator>(op.second));
      matrixChanging[r][c] = true;
    }
  }

  template <class Traits>
  void ProblemStatSeq<Traits>::
  addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> > ops)
  {
    for (auto op : ops) {
      int r = op.first.first;
      int c = op.first.second;
      matrixOperators[r][c].emplace_back(std::make_shared<ElementOperator>(op.second.first));
      matrixChanging[r][c] = matrixChanging[r][c] || op.second.second;
    }
  }

  template <class Traits>
    template <class OperatorType>
  void ProblemStatSeq<Traits>::
  addVectorOperator(OperatorType const& op, int i, double* factor, double* estFactor)
  {
    vectorOperators[i].emplace_back(std::make_shared<OperatorType>(op), factor, estFactor);
    vectorChanging[i] = true;
  }

  template <class Traits>
  void ProblemStatSeq<Traits>::
  addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i, double* factor, double* estFactor)
  {
    vectorBoundaryOperators[i].emplace_back(std::make_shared<IntersectionOperator>(op), factor, estFactor, b);
    vectorChanging[i] = true;
  }


  template <class Traits>
  void ProblemStatSeq<Traits>::
  addVectorOperator(std::map< int, ElementOperator> ops)
  {
    for (auto op : ops) {
      vectorOperators[op.first].emplace_back(std::make_shared<ElementOperator>(op.second));
      vectorChanging[op.first] = true;
    }
  }

  template <class Traits>
  void ProblemStatSeq<Traits>::
  addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops)
  {
    for (auto op : ops) {
      vectorOperators[op.first].emplace_back(std::make_shared<ElementOperator>(op.second.first));
      vectorChanging[op.first] = vectorChanging[op.first] || op.second.second;
    }
  }

  /// Adds a Dirichlet boundary condition
  template <class Traits>
    template <class Predicate, class Values>
  void ProblemStatSeq<Traits>::
  addDirichletBC(Predicate const& predicate, int row, int col, Values const& values)
  {
    static_assert( Concepts::Callable<Predicate, WorldVector>,
      "Function passed to addDirichletBC for predicate does not model the Callable<WorldVector> concept");

    static_assert( Concepts::Callable<Values, WorldVector>,
      "Function passed to addDirichletBC for values does not model the Callable<WorldVector> concept");

    test_exit(row >= 0 && row < nComponents, "row number out of range: ", row);
    test_exit(col >= 0 && col < nComponents, "col number out of range: ", col);

//     boundaryConditionSet = true;

    using BcType = DirichletBC<WorldVector>;
    dirichletBc[row][col].emplace_back( new BcType{predicate, values} );
  }


  template <class Traits>
  void ProblemStatSeq<Traits>::
  solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
  {
    Timer t;

    SolverInfo solverInfo(name + "->solver");
    solverInfo.setCreateMatrixData(createMatrixData);
    solverInfo.setStoreMatrixData(storeMatrixData);

    solution->compress();
    linearSolver->solve(systemMatrix->getMatrix(),
                        solution->getVector(), rhs->getVector(),
                        solverInfo);

    if (solverInfo.getInfo() > 0) {
      msg("solution of discrete system needed ", t.elapsed(), " seconds");

      if (solverInfo.getAbsResidual() >= 0.0) {
        if (solverInfo.getRelResidual() >= 0.0)
          msg("Residual norm: ||b-Ax|| = ", solverInfo.getAbsResidual(),
                           ", ||b-Ax||/||b|| = ", solverInfo.getRelResidual());
        else
          msg("Residual norm: ||b-Ax|| = ", solverInfo.getAbsResidual());
      }
    }

    if (solverInfo.doBreak()) {
      std::stringstream tol_str;
      if (solverInfo.getAbsTolerance() > 0
          && solverInfo.getAbsResidual() > solverInfo.getAbsTolerance())
        tol_str << "absTol = " << solverInfo.getAbsTolerance() << " ";
      if (solverInfo.getRelTolerance() > 0
          && solverInfo.getRelResidual() > solverInfo.getRelTolerance())
        tol_str << "relTol = " << solverInfo.getRelTolerance() << " ";
      error_exit("Tolerance ", tol_str.str(), " could not be reached!");
    }
  }


  template <class Traits>
  void ProblemStatSeq<Traits>::
  buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag, bool asmMatrix_, bool asmVector_)
  {
    Timer t;

    forEach(range_<0, nComponents>, [this](auto const _r) {
      this->getFeSpace(_r).update(this->leafGridView());
    });


    std::size_t nnz = 0;
    forEach(range_<0, nComponents>, [&,this](auto const _r)
    {
      static const int R = decltype(_r)::value;
      msg(this->getFeSpace(_r).size(), " DOFs for FeSpace[", R, "]");

      bool asmVector = asmVector_ && (!vectorAssembled[R] || vectorChanging[R]);

      if (asmVector) {
	      rhs->compress(_r);
        rhs->getVector(_r) = 0.0;
      }

      forEach(range_<0, nComponents>, [&,this](auto const _c)
      {
        static const int C = decltype(_c)::value;

        using MatrixData = typename ProblemStatSeq<Traits>::template MatrixData<R, C>;
        using VectorData = typename ProblemStatSeq<Traits>::template VectorData<R>;

        // The DOFMatrix which should be assembled
        auto& dofmatrix    = (*systemMatrix)(_r, _c);
        auto& solution_vec = (*solution)[_c];
        auto& rhs_vec      = (*rhs)[_r];

        bool assembleMatrix = asmMatrix_ && (!matrixAssembled[R][C] || matrixChanging[R][C]);
        bool assembleVector = asmVector  && R == C;

        if (assembleMatrix) {
          // init boundary condition
          for (int c = 0; c < nComponents; ++c)
            for (auto bc : dirichletBc[R][c])
              bc->init(c == C, dofmatrix, solution_vec, rhs_vec);
        }

        auto mat = MatrixData{dofmatrix, matrixOperators[R][C], matrixBoundaryOperators[R][C], matrixIntersectionOperators[R][C], assembleMatrix};
        auto vec = VectorData{rhs_vec, vectorOperators[R], vectorBoundaryOperators[R], vectorIntersectionOperators[R], assembleVector};

        // assemble the DOFMatrix block and the corresponding rhs vector, of r==c

        dofmatrix.init(assembleMatrix);
        this->assemble(mat, vec, this->leafGridView());
        dofmatrix.finish();

        if (assembleMatrix)
          matrixAssembled[R][C] = true;
        if (assembleVector)
          vectorAssembled[R] = true;

        if (assembleMatrix) {
          // finish boundary condition
          for (int c = 0; c < nComponents; ++c)
            for (auto bc : dirichletBc[R][c])
              bc->finish(c == C, dofmatrix, solution_vec, rhs_vec);

          nnz += dofmatrix.getNnz();
        }
      });
    });

    msg("fillin of assembled matrix: ", nnz);

    msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
  }


  template <class Traits>
  void ProblemStatSeq<Traits>::writeFiles(AdaptInfo& adaptInfo, bool force)
  {
    Timer t;
    filewriter->write(adaptInfo.getTime(), *solution);
    msg("writeFiles needed ", t.elapsed(), " seconds");
  }


  template <class Traits>
    template <class LhsData, class RhsData, class GV>
  void ProblemStatSeq<Traits>::assemble(LhsData lhs, RhsData rhs, GV const& gridView)
  {
    auto const& rowFeSpace = lhs.matrix.getRowFeSpace();
    auto const& colFeSpace = lhs.matrix.getColFeSpace();

    if ((lhs.operators.empty() || !lhs.assemble) &&
        (rhs.operators.empty() || !rhs.assemble))
      return; // nothing to do

    for (auto scaledOp : lhs.operators)
      scaledOp.op->init(rowFeSpace, colFeSpace);
    for (auto scaledOp : rhs.operators)
      scaledOp.op->init(rowFeSpace, colFeSpace);

    auto rowLocalView = rowFeSpace.localView();
    auto rowIndexSet  = rowFeSpace.localIndexSet();

    auto colLocalView = colFeSpace.localView();
    auto colIndexSet  = colFeSpace.localIndexSet();

    for (auto const& element : elements(gridView)) {
      rowLocalView.bind(element);
      rowIndexSet.bind(rowLocalView);

      colLocalView.bind(element);
      colIndexSet.bind(colLocalView);

      if (lhs.assemble) {
        ElementMatrix elementMatrix;
        bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix,
                                    lhs.operators, lhs.boundary_operators, lhs.intersection_operators);
        if (add)
          addElementMatrix(lhs.matrix, rowIndexSet, colIndexSet, elementMatrix);
      }

      if (rhs.assemble) {
        ElementVector elementVector;
        bool add = getElementVector(rowLocalView, elementVector,
                                    rhs.operators, rhs.boundary_operators, rhs.intersection_operators);
        if (add)
          addElementVector(rhs.vector, rowIndexSet, elementVector);
      }
    }
  }


  template <class Traits>
    template <class RowView, class ColView>
  bool ProblemStatSeq<Traits>::getElementMatrix(RowView const& rowLocalView,
                                                ColView const& colLocalView,
                                                ElementMatrix& elementMatrix,
                                                std::list<Scaled<ElementOperator>>& operators,
                                                std::list<Scaled<IntersectionOperator>>& boundary_operators,
                                                std::list<Scaled<IntersectionOperator>>& intersection_operators)
  {
    auto const nRows = rowLocalView.tree().finiteElement().size();
    auto const nCols = colLocalView.tree().finiteElement().size();

    auto const& element = rowLocalView.element();
    auto const& gridView = rowLocalView.globalBasis().gridView();

    // fills the entire matrix with zeroes
    elementMatrix.change_dim(nRows, nCols);
    set_to_zero(elementMatrix);

    bool add = false;

    auto assemble_operators = [&](auto& e, auto& operator_list) {
      for (auto scaled : operator_list) {
        bool add_op = scaled.op->getElementMatrix(e, rowLocalView, colLocalView, elementMatrix, scaled.factor);
        add = add || add_op;
      }
    };

    // assemble element operators
    assemble_operators(element, operators);

    // assemble intersection operators
    if (!boundary_operators.empty() || !intersection_operators.empty()) {
      for (auto const& intersection : intersections(gridView, element)) {
        if (intersection.boundary())
          assemble_operators(intersection, boundary_operators);
        else
          assemble_operators(intersection, intersection_operators); // TODO: assemble intersection operators also on the boundary?
      }
    }

    return add;
  }


  template <class Traits>
    template <class RowView>
  bool ProblemStatSeq<Traits>::getElementVector(RowView const& rowLocalView,
                                                ElementVector& elementVector,
                                                std::list<Scaled<ElementOperator>>& operators,
                                                std::list<Scaled<IntersectionOperator>>& boundary_operators,
                                                std::list<Scaled<IntersectionOperator>>& intersection_operators)
  {
    auto const nRows = rowLocalView.tree().finiteElement().size();

    auto const& element = rowLocalView.element();
    auto const& gridView = rowLocalView.globalBasis().gridView();

    // Set all entries to zero
    elementVector.change_dim(nRows);
    set_to_zero(elementVector);

    bool add = false;

    auto assemble_operators = [&](auto& e, auto& operator_list) {
      for (auto scaled : operator_list) {
        bool add_op = scaled.op->getElementVector(e, rowLocalView, elementVector, scaled.factor);
        add = add || add_op;
      }
    };

    // assemble element operators
    assemble_operators(element, operators);

    // assemble intersection operators
    if (!boundary_operators.empty() || !intersection_operators.empty()) {
      for (auto const& intersection : intersections(gridView, element)) {
        if (intersection.boundary())
          assemble_operators(intersection, boundary_operators);
        else
          assemble_operators(intersection, intersection_operators); // TODO: assemble intersection operators also on the boundary?
      }
    }

    return add;
  }


  template <class Traits>
    template <class Matrix, class RowIndexSet, class ColIndexSet>
  void ProblemStatSeq<Traits>::addElementMatrix(Matrix& dofmatrix,
                                                RowIndexSet const& rowIndexSet,
                                                ColIndexSet const& colIndexSet,
                                                ElementMatrix const& elementMatrix)
  {
    for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
      // The global index of the i−th vertex of the element
      auto const row = rowIndexSet.index(i);
      for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
        // The global index of the j−th vertex of the element
        auto const col = colIndexSet.index(j);
        dofmatrix(row, col) += elementMatrix[i][j];
      }
    }
  }


  template <class Traits>
    template <class Vector, class IndexSet>
  void ProblemStatSeq<Traits>::addElementVector(Vector& dofvector,
                                                IndexSet const& indexSet,
                                                ElementVector const& elementVector)
  {
    for (std::size_t i = 0; i < size(elementVector); ++i) {
      // The global index of the i-th vertex
      auto const row = indexSet.index(i);
      dofvector[row] += elementVector[i];
    }
  }

} // end namespace AMDiS