project_Hgrad_to_Hdiv.hpp

/******************************************************************************
 * Copyright (C) Siarhei Uzunbajakau, 2023.
 *
 * This program is free software. You can use, modify, and redistribute it under
 * the terms of the GNU Lesser General Public License as published by the Free
 * Software Foundation, either version 3 or (at your option) any later version.
 * This program is distributed without any warranty.
 *
 * Refer to COPYING.LESSER for more details.
 ******************************************************************************/
 
#ifndef ProjectHgradToHdiv_H__
#define ProjectHgradToHdiv_H__
 
#include <deal.II/base/timer.h>
#include <deal.II/base/work_stream.h>
 
#include <deal.II/grid/tria.h>
 
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
 
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_direct.h>
 
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_control.h>
 
#include <deal.II/fe/fe_raviart_thomas.h>
 
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
 
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
 
#include <fstream>
#include <iomanip>
#include <ios>
#include <iostream>
#include <string>
 
#include "constants.hpp"
#include "static_scalar_input.hpp"
 
#define VE scratch_data.ve
 
#define TMR(__name) TimerOutput::Scope timer_section(timer, __name)
 
using namespace dealii;
 
namespace StaticScalarSolver {
template<int dim, int stage = 1>
class ProjectHgradToHdiv
{
public:
  ProjectHgradToHdiv() = delete;
 
  ProjectHgradToHdiv(unsigned int p,
                     unsigned int mapping_degree,
                     const Triangulation<dim>& triangulation_Hgrad,
                     const DoFHandler<dim>& dof_handler_Hgrad,
                     const Vector<double>& solution_Hgrad,
                     std::string fname = "Hdiv",
                     const Function<dim>* exact_solution = nullptr,
                     bool axisymmetric = false,
                     bool vector_potential = false,
                     bool print_time_tables = false,
                     bool project_exact_solution = false,
                     bool log_cg_convergence = false,
                     bool write_higher_order_cells = false);
 
  double get_L2_norm() { return L2_norm; };
 
  double get_Linfty_norm() { return Linfty_norm; }
 
  unsigned int get_n_cells() const
  {
    return static_cast<unsigned int>(triangulation_Hgrad.n_active_cells());
  }
 
  unsigned int get_n_dofs() const
  {
    return static_cast<unsigned int>(dof_handler_Hdiv.n_dofs());
  }
 
  void clear()
  {
    system_matrix.clear();
    system_rhs.reinit(0);
  }
 
  const Triangulation<dim>& get_tria() const { return triangulation_Hgrad; }
 
  const DoFHandler<dim>& get_dof_handler() const { return dof_handler_Hdiv; }
 
  const Vector<double>& get_solution() const { return solution_Hdiv; }
 
  void save_matrix_and_rhs_to_csv(std::string fname) const;
 
private:
  void setup();
  void assemble();
  void solve();
  void save() const;
  void compute_error_norms();
  void project_exact_solution_fcn();
 
  const std::string fname;
 
  const DoFHandler<dim>& dof_handler_Hgrad;
  const Vector<double>& solution_Hgrad;
 
  const Triangulation<dim>& triangulation_Hgrad;
  const FE_RaviartThomas<dim> fe_Hdiv;
  DoFHandler<dim> dof_handler_Hdiv;
 
  SparsityPattern sparsity_pattern;
  SparseMatrix<double> system_matrix;
 
  Vector<double> solution_Hdiv;
  Vector<double> system_rhs;
 
  Vector<double> projected_exact_solution;
 
  AffineConstraints<double> constraints;
 
  const Function<dim>* exact_solution;
 
  const unsigned int mapping_degree;
  const bool axisymmetric;
  const bool vector_potential;
  const bool project_exact_solution;
  const bool log_cg_convergence;
  const bool write_higher_order_cells;
 
  Vector<double> L2_per_cell;
  double L2_norm;
 
  Vector<double> Linfty_per_cell;
  double Linfty_norm;
 
  // ----------------------------------------------------------------------------
  // These structures and functions are related to the Work Stream algorithm.
  // See article "WorkStream – A Design Pattern for Multicore-Enabled Finite
  // Element Computations." by BRUNO TURCKSIN, MARTIN KRONBICHLER,
  // WOLFGANG BANGERTH for more details.
  // ----------------------------------------------------------------------------
 
  using IteratorTuple =
    std::tuple<typename DoFHandler<dim>::active_cell_iterator,
               typename DoFHandler<dim>::active_cell_iterator>;
 
  using IteratorPair = SynchronousIterators<IteratorTuple>;
 
  struct AssemblyScratchData
  {
    AssemblyScratchData(const FiniteElement<dim>& fe,
                        const DoFHandler<dim>& dof_handr_Hgrad,
                        const Vector<double>& dofs_Hgrad,
                        bool axisymmetric,
                        bool vector_potential,
                        unsigned int mapping_degree);
 
    AssemblyScratchData(const AssemblyScratchData& scratch_data);
 
    MappingQ<dim> mapping;
    Constants::QuadratureTableVector<dim> qt;
    FEValues<dim> fe_values_Hdiv;
    FEValues<dim> fe_values_Hgrad;
 
    const unsigned int dofs_per_cell;
    const unsigned int n_q_points;
 
    TheCoefficient<dim, stage> the_coefficient;
    std::vector<double> the_coefficient_list;
 
    std::vector<Tensor<1, dim>> vector_gradients;
 
    // Two-dimensional vector curl of an out-of-plane (oop) vector.
    std::vector<Tensor<1, dim>> nabla_xV_oopvector;
 
    const FEValuesExtractors::Vector ve;
 
    const DoFHandler<dim>& dof_hand_Hgrad;
    const Vector<double>& dofs_Hgrad;
 
    const bool axisymmetric;
    const bool vector_potential;
 
    double axi_mult; // Equals the distance to the axis of rotation symmetry, r,
                     // in the recipes for axisymmetric projections. Equals 1.0
                     // in all other recipes. All integrands are multiplied by
                     // this multiplier.
  };
 
  struct AssemblyCopyData
  {
    FullMatrix<double> cell_matrix;
    Vector<double> cell_rhs;
    std::vector<types::global_dof_index> local_dof_indices;
  };
 
  void system_matrix_local(const IteratorPair& IP,
                           AssemblyScratchData& scratch_data,
                           AssemblyCopyData& copy_data);
 
  void copy_local_to_global(const AssemblyCopyData& copy_data);
 
  //-----------------------------------------------------------------------------
  //-----------------------------------------------------------------------------
  //-----------------------------------------------------------------------------
};
 
template<int dim, int stage>
ProjectHgradToHdiv<dim, stage>::ProjectHgradToHdiv(
  unsigned int p,
  unsigned int mapping_degree,
  const Triangulation<dim>& triangulation_Hgrad,
  const DoFHandler<dim>& dof_handler_Hgrad,
  const Vector<double>& solution_Hgrad,
  std::string fname,
  const Function<dim>* exact_solution,
  bool axisymmetric,
  bool vector_potential,
  bool print_time_tables,
  bool project_exact_solution,
  bool log_cg_convergence,
  bool write_higher_order_cells)
  : fname(fname)
  , dof_handler_Hgrad(dof_handler_Hgrad)
  , solution_Hgrad(solution_Hgrad)
  , triangulation_Hgrad(triangulation_Hgrad)
  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  // The public attribute fe.degree is the maximal polynomial degree of a
  // shape function in a single coordinate direction, not the degree of
  // the finite element. For FE_Nedelec and FE_RaviartThomas degree of
  // the finite element is: degree_of_element = fe.degree - 1.
  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  //    fe_Hdiv(dof_handler_Hgrad.get_fe().degree - 1 ),
  , fe_Hdiv(p)
  , exact_solution(exact_solution)
  , mapping_degree(mapping_degree)
  , axisymmetric(axisymmetric)
  , vector_potential(vector_potential)
  , project_exact_solution(project_exact_solution)
  , log_cg_convergence(log_cg_convergence)
  , write_higher_order_cells(write_higher_order_cells)
{
  if (axisymmetric) {
    Assert(
      dim == 2,
      ExcMessage("The setting axisymmetric=true is only allowed if dim=2."));
  }
 
  if (vector_potential) {
    Assert(dim == 2,
           ExcMessage(
             "The setting vector_potential=true can only be used if dim=2."));
  }
 
  TimerOutput::OutputFrequency tf =
    (print_time_tables) ? TimerOutput::summary : TimerOutput::never;
 
  TimerOutput timer(std::cout, tf, TimerOutput::cpu_and_wall_times_grouped);
 
  {
    TMR("Setup");
    setup();
  }
  {
    TMR("Assemble");
    assemble();
  }
  {
    TMR("Solve");
    solve();
  }
 
  if (exact_solution) {
    {
      TMR("Compute error norms");
      compute_error_norms();
    }
 
    if (project_exact_solution) {
      {
        TMR("Project exact solution");
        project_exact_solution_fcn();
      }
    }
  }
 
  {
    TMR("Save");
    save();
  }
}
 
template<int dim, int stage>
void
ProjectHgradToHdiv<dim, stage>::save_matrix_and_rhs_to_csv(
  std::string fname) const
{
  std::ofstream ofs_matrix(fname + "_matrix.csv");
  std::ofstream ofs_rhs(fname + "_rhs.csv");
 
  for (unsigned int i = 0; i < system_matrix.m(); ++i) {
    ofs_rhs << system_rhs(i);
    if (i < (system_matrix.m() - 1))
      ofs_rhs << "\n";
 
    for (unsigned int j = 0; j < system_matrix.n(); ++j) {
      ofs_matrix << std::scientific << std::setprecision(16)
                 << system_matrix.el(i, j);
 
      if (j < (system_matrix.m() - 1))
        ofs_matrix << ", ";
    }
    if (i < (system_matrix.m() - 1))
      ofs_matrix << "\n";
  }
 
  ofs_rhs.close();
  ofs_matrix.close();
}
 
template<int dim, int stage>
void
ProjectHgradToHdiv<dim, stage>::setup()
{
  constraints.close();
 
  dof_handler_Hdiv.reinit(triangulation_Hgrad);
  dof_handler_Hdiv.distribute_dofs(fe_Hdiv);
 
  DynamicSparsityPattern dsp(dof_handler_Hdiv.n_dofs(),
                             dof_handler_Hdiv.n_dofs());
  DoFTools::make_sparsity_pattern(dof_handler_Hdiv, dsp, constraints, false);
 
  sparsity_pattern.copy_from(dsp);
  system_matrix.reinit(sparsity_pattern);
  solution_Hdiv.reinit(dof_handler_Hdiv.n_dofs());
  system_rhs.reinit(dof_handler_Hdiv.n_dofs());
 
  if (project_exact_solution)
    projected_exact_solution.reinit(dof_handler_Hdiv.n_dofs());
 
  if (exact_solution) {
    L2_per_cell.reinit(triangulation_Hgrad.n_active_cells());
    Linfty_per_cell.reinit(triangulation_Hgrad.n_active_cells());
  }
}
 
template<int dim, int stage>
void
ProjectHgradToHdiv<dim, stage>::assemble()
{
  WorkStream::run(IteratorPair(IteratorTuple(dof_handler_Hdiv.begin_active(),
                                             dof_handler_Hgrad.begin_active())),
                  IteratorPair(IteratorTuple(dof_handler_Hdiv.end(),
                                             dof_handler_Hgrad.end())),
                  *this,
                  &ProjectHgradToHdiv::system_matrix_local,
                  &ProjectHgradToHdiv::copy_local_to_global,
                  AssemblyScratchData(fe_Hdiv,
                                      dof_handler_Hgrad,
                                      solution_Hgrad,
                                      axisymmetric,
                                      vector_potential,
                                      mapping_degree),
                  AssemblyCopyData());
}
 
template<int dim, int stage>
ProjectHgradToHdiv<dim, stage>::AssemblyScratchData::AssemblyScratchData(
  const FiniteElement<dim>& fe,
  const DoFHandler<dim>& dof_hand_Hgrad,
  const Vector<double>& dofs_Hgrad,
  bool axisymmetric,
  bool vector_potential,
  unsigned int mapping_degree)
  : mapping(mapping_degree)
  , qt(fe.degree - 1)
  , fe_values_Hdiv(mapping,
                   fe,
                   QGauss<dim>(qt.sim()),
                   update_values | update_quadrature_points | update_JxW_values)
  , fe_values_Hgrad(mapping,
                    dof_hand_Hgrad.get_fe(),
                    QGauss<dim>(qt.sim()),
                    update_gradients)
  , dofs_per_cell(fe_values_Hdiv.dofs_per_cell)
  , n_q_points(fe_values_Hdiv.get_quadrature().size())
  , the_coefficient_list(n_q_points)
  , vector_gradients(n_q_points, Tensor<1, dim>())
  , nabla_xV_oopvector(n_q_points, Tensor<1, dim>())
  , ve(0)
  , dof_hand_Hgrad(dof_hand_Hgrad)
  , dofs_Hgrad(dofs_Hgrad)
  , axisymmetric(axisymmetric)
  , vector_potential(vector_potential)
  , axi_mult(1.0)
{
}
 
template<int dim, int stage>
ProjectHgradToHdiv<dim, stage>::AssemblyScratchData::AssemblyScratchData(
  const AssemblyScratchData& scratch_data)
  : mapping(scratch_data.mapping.get_degree())
  , qt(scratch_data.qt)
  , fe_values_Hdiv(mapping,
                   scratch_data.fe_values_Hdiv.get_fe(),
                   scratch_data.fe_values_Hdiv.get_quadrature(),
                   update_values | update_quadrature_points | update_JxW_values)
  , fe_values_Hgrad(mapping,
                    scratch_data.fe_values_Hgrad.get_fe(),
                    scratch_data.fe_values_Hgrad.get_quadrature(),
                    update_gradients)
  , dofs_per_cell(fe_values_Hdiv.dofs_per_cell)
  , n_q_points(fe_values_Hdiv.get_quadrature().size())
  , the_coefficient_list(n_q_points)
  , vector_gradients(n_q_points, Tensor<1, dim>())
  , nabla_xV_oopvector(n_q_points, Tensor<1, dim>())
  , ve(0)
  , dof_hand_Hgrad(scratch_data.dof_hand_Hgrad)
  , dofs_Hgrad(scratch_data.dofs_Hgrad)
  , axisymmetric(scratch_data.axisymmetric)
  , vector_potential(scratch_data.vector_potential)
  , axi_mult(1.0)
{
}
 
template<int dim, int stage>
void
ProjectHgradToHdiv<dim, stage>::system_matrix_local(
  const IteratorPair& IP,
  AssemblyScratchData& scratch_data,
  AssemblyCopyData& copy_data)
{
  // See the color boxes
  // (1) Recipe for projections from H(grad) to H(div) nr. 6 and 7
  // (2) Recipe for projections from H(grad) to H(div) nr. 8 and 9 (planar)
  // (3) Recipe for projections from H(grad) to H(div) nr. 8 and 9 (axisym.)
  // (4) Recipe for projections from H(grad) to H(div) nr. 10 and 11 (planar)
  // (5) Recipe for projections from H(grad) to H(div) nr. 10 (axisym.)
  //
  // The comments below refer to these recipes by number, i.e., recipe (1),
  // recipe (2), etc.
 
  copy_data.cell_matrix.reinit(scratch_data.dofs_per_cell,
                               scratch_data.dofs_per_cell);
 
  copy_data.cell_rhs.reinit(scratch_data.dofs_per_cell);
 
  copy_data.local_dof_indices.resize(scratch_data.dofs_per_cell);
 
  scratch_data.fe_values_Hdiv.reinit(std::get<0>(*IP));
  scratch_data.fe_values_Hgrad.reinit(std::get<1>(*IP));
 
  scratch_data.fe_values_Hgrad.get_function_gradients(
    scratch_data.dofs_Hgrad, scratch_data.vector_gradients);
 
  if (!scratch_data.vector_potential)
    scratch_data.the_coefficient.value_list(
      scratch_data.fe_values_Hdiv.get_quadrature_points(),
      std::get<0>(*IP)->material_id(),
      std::get<0>(*IP)->user_index(),
      scratch_data.the_coefficient_list);
 
  for (unsigned int q_index = 0; q_index < scratch_data.n_q_points; ++q_index) {
 
    scratch_data.axi_mult = 1.0;
 
    if ((scratch_data.axisymmetric) && (!scratch_data.vector_potential))
      scratch_data.axi_mult =
        scratch_data.fe_values_Hdiv.quadrature_point(q_index)[0];
 
    for (unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i) {
      for (unsigned int j = 0; j < scratch_data.dofs_per_cell; ++j) {
        // Mass matrix is the same in all recipes.
        copy_data.cell_matrix(i, j) +=
          scratch_data.axi_mult *                             // 1.0 or r
          scratch_data.fe_values_Hdiv[VE].value(i, q_index) * // curl N_i
          scratch_data.fe_values_Hdiv[VE].value(j, q_index) * // curl N_j
          scratch_data.fe_values_Hdiv.JxW(q_index);           // dV (dS in 2D)
      }
 
      if (scratch_data.vector_potential) // A and A'.
      {
        // The option vector_potential = true exists only in 2D
        scratch_data.nabla_xV_oopvector[q_index][0] =
          scratch_data.vector_gradients[q_index][1];
 
        scratch_data
          .nabla_xV_oopvector[q_index][1] = // nabla_xV_oopvector is the
          -scratch_data.vector_gradients[q_index][0]; // two-dimensional vector
                                                      // field in H(div),
                                                      // 'curl_v A' for short.
        double tmp;
        tmp = scratch_data.nabla_xV_oopvector[q_index] *          // curl_v A
              scratch_data.fe_values_Hdiv[VE].value(i, q_index) * // N_i
              scratch_data.fe_values_Hdiv.JxW(q_index);           // dS
 
        if (scratch_data.axisymmetric) // A'.
        {                              // Integral b_i in recipe (5).
          copy_data.cell_rhs(i) -= tmp;
 
          // The recipe (5) yields B'=rB, where r is the distance to the axis of
          // rotation symmetry. If you do not like to have the scaled magnetic
          // field, B', and would rather have the magnetic field itself, B,
          // replace the line above with the following line. Note, that a
          // quadrature point can, in general, be at the origin. If so, the line
          // below will yield an error.
 
          //          copy_data.cell_rhs(i) -=
          // tmp/scratch_data.fe_values_Hdiv.quadrature_point(q_index)[0];
        } else // A.
        {      // Integral b_i in recipe (4).
          copy_data.cell_rhs(i) += tmp;
        }
 
      } else // Phi, Psi, Theta.
      {
        // Integral b_i in recipes (1), (2), and (3).
        copy_data.cell_rhs(i) -=
          scratch_data.axi_mult *                             // 1.0 or r
          scratch_data.the_coefficient_list[q_index] *        // epsilon
          scratch_data.vector_gradients[q_index] *            // grad PHI
          scratch_data.fe_values_Hdiv[VE].value(i, q_index) * // N_i
          scratch_data.fe_values_Hdiv.JxW(q_index);           // dV (dS in 2D)
      }
    }
  }
 
  std::get<0>(*IP)->get_dof_indices(copy_data.local_dof_indices);
}
 
template<int dim, int stage>
void
ProjectHgradToHdiv<dim, stage>::copy_local_to_global(
  const AssemblyCopyData& copy_data)
{
  constraints.distribute_local_to_global(copy_data.cell_matrix,
                                         copy_data.cell_rhs,
                                         copy_data.local_dof_indices,
                                         system_matrix,
                                         system_rhs);
}
 
template<int dim, int stage>
void
ProjectHgradToHdiv<dim, stage>::solve()
{
  SolverControl control(
    1000 * system_rhs.size(), 1e-12 * system_rhs.l2_norm(), false, false);
 
  if (log_cg_convergence)
    control.enable_history_data();
 
  GrowingVectorMemory<Vector<double>> memory;
  SolverCG<Vector<double>> cg(control, memory);
 
  PreconditionJacobi<SparseMatrix<double>> preconditioner;
  preconditioner.initialize(system_matrix, 1.0);
 
  cg.solve(system_matrix, solution_Hdiv, system_rhs, preconditioner);
 
  if (log_cg_convergence) {
    const std::vector<double> history_data = control.get_history_data();
 
    std::ofstream ofs(fname + "_cg_convergence.csv");
 
    unsigned int i = 1;
    for (auto item : history_data) {
      ofs << i << ", " << item << "\n";
      i++;
    }
 
    ofs.close();
  }
}
 
template<int dim, int stage>
void
ProjectHgradToHdiv<dim, stage>::save() const
{
  std::vector<std::string> solution_names(dim, "VectorField");
  std::vector<DataComponentInterpretation::DataComponentInterpretation>
    interpretation(dim,
                   DataComponentInterpretation::component_is_part_of_vector);
 
  DataOut<dim> data_out;
 
  data_out.add_data_vector(
    dof_handler_Hdiv, solution_Hdiv, solution_names, interpretation);
 
  if (project_exact_solution) {
    std::vector<std::string> solution_names_ex(dim, "VectorFieldExact");
 
    data_out.add_data_vector(dof_handler_Hdiv,
                             projected_exact_solution,
                             solution_names_ex,
                             interpretation);
  }
 
  data_out.add_data_vector(L2_per_cell, "L2norm");
  data_out.add_data_vector(Linfty_per_cell, "LinftyNorm");
 
  std::ofstream ofs;
 
  if (write_higher_order_cells) {
    DataOutBase::VtkFlags flags;
    flags.write_higher_order_cells = true;
    data_out.set_flags(flags);
 
    const MappingQ<dim> mapping(mapping_degree);
 
    data_out.build_patches(mapping,
                           fe_Hdiv.degree + 2,
                           DataOut<dim>::CurvedCellRegion::curved_inner_cells);
 
    ofs.open(fname + ".vtu");
    data_out.write_vtu(ofs);
 
  } else {
 
    data_out.build_patches();
 
    ofs.open(fname + ".vtk");
    data_out.write_vtk(ofs);
  }
 
  ofs.close();
}
 
template<int dim, int stage>
void
ProjectHgradToHdiv<dim, stage>::compute_error_norms()
{
  Weight<dim, stage> weight;
  const Function<dim, double>* mask = &weight;
 
  Constants::QuadratureTableVector<dim> qt(dof_handler_Hdiv.get_fe().degree -
                                           1);
  QGauss<dim> quadrature(qt.enorm());
 
  VectorTools::integrate_difference(MappingQ<dim>(mapping_degree),
                                    dof_handler_Hdiv,
                                    solution_Hdiv,
                                    *exact_solution,
                                    L2_per_cell,
                                    quadrature,
                                    VectorTools::L2_norm,
                                    mask);
 
  L2_norm = VectorTools::compute_global_error(
    triangulation_Hgrad, L2_per_cell, VectorTools::L2_norm);
 
  VectorTools::integrate_difference(MappingQ<dim>(mapping_degree),
                                    dof_handler_Hdiv,
                                    solution_Hdiv,
                                    *exact_solution,
                                    Linfty_per_cell,
                                    QGauss<dim>(1),
                                    VectorTools::Linfty_norm,
                                    mask // & B_mask
  );
 
  Linfty_norm = Linfty_per_cell.linfty_norm();
}
 
template<int dim, int stage>
void
ProjectHgradToHdiv<dim, stage>::project_exact_solution_fcn()
{
  Constants::QuadratureTableVector<dim> qt(fe_Hdiv.degree - 1);
 
  AffineConstraints<double> constraints_empty;
  constraints_empty.close();
 
  VectorTools::project(MappingQ<dim>(mapping_degree),
                       dof_handler_Hdiv,
                       constraints_empty,
                       QGauss<dim>(qt.sim()),
                       *exact_solution,
                       projected_exact_solution);
}
 
} // namespace StaticScalarSolver
 
#endif