ProblemStat.inc.hpp

#pragma once

#include "AdaptInfo.hpp"
#include "Loops.hpp"
#include "Timer.hpp"

namespace AMDiS
{
  template <class Traits>
  void ProblemStatSeq<Traits>::initialize(Flag initFlag,
                                          Self* adoptProblem,
                                          Flag adoptFlag) 
  {
    // === create meshes ===
    if (mesh) {
      AMDIS_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)
      AMDIS_WARNING("no mesh created");
    
    
    int globalRefinements = 0;
    Parameters::get(meshName + "->global refinements", globalRefinements);
    if (globalRefinements > 0) {
      mesh->globalRefine(globalRefinements);
      AMDIS_MSG("After refinement:");
      AMDIS_MSG("#elements = "    << mesh->size(0));
      AMDIS_MSG("#faces/edges = " << mesh->size(1));
      AMDIS_MSG("#vertices = "    << mesh->size(dim));
      AMDIS_MSG("");
    }
    
      
    
    // === create fespace ===
    if (feSpaces) {
      AMDIS_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)
      AMDIS_WARNING("no feSpaces created\n");
    
    AMDIS_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) {
      AMDIS_WARNING("solver already created\n");
    }
    else {
      if (initFlag.isSet(INIT_SOLVER))
        createSolver();

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

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

    // === create file writer ===
    if (initFlag.isSet(INIT_FILEWRITER))
      createFileWriter();
    
    solution->compress();
  }


  template <class Traits>
  void ProblemStatSeq<Traits>::addMatrixOperator(OperatorType& op, 
                                                 int i, int j,
                                                 double* factor, 
                                                 double* estFactor)
  {
    // TODO: currently the factors are ignored
    assert( estFactor == NULL );
    
    matrixOperators[{i,j}].emplace_back(new OperatorType{op});
    matrixFactors[{i,j}].push_back(factor);
  }

  
  template <class Traits>
  void ProblemStatSeq<Traits>::addVectorOperator(OperatorType& op, 
                                                 int i,
                                                 double* factor, 
                                                 double* estFactor)
  {
    // TODO: currently the est-factors are ignored
    assert( estFactor == NULL );
    
    vectorOperators[i].emplace_back(new OperatorType{op});
    vectorFactors[i].push_back(factor);
  }

  
  
  /// 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(Dune::Functions::Concept::isCallable<Predicate, WorldVector>(), 
      "Function passed to addDirichletBC for predicate does not model the Callable<WorldVector> concept");
    
    static_assert(Dune::Functions::Concept::isCallable<Values, WorldVector>(), 
      "Function passed to addDirichletBC for values does not model the Callable<WorldVector> concept");

    AMDIS_TEST_EXIT(row >= 0 && row < nComponents, "row number out of range: " << row);
    AMDIS_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);
    
    AMDIS_MSG("solution of discrete system needed " << t.elapsed() << " seconds");
    
    if (solverInfo.getInfo() > 0) {
      if (solverInfo.getAbsResidual() >= 0.0) {
        if (solverInfo.getRelResidual() >= 0.0)
          AMDIS_MSG("Residual norm: ||b-Ax|| = " << solverInfo.getAbsResidual()
                        << ", ||b-Ax||/||b|| = " << solverInfo.getRelResidual());
        else
          AMDIS_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() << " ";
      AMDIS_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) 
    {
      AMDIS_MSG(this->getFeSpace(_r).size() << " DOFs for FeSpace[" << _r << "]");
      
      if (asmVector) {
	rhs->compress(_r);
        rhs->getVector(_r) = 0.0;
      }

      For<0, nComponents>::loop([this, &nnz, asmMatrix, asmVector, _r](auto const _c) 
      {
        // The DOFMatrix which should be assembled
        auto& dofmatrix    = (*systemMatrix)(_r, _c);
        auto& solution_vec = (*solution)[_c];
        auto& rhs_vec      = (*rhs)[_r];
	  
        int r = 0, c = 0;
	if (asmMatrix) {

	  // init boundary condition
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
              for (auto bc : bc_list.second)
                bc->init(c == int(_c), dofmatrix, solution_vec, rhs_vec);
	    }
	  }
	  
	  // assemble the DOFMatrix block and the corresponding rhs vector, of r==c
	  this->assemble({int(_r), int(_c)}, &dofmatrix, (_r == _c && asmVector ? &rhs_vec : NULL));

	  // finish boundary condition
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
              for (auto bc : bc_list.second)
                bc->finish(c == int(_c), dofmatrix, solution_vec, rhs_vec);
	    }
	  }
	  
	  nnz += dofmatrix.getNnz();
	}
      });
    });
    
    AMDIS_MSG("fillin of assembled matrix: " << nnz);
    
    AMDIS_MSG("buildAfterCoarsen needed " << t.elapsed() << " seconds");
  }
  

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


  template <class Traits>
    template <class Matrix, class Vector>
  void ProblemStatSeq<Traits>::assemble(std::pair<int, int> row_col,
                                        Matrix* dofmatrix,
                                        Vector* rhs)
  {
    auto const& rowFeSpace = dofmatrix->getRowFeSpace();
    auto const& colFeSpace = dofmatrix->getColFeSpace();
    
    dofmatrix->init();
    auto matrixOp = matrixOperators[row_col];
    auto matrixOpFac = matrixFactors[row_col];
    auto vectorOp = vectorOperators[row_col.first];
    auto vectorOpFac = vectorFactors[row_col.first];
    
    for (auto op : matrixOp)
      op->init(rowFeSpace, colFeSpace);
    for (auto op : vectorOp)
      op->init(rowFeSpace, colFeSpace);
            
    auto rowLocalView = rowFeSpace.localView();
    auto rowIndexSet = rowFeSpace.localIndexSet();
    
    auto colLocalView = colFeSpace.localView();
    auto colIndexSet = colFeSpace.localIndexSet();
    
    for (const auto& element : elements(*meshView)) {
      rowLocalView.bind(element);
      rowIndexSet.bind(rowLocalView);
      
      colLocalView.bind(element);
      colIndexSet.bind(colLocalView);
      
      if (dofmatrix) {
        ElementMatrix elementMatrix;
        bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix, 
                                    matrixOp, matrixOpFac);
        if (add)
          addElementMatrix(*dofmatrix, rowIndexSet, colIndexSet, elementMatrix);
      }
      
      if (rhs) {
        ElementVector elementVector;
        bool add = getElementVector(rowLocalView, elementVector, 
                                    vectorOp, vectorOpFac);
        if (add)
          addElementVector(*rhs, rowIndexSet, elementVector);
      }
    }
    dofmatrix->finish();
  }
  
  
  template <class Traits>
    template <class RowView, class ColView>
  bool ProblemStatSeq<Traits>::getElementMatrix(RowView const& rowLocalView,
                                                ColView const& colLocalView,
                                                ElementMatrix& elementMatrix,
                                                std::list<std::shared_ptr<OperatorType>>& operators,
                                                std::list<double*> const& factors)
  {
    const auto nRows = rowLocalView.tree().finiteElement().size();
    const auto nCols = colLocalView.tree().finiteElement().size();

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

    // fills the entire matrix with zeroes
    elementMatrix.setSize(nRows, nCols);
    elementMatrix = 0.0;

    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<std::shared_ptr<OperatorType>>& operators,
                                                std::list<double*> const& factors)
  {
    const auto nRows = rowLocalView.tree().finiteElement().size();

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

    // Set all entries to zero
    elementVector.resize(nRows);
    elementVector = 0.0;
    
    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 < elementMatrix.N(); ++i) {
      // The global index of the i−th vertex of the element
      const auto row = rowIndexSet.index(i);
      for (size_t j = 0; j < elementMatrix.M(); ++j) {
        // The global index of the j−th vertex of the element
        const auto 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 < elementVector.size(); ++i) {
      // The global index of the i-th vertex
      const auto row = indexSet.index(i);
      dofvector[row] += elementVector[i];
    }
  }

} // end namespace AMDiS