directlinearsolver.h

Go to the documentation of this file.

#ifndef DIRECT_LINEAR_SOLVER_H_
#define DIRECT_LINEAR_SOLVER_H_

#include <deal.II/lac/vector.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/numerics/vector_tools.h>

#include <vector>

#include "parameterreader.h"

namespace DOpE
{
  template <typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
    class DirectLinearSolverWithMatrix
  {
  public:
    DirectLinearSolverWithMatrix(ParameterReader &param_reader);
    ~DirectLinearSolverWithMatrix();
   
    static void declare_params(ParameterReader &param_reader);
 
    template<typename PROBLEM>
      void ReInit(PROBLEM& pde);

    template<typename PROBLEM, typename INTEGRATOR>
      void Solve(PROBLEM& pde, INTEGRATOR& integr, VECTOR &rhs, VECTOR &solution, bool force_matrix_build=false);
     
  protected:

  private:
    SPARSITYPATTERN sparsity_pattern_;
    MATRIX matrix_;

    dealii::SparseDirectUMFPACK* A_direct_;

  };

/*********************************Implementation************************************************/
 
  template <typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
    void DirectLinearSolverWithMatrix<SPARSITYPATTERN,MATRIX,VECTOR>::declare_params(ParameterReader &/*param_reader*/) 
  {
  }

  /******************************************************/

  template <typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
    DirectLinearSolverWithMatrix<SPARSITYPATTERN,MATRIX,VECTOR>::DirectLinearSolverWithMatrix(
      ParameterReader &/*param_reader*/) 
  {
    A_direct_ = NULL;
  }

/******************************************************/

template <typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
 DirectLinearSolverWithMatrix<SPARSITYPATTERN,MATRIX,VECTOR>::~DirectLinearSolverWithMatrix()
{
  if(A_direct_ != NULL)
  {
    delete A_direct_;
  }
}

/******************************************************/

template <typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
  template<typename PROBLEM>
  void  DirectLinearSolverWithMatrix<SPARSITYPATTERN,MATRIX,VECTOR>::ReInit(PROBLEM& pde)
{
  matrix_.clear();
  pde.ComputeSparsityPattern(sparsity_pattern_);
  matrix_.reinit(sparsity_pattern_);
  
  if(A_direct_ != NULL)
  {
    delete A_direct_;
    A_direct_= NULL;
  }
}

/******************************************************/

template <typename SPARSITYPATTERN, typename MATRIX, typename VECTOR>
  template<typename PROBLEM, typename INTEGRATOR>
  void DirectLinearSolverWithMatrix<SPARSITYPATTERN,MATRIX,VECTOR>::Solve(PROBLEM& pde,
                                                                          INTEGRATOR& integr,
                                                                          VECTOR &rhs, 
                                                                          VECTOR &solution, 
                                                                          bool force_matrix_build)
{
  if(force_matrix_build)
  { 
    integr.ComputeMatrix (pde,matrix_);    
  }

  integr.ApplyNewtonBoundaryValues(pde,matrix_,rhs,solution);
  
  if(A_direct_ == NULL)
  {
    A_direct_ = new dealii::SparseDirectUMFPACK;
    A_direct_->initialize(matrix_);
  }
  else if(force_matrix_build)
  {
     A_direct_->factorize(matrix_);
  }
 
  dealii::Vector<double> sol;
  sol = rhs;
  A_direct_->solve(sol);
  solution = sol;

  pde.GetDoFConstraints().distribute(solution);

}


}
#endif