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>
  void ProblemStatSeq<Traits>::addMatrixOperator(shared_ptr<OperatorType> op,
                                                 int i, int j,
                                                 double* factor,
                                                 double* estFactor)
  {
    // TODO: currently the factors are ignored
    assert( estFactor == NULL );

    matrixOperators[{i,j}].push_back(op);
    matrixFactors[{i,j}].push_back(factor);
    matrixChanging[{i,j}] = true;
  }


  template <class Traits>
  void ProblemStatSeq<Traits>::addMatrixOperator(OperatorType& op,
                                                 int i, int j,
                                                 double* factor,
                                                 double* estFactor)
  {
    addMatrixOperator(make_shared<OperatorType>(op), i, j, factor, estFactor);
  }

  template <class Traits>
  void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, shared_ptr<OperatorType> > ops)
  {
    for (auto op : ops)
      addMatrixOperator(op.second, op.first.first, op.first.second);
  }

  template <class Traits>
  void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, std::pair<shared_ptr<OperatorType>,bool> > ops)
  {
    for (auto op : ops) {
      auto row_col = op.first;
      matrixOperators[row_col].push_back(op.second.first);
      matrixFactors[row_col].push_back(NULL);
      matrixChanging[row_col] = matrixChanging[row_col] || op.second.second;
    }
  }

  template <class Traits>
  void ProblemStatSeq<Traits>::addVectorOperator(shared_ptr<OperatorType> op,
                                                 int i,
                                                 double* factor,
                                                 double* estFactor)
  {
    // TODO: currently the est-factors are ignored
    assert( estFactor == NULL );

    vectorOperators[i].push_back(op);
    vectorFactors[i].push_back(factor);
    vectorChanging[i] = true;
  }

  template <class Traits>
  void ProblemStatSeq<Traits>::addVectorOperator(OperatorType& op,
                                                 int i,
                                                 double* factor,
                                                 double* estFactor)
  {
    addVectorOperator(make_shared<OperatorType>(op), i, factor, estFactor);
  }

  template <class Traits>
  void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, shared_ptr<OperatorType> > ops)
  {
    for (auto op : ops)
      addVectorOperator(op.second, op.first);
  }

  template <class Traits>
  void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, std::pair<shared_ptr<OperatorType>, bool> > ops)
  {
    for (auto op : ops) {
      vectorOperators[op.first].push_back(op.second.first);
      vectorFactors[op.first].push_back(NULL);
      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;

    For<0, nComponents>::loop([this](auto const _r) {
	    this->getNodeFactory(_r).initializeIndices();
    });


    size_t nnz = 0;
    For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector_](auto const r)
    {
      msg(this->getFeSpace(r).size(), " DOFs for FeSpace[", r, "]");

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

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

      For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector](auto const c)
      {
        static const int _c = decltype(c)::value;
        static const int _r = decltype(r)::value;
        index_<_r> r{};

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

        auto row_col = std::make_pair(_r, _c);
        bool assembleMatrix = asmMatrix_ && (!matrixAssembled[row_col] || matrixChanging[row_col]);
        bool assembleVector = asmVector  && _r == _c;

        int r_ = 0, c_ = 0;
        if (assembleMatrix) {
          // init boundary condition
          for (auto bc_list : dirichletBc) {
            std::tie(r_, c_) = bc_list.first;
            if (r_ == _r) {
                    for (auto bc : bc_list.second)
                      bc->init(c_ == _c, dofmatrix, solution_vec, rhs_vec);
            }
          }
        }

        auto mat = MatrixData{dofmatrix, matrixOperators[row_col], matrixFactors[row_col], assembleMatrix};
        auto vec = VectorData{rhs_vec,   vectorOperators[_r], vectorFactors[_r], assembleVector};

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

        dofmatrix.init(assembleMatrix);
        this->assemble(mat, vec, elements(*meshView));
        dofmatrix.finish();

        if (assembleMatrix)
          matrixAssembled[row_col] = true;
        if (assembleVector)
          vectorAssembled[_r] = true;

        if (assembleMatrix) {
          // finish boundary condition
          for (auto bc_list : dirichletBc) {
            std::tie(r_, c_) = bc_list.first;
            if (r_ == _r) {
                    for (auto bc : bc_list.second)
                      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 Elements>
  void ProblemStatSeq<Traits>::assemble(LhsData lhs, RhsData rhs,
                                        Elements const& elements)
  {
    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 op : lhs.operators)
      op->init(rowFeSpace, colFeSpace);
    for (auto op : rhs.operators)
      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) {
      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.factors);
        if (add)
          addElementMatrix(lhs.matrix, rowIndexSet, colIndexSet, elementMatrix);
      }

      if (rhs.assemble) {
        ElementVector elementVector;
        bool add = getElementVector(rowLocalView, elementVector,
                                    rhs.operators, rhs.factors);
        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<shared_ptr<OperatorType>>& operators,
                                                std::list<double*> const& factors)
  {
    auto const nRows = rowLocalView.tree().finiteElement().size();
    auto const nCols = colLocalView.tree().finiteElement().size();

    test_exit_dbg(operators.size() == factors.size(),
                  "Nr. of operators and factors must be consistent!");

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

    bool add = false;

    auto op_it = operators.begin();
    auto fac_it = factors.begin();
    for (; op_it != operators.end(); ++op_it, ++fac_it) {
      bool add_op = (*op_it)->getElementMatrix(rowLocalView, colLocalView, elementMatrix, *fac_it);
      add = add || add_op;
    }

    return add;
  }


  template <class Traits>
    template <class RowView>
  bool ProblemStatSeq<Traits>::getElementVector(RowView const& rowLocalView,
                                                ElementVector& elementVector,
                                                std::list<shared_ptr<OperatorType>>& operators,
                                                std::list<double*> const& factors)
  {
    auto const nRows = rowLocalView.tree().finiteElement().size();

    test_exit_dbg(operators.size() == factors.size(),
                  "Nr. of operators and factors must be consistent!");

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

    bool add = false;

    auto op_it = operators.begin();
    auto fac_it = factors.begin();
    for (; op_it != operators.end(); ++op_it, ++fac_it) {
      bool add_op = (*op_it)->getElementVector(rowLocalView, elementVector, *fac_it);
      add = add || add_op;
    }

    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 (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 (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 (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