amandus.h

Go to the documentation of this file.

/**********************************************************************
 *  Copyright (C) 2011 - 2014 by the authors
 *  Distributed under the MIT License
 *
 * See the files AUTHORS and LICENSE in the project root directory
 *
 **********************************************************************/

#ifndef amandus_amandus_h
#define amandus_amandus_h

#include <deal.II/lac/arpack_solver.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/precondition_block.h>
#include <deal.II/lac/relaxation_block.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_richardson.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparse_matrix.h>

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_refinement.h>

#include <deal.II/dofs/dof_tools.h>

#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/integration_info.h>
#include <deal.II/meshworker/loop.h>
#include <deal.II/meshworker/output.h>
#include <deal.II/meshworker/simple.h>

#include <deal.II/multigrid/mg_coarse.h>
#include <deal.II/multigrid/mg_matrix.h>
#include <deal.II/multigrid/mg_smoother.h>
#include <deal.II/multigrid/mg_tools.h>
#include <deal.II/multigrid/mg_transfer.h>
#include <deal.II/multigrid/multigrid.h>

#include <deal.II/algorithms/operator.h>
#include <deal.II/base/flow_function.h>
#include <deal.II/base/function_lib.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>

#include <amandus/integrator.h>

#include <fstream>
#include <iostream>

DeclExceptionMsg(ExcNeedArpack,
                 "To use this functionality, deal.II must be configured with Arpack.");

class AmandusParameters : public dealii::ParameterHandler
{
public:
  AmandusParameters();

  void read(int argc, const char** argv);
};

template <int dim>
class AmandusApplicationSparse : public dealii::Subscriptor
{
public:
  typedef dealii::MeshWorker::IntegrationInfo<dim> CellInfo;

  AmandusApplicationSparse(dealii::Triangulation<dim>& triangulation,
                           const dealii::FiniteElement<dim>& fe, bool use_umfpack = false);

  virtual void parse_parameters(dealii::ParameterHandler& param);

  void set_number_of_matrices(unsigned int n);

  void set_boundary(dealii::types::boundary_id index,
                    dealii::ComponentMask mask = dealii::ComponentMask());

  void set_meanvalue(dealii::ComponentMask mask = dealii::ComponentMask());

  virtual void setup_vector(dealii::Vector<double>& v) const;

  virtual void update_vector_inhom_boundary(dealii::Vector<double>& v,
                                            const dealii::Function<dim>& inhom_boundary,
                                            bool projection = false) const;

  virtual void setup_system();

  void assemble_right_hand_side(dealii::AnyData& out, const dealii::AnyData& in,
                                const AmandusIntegrator<dim>& integrator) const;

  void refine_mesh(const bool global = false);

  const dealii::DoFHandler<dim>& dofs() const;

  const dealii::ConstraintMatrix& constraints() const;

  const dealii::ConstraintMatrix& hanging_nodes() const;

  virtual void setup_constraints();
  void assemble_matrix(const dealii::AnyData& in, const AmandusIntegrator<dim>& integrator);
  virtual void assemble_mg_matrix(const dealii::AnyData& in,
                                  const AmandusIntegrator<dim>& integrator);

  const dealii::Vector<double>& indicators() const;

  double estimate(const dealii::AnyData& in, AmandusIntegrator<dim>& integrator);
  void error(dealii::BlockVector<double>& out, const dealii::AnyData& in,
             const AmandusIntegrator<dim>& integrator);

  void error(dealii::BlockVector<double>& out, const dealii::AnyData& in,
             const ErrorIntegrator<dim>& integrator);

  void error(const dealii::AnyData& in, const AmandusIntegrator<dim>& integrator,
             unsigned int num_errs);
  virtual void solve(dealii::Vector<double>& sol, const dealii::Vector<double>& rhs);

  virtual void arpack_solve(std::vector<std::complex<double>>& eigenvalues,
                            std::vector<dealii::Vector<double>>& eigenvectors);

  void output_results(unsigned int refinement_cycle, const dealii::AnyData* data = 0) const;

  void verify_residual(dealii::AnyData& out, const dealii::AnyData& in,
                       const AmandusIntegrator<dim>& integrator) const;

  dealii::ReductionControl control;

  dealii::SmartPointer<dealii::ParameterHandler> param;

  typename dealii::Triangulation<dim>::Signals& signals;

protected:
  dealii::SmartPointer<dealii::Triangulation<dim>, AmandusApplicationSparse<dim>> triangulation;

  const dealii::MappingQ1<dim> mapping;

  dealii::SmartPointer<const dealii::FiniteElement<dim>, AmandusApplicationSparse<dim>> fe;

  dealii::DoFHandler<dim> dof_handler;

  std::map<dealii::types::boundary_id, dealii::ComponentMask> boundary_masks;

  dealii::ComponentMask meanvalue_mask;

  dealii::ConstraintMatrix constraint_matrix;

  dealii::ConstraintMatrix hanging_node_constraints;

  dealii::SparsityPattern sparsity;
  std::vector<dealii::SparseMatrix<double>> matrix;
  const bool use_umfpack;
  dealii::SparseDirectUMFPACK inverse;

  dealii::BlockVector<double> estimates;

  std::vector<dealii::DataComponentInterpretation::DataComponentInterpretation> output_data_types;
};

template <int dim, typename RELAXATION = dealii::RelaxationBlockSSOR<dealii::SparseMatrix<double>>>
class AmandusApplication : public AmandusApplicationSparse<dim>
{
public:
  typedef dealii::MeshWorker::IntegrationInfo<dim> CellInfo;

  AmandusApplication(dealii::Triangulation<dim>& triangulation,
                     const dealii::FiniteElement<dim>& fe);

  virtual void parse_parameters(dealii::ParameterHandler& param);

  void set_advection_directions(std::vector<dealii::Tensor<1, dim>>& adv_dir);

  void setup_system();

  void setup_constraints();

  void assemble_mg_matrix(const dealii::AnyData& in, const AmandusIntegrator<dim>& integrator);

  void solve(dealii::Vector<double>& sol, const dealii::Vector<double>& rhs);

  virtual void arpack_solve(std::vector<std::complex<double>>& eigenvalues,
                            std::vector<dealii::Vector<double>>& eigenvectors);

  bool vertex_patches = true;
  bool boundary_patches = false;

protected:
  dealii::MGConstrainedDoFs mg_constraints;

  dealii::MGLevelObject<dealii::SparsityPattern> mg_sparsity;
  dealii::MGLevelObject<dealii::SparsityPattern> mg_sparsity_fluxes;
  dealii::MGLevelObject<dealii::SparseMatrix<double>> mg_matrix;

  dealii::MGLevelObject<dealii::SparseMatrix<double>> mg_matrix_down;
  dealii::MGLevelObject<dealii::SparseMatrix<double>> mg_matrix_up;
  dealii::MGLevelObject<dealii::SparseMatrix<double>> mg_matrix_flux_down;
  dealii::MGLevelObject<dealii::SparseMatrix<double>> mg_matrix_flux_up;

  dealii::MGTransferPrebuilt<dealii::Vector<double>> mg_transfer;

  dealii::FullMatrix<double> coarse_matrix;
  dealii::MGCoarseGridSVD<double, dealii::Vector<double>> mg_coarse;

  dealii::MGLevelObject<typename RELAXATION::AdditionalData> smoother_data;
  dealii::mg::SmootherRelaxation<RELAXATION, dealii::Vector<double>> mg_smoother;

  std::vector<dealii::Tensor<1, dim>> advection_directions;

  bool log_smoother_statistics = false;
  bool right_preconditioning = true;
  bool use_default_residual = true;
  double smoother_relaxation = 1.0;
};

template <int dim>
class AmandusUMFPACK : public AmandusApplicationSparse<dim>
{
public:
  AmandusUMFPACK(dealii::Triangulation<dim>& triangulation, const dealii::FiniteElement<dim>& fe);
};

template <int dim>
class AmandusResidual : public dealii::Algorithms::OperatorBase
{
public:
  AmandusResidual(const AmandusApplicationSparse<dim>& application,
                  AmandusIntegrator<dim>& integrator);

  virtual void operator()(dealii::AnyData& out, const dealii::AnyData& in);

protected:
  dealii::SmartPointer<const AmandusApplicationSparse<dim>, AmandusResidual<dim>> application;
  dealii::SmartPointer<AmandusIntegrator<dim>, AmandusResidual<dim>> integrator;
};

template <int dim>
class AmandusSolve : public dealii::Algorithms::OperatorBase
{
public:
  AmandusSolve(AmandusApplicationSparse<dim>& application, AmandusIntegrator<dim>& integrator);
  virtual void operator()(dealii::AnyData& out, const dealii::AnyData& in);

private:
  dealii::SmartPointer<AmandusApplicationSparse<dim>, AmandusSolve<dim>> application;
  dealii::SmartPointer<AmandusIntegrator<dim>, AmandusSolve<dim>> integrator;
};

template <int dim>
inline void
AmandusApplicationSparse<dim>::set_number_of_matrices(unsigned int n)
{
  matrix.resize(n);
}

template <int dim>
inline const dealii::DoFHandler<dim>&
AmandusApplicationSparse<dim>::dofs() const
{
  return dof_handler;
}

template <int dim>
inline const dealii::ConstraintMatrix&
AmandusApplicationSparse<dim>::constraints() const
{
  return constraint_matrix;
}

template <int dim>
inline const dealii::ConstraintMatrix&
AmandusApplicationSparse<dim>::hanging_nodes() const
{
  return hanging_node_constraints;
}

template <int dim>
inline const dealii::Vector<double>&
AmandusApplicationSparse<dim>::indicators() const
{
  return estimates.block(0);
}

template <int dim, typename RELAXATION>
inline void AmandusApplication<dim, RELAXATION>::set_advection_directions(
  std::vector<dealii::Tensor<1, dim>>& adv_dir)
{
  advection_directions = adv_dir;
}

#endif