static_vector_solver_i.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 StaticVectorSolverI_H__
#define StaticVectorSolverI_H__
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/types.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_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
 
#include <deal.II/fe/fe_nedelec.h>
#include <deal.II/fe/fe_values.h>
 
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>
 
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/vector_tools_common.h>
#include <deal.II/numerics/vector_tools_project.h>
 
#include <fstream>
#include <iomanip>
#include <ios>
#include <iostream>
#include <map>
 
#include "constants.hpp"
#include "settings.hpp"
#include "static_vector_input.hpp"
 
#define VE scratch_data.ve
 
#define TMR(__name) TimerOutput::Scope timer_section(timer, __name)
 
using namespace dealii;
 
namespace StaticVectorSolver {
template<int dim, int stage = 1>
class Solver1
{
public:
  Solver1() = delete;
 
  Solver1(unsigned int p,
          unsigned int mapping_degree,
          unsigned int type_of_pde_rhs = 3,
          double eta_squared = 0.0,
          std::string fname = "data",
          const Function<dim>* exact_solution = nullptr,
          bool print_time_tables = false,
          bool project_exact_solution = false,
          bool write_higher_order_cells = false)
    : fe(p)
    , mapping_degree(mapping_degree)
    , type_of_pde_rhs(type_of_pde_rhs)
    , eta_squared(eta_squared)
    , fname(fname)
    , exact_solution(exact_solution)
    , print_time_tables(print_time_tables)
    , project_exact_solution(project_exact_solution)
    , write_higher_order_cells(write_higher_order_cells)
  {
    Assert(((dim == 2) || (dim == 3)), ExcInternalError());
    Assert(p < 5, ExcInternalError());
    Assert(type_of_pde_rhs < 4, ExcInternalError());
  }
 
  virtual void make_mesh() = 0;
 
  virtual void fill_dirichlet_stack() = 0;
 
  virtual void solve() = 0;
 
  void setup();
 
  void assemble();
 
  void compute_error_norms();
 
  void project_exact_solution_fcn();
 
  void save() const;
 
  void save_matrix_and_rhs_to_csv(std::string fname) const;
 
  void clear()
  {
    system_matrix.clear();
    system_rhs.reinit(0);
  }
 
  const Triangulation<dim>& get_tria() const { return triangulation; }
 
  const DoFHandler<dim>& get_dof_handler() const { return dof_handler; }
 
  const Vector<double>& get_solution() const { return solution; }
 
  unsigned int get_n_cells() const
  {
    return static_cast<unsigned int>(triangulation.n_active_cells());
  }
 
  unsigned int get_n_dofs() const
  {
    return static_cast<unsigned int>(dof_handler.n_dofs());
  }
 
  unsigned int get_rhs_type() const { return type_of_pde_rhs; }
 
  double get_L2_norm() const { return L2_norm; }
 
  double get_Linfty_norm() const { return Linfty_norm; }
 
  unsigned int get_mapping_degree() const { return mapping_degree; }
 
  void run()
  {
    TimerOutput::OutputFrequency tf =
      (print_time_tables) ? TimerOutput::summary : TimerOutput::never;
 
    TimerOutput timer(std::cout, tf, TimerOutput::cpu_and_wall_times_grouped);
 
    {
      TMR("Make mesh");
      make_mesh();
    }
    {
      TMR("Fill Dirichlet stack");
      fill_dirichlet_stack();
    }
    {
      TMR("Setup");
      setup();
    }
    {
      TMR("Assemble");
      assemble();
    }
    {
      TMR("Solve");
      solve();
    }
 
    if (exact_solution) {
      if (project_exact_solution) {
        TMR("Project exact solution");
        project_exact_solution_fcn();
      }
 
      {
        TMR("Compute error norms");
        compute_error_norms();
      }
    }
 
    {
      TMR("Save");
      save();
    }
  };
 
  virtual ~Solver1() = default;
 
protected:
  std::map<types::boundary_id, const Function<dim>*> dirichlet_stack;
 
  Triangulation<dim> triangulation;
 
  const FE_Nedelec<dim> fe;
 
  DoFHandler<dim> dof_handler;
 
  Vector<double> solution;
 
  Vector<double> projected_exact_solution;
 
  AffineConstraints<double> constraints;
 
  SparsityPattern sparsity_pattern;
 
  SparseMatrix<double> system_matrix;
 
  Vector<double> system_rhs;
 
  double L2_norm;
 
  double Linfty_norm;
 
private:
  const unsigned int mapping_degree;
  const unsigned int type_of_pde_rhs;
  const double eta_squared;
  const std::string fname;
  const Function<dim>* exact_solution;
  const bool print_time_tables;
  const bool project_exact_solution;
  const bool write_higher_order_cells;
 
  Vector<float> L2_per_cell;
  Vector<float> Linfty_per_cell;
 
  // ----------------------------------------------------------------------------
  // 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.
  // ----------------------------------------------------------------------------
  struct AssemblyScratchData
  {
    AssemblyScratchData(const FiniteElement<dim>& fe,
                        unsigned int mapping_degree,
                        unsigned int type_of_pde_rhs,
                        double eta_squared);
 
    AssemblyScratchData(const AssemblyScratchData& scratch_data);
 
    TheCoefficient<dim, stage> the_coefficient;
    PdeRhs<dim, stage> pde_rhs;
    Gamma<dim, stage> gamma;
    RobinRhs<dim, stage> robin_rhs;
    FreeSurfaceCurrent<dim, stage> free_surface_current;
 
    MappingQ<dim> mapping;
    Constants::QuadratureTableVector<dim> qt;
    FEValues<dim> fe_values;
    FEFaceValues<dim> fe_face_values;
 
    const unsigned int dofs_per_cell;
    const unsigned int n_q_points;
    const unsigned int n_q_points_face;
 
    std::vector<double> the_coefficient_list;
    std::vector<Tensor<1, dim>> pde_rhs_list;
    std::vector<Tensor<1, dim>> pde_rhs_list_face;
    std::vector<double> gamma_list;
    std::vector<Tensor<1, dim>> robin_rhs_list;
    std::vector<Tensor<1, dim>> free_surface_current_list;
 
    const unsigned int type_of_pde_rhs;
    const double eta_squared;
    const FEValuesExtractors::Vector ve;
    bool do_robin;
    bool do_K;
    bool do_T_on_boundary;
  };
 
  struct AssemblyCopyData
  {
    FullMatrix<double> cell_matrix;
    Vector<double> cell_rhs;
    std::vector<types::global_dof_index> local_dof_indices;
  };
 
  void system_matrix_local(
    const typename DoFHandler<dim>::active_cell_iterator& cell,
    AssemblyScratchData& scratch_data,
    AssemblyCopyData& copy_data);
 
  void copy_local_to_global(const AssemblyCopyData& copy_data);
  //-----------------------------------------------------------------------------
  //-----------------------------------------------------------------------------
  //-----------------------------------------------------------------------------
};
 
template<int dim, int stage>
void
Solver1<dim, stage>::setup()
{
  dof_handler.reinit(triangulation);
  dof_handler.distribute_dofs(fe);
 
  constraints.clear();
  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-but-set-variable"
  for (auto item : dirichlet_stack) {
    Assert(item.first % 2 == 1, ExcInternalError());
  }
#pragma GCC diagnostic pop
 
  for (auto item : dirichlet_stack)
    VectorTools::project_boundary_values_curl_conforming_l2(
      dof_handler,
      0,            // first vector component
      *item.second, // boundary function
      item.first,   // boundary id
      constraints,  // constraints
      MappingQ<dim>(mapping_degree));
 
  constraints.close();
 
  DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
  DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
 
  sparsity_pattern.copy_from(dsp);
  system_matrix.reinit(sparsity_pattern);
  solution.reinit(dof_handler.n_dofs());
  system_rhs.reinit(dof_handler.n_dofs());
 
  if (project_exact_solution)
    projected_exact_solution.reinit(dof_handler.n_dofs());
 
  if (exact_solution) {
    L2_per_cell.reinit(triangulation.n_active_cells());
    Linfty_per_cell.reinit(triangulation.n_active_cells());
  }
}
 
template<int dim, int stage>
void
Solver1<dim, stage>::assemble()
{
  WorkStream::run(
    dof_handler.begin_active(),
    dof_handler.end(),
    *this,
    &Solver1::system_matrix_local,
    &Solver1::copy_local_to_global,
    AssemblyScratchData(fe, mapping_degree, type_of_pde_rhs, eta_squared),
    AssemblyCopyData());
}
 
template<int dim, int stage>
Solver1<dim, stage>::AssemblyScratchData::AssemblyScratchData(
  const FiniteElement<dim>& fe,
  unsigned int mapping_degree,
  unsigned int type_of_pde_rhs,
  double eta_squared)
  : mapping(mapping_degree)
  , qt(fe.degree - 1)
  , fe_values(mapping,
              fe,
              QGauss<dim>(qt.sim()),
              update_gradients | update_values | update_quadrature_points |
                update_JxW_values)
  , fe_face_values(mapping,
                   fe,
                   QGauss<dim - 1>(qt.sim()),
                   update_values | update_normal_vectors |
                     update_quadrature_points | update_JxW_values)
  , dofs_per_cell(fe_values.dofs_per_cell)
  , n_q_points(fe_values.get_quadrature().size())
  , n_q_points_face(fe_face_values.get_quadrature().size())
  , the_coefficient_list(n_q_points)
  , pde_rhs_list(n_q_points, Tensor<1, dim>())
  , pde_rhs_list_face(n_q_points_face, Tensor<1, dim>())
  , gamma_list(n_q_points_face)
  , robin_rhs_list(n_q_points_face, Tensor<1, dim>())
  , free_surface_current_list(n_q_points_face, Tensor<1, dim>())
  , type_of_pde_rhs(type_of_pde_rhs)
  , eta_squared(eta_squared)
  , ve(0)
{
}
 
template<int dim, int stage>
Solver1<dim, stage>::AssemblyScratchData::AssemblyScratchData(
  const AssemblyScratchData& scratch_data)
  : mapping(scratch_data.mapping.get_degree())
  , qt(scratch_data.qt)
  , fe_values(mapping,
              scratch_data.fe_values.get_fe(),
              scratch_data.fe_values.get_quadrature(),
              update_gradients | update_values | update_quadrature_points |
                update_JxW_values)
  , fe_face_values(mapping,
                   scratch_data.fe_face_values.get_fe(),
                   scratch_data.fe_face_values.get_quadrature(),
                   update_values | update_normal_vectors |
                     update_quadrature_points | update_JxW_values)
  , dofs_per_cell(fe_values.dofs_per_cell)
  , n_q_points(fe_values.get_quadrature().size())
  , n_q_points_face(fe_face_values.get_quadrature().size())
  , the_coefficient_list(n_q_points)
  , pde_rhs_list(n_q_points, Tensor<1, dim>())
  , pde_rhs_list_face(n_q_points_face, Tensor<1, dim>())
  , gamma_list(n_q_points_face)
  , robin_rhs_list(n_q_points_face, Tensor<1, dim>())
  , free_surface_current_list(n_q_points_face, Tensor<1, dim>())
  , type_of_pde_rhs(scratch_data.type_of_pde_rhs)
  , eta_squared(scratch_data.eta_squared)
  , ve(0)
{
}
 
template<int dim, int stage>
void
Solver1<dim, stage>::system_matrix_local(
  const typename DoFHandler<dim>::active_cell_iterator& cell,
  AssemblyScratchData& scratch_data,
  AssemblyCopyData& copy_data)
{
  // See the following boxes:
  // (1) Recipe for static vector solver in 3D
  // (2) Recipe for static vector solver in 2D
  // (3) Recipe for static vector solver in 3D (current vect. potential)
 
  // The comments below refer to these recipes by number, i.e., recipe (1),
  // recipe (2), and recipe (3).
 
  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.reinit(cell);
 
  scratch_data.the_coefficient.value_list(
    scratch_data.fe_values.get_quadrature_points(),
    cell->material_id(),
    cell->user_index(),
    scratch_data.the_coefficient_list);
 
  scratch_data.pde_rhs.value_list(
    scratch_data.fe_values.get_quadrature_points(),
    cell->material_id(),
    cell->user_index(),
    scratch_data.pde_rhs_list);
 
  // The curl of a Nedelec shape function is Tensor<1,3> in three dimensions and
  // Tensor<1,1> in two dimensions. It is somewhat difficult to fit this into
  // the "dim" class template paradigm of deal.II. Consequently, we have to use
  // the if (dim==2) and if (dim==3) filters or instantiate this function
  // template explicitly for two and three dimensions. We choose the former.
  //
  // If the current vector potential, T, is used in recipes (1) and (2), it must
  // be implemented by the same class that implements the right-hand side of the
  // curl-curl equation. Neglect the fact that the right-hand side is the curl
  // of the current vector potential, not the current vector potential itself.
  // So, "pde_rhs" is a misnomer if the current vector potential is used on the
  // right-hand side of the PDE. The current vector potential must be
  // Tensor<1,3> in three dimensions. In two dimensions it must fill the first
  // component of Tensor<1,2>. The second component is ignored:
  //
  // Tensor<1,2> a;
  // a[0] = T;
  // a[1] = whatever;
  //
  // The same holds for the free-current density, J_f, if the right-hand side of
  // the PDE is the curl of J_f, i.e., recipe (3).
 
  for (unsigned int q_index = 0; q_index < scratch_data.n_q_points; ++q_index) {
    for (unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i) {
      for (unsigned int j = 0; j < scratch_data.dofs_per_cell; ++j) {
        copy_data.cell_matrix(i, j) +=
          ( // Integral I_a1+I_a3 in recipes (1), (2), and (3).
            (1.0 / scratch_data.the_coefficient_list[q_index]) * // 1 / mu
              scratch_data.fe_values[VE].curl(i, q_index) *      // curl N_i
              scratch_data.fe_values[VE].curl(j, q_index)        // curl N_j
            + scratch_data.eta_squared *                         // eta^2
                scratch_data.fe_values[VE].value(i, q_index) *   // N_i
                scratch_data.fe_values[VE].value(j, q_index)     // N_j
            ) *
          scratch_data.fe_values.JxW(q_index); // dV (dS in 2D)
      }
 
      switch (scratch_data.type_of_pde_rhs) {
        case 0:
          // Integral I_b3 in recipes (1) and (2) with J_f=0.
          copy_data.cell_rhs(i) = 0.0;
          break;
        case 1:
          // Integral I_b3 in recipes (1) and (2).
          copy_data.cell_rhs(i) +=
            scratch_data.pde_rhs_list[q_index] *           // J_f
            scratch_data.fe_values[VE].value(i, q_index) * // N_i
            scratch_data.fe_values.JxW(q_index);           // dV (dS in 2D)
          break;
        case 2:
        case 3:
          if (dim == 2) { // Integral I_b3-1 in recipe (2)
            copy_data.cell_rhs(i) +=
              scratch_data.pde_rhs_list[q_index][0] *          // T
              scratch_data.fe_values[VE].curl(i, q_index)[0] * // curl_s N_i
              scratch_data.fe_values.JxW(q_index);             // dS
          } else if (dim == 3) { // Integral I_b3-1 in recipe (1) and (3).
            copy_data.cell_rhs(i) +=
              (scratch_data.pde_rhs_list[q_index][0] *
                 scratch_data.fe_values[VE].curl(i, q_index)[0] +
               scratch_data.pde_rhs_list[q_index][1] *
                 scratch_data.fe_values[VE].curl(i, q_index)[1] +
               scratch_data.pde_rhs_list[q_index][2] *
                 scratch_data.fe_values[VE].curl(i, q_index)[2]) *
              scratch_data.fe_values.JxW(
                q_index); // If recipe (1): T.(curl N_i)dV.
                          // If recipe (3): J_f.(curl N_i)dV.
            // So, the scratch_data.pde_rhs_list must contain the values of T,
            // or the values of J_f, depending on the context.
          } else {
            Assert(false, ExcInternalError());
          }
          break;
        default:
          Assert(false, ExcInternalError());
          break;
      }
    }
  }
 
  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f) {
    scratch_data.do_robin = ((cell->face(f)->at_boundary()) &&
                             (cell->face(f)->boundary_id() % 2 == 0) &&
                             (cell->face(f)->boundary_id() != 0));
 
    scratch_data.do_K =
      ((cell->user_index() > 0) && (cell->face(f)->user_index() > 0));
 
    scratch_data.do_T_on_boundary =
      ((cell->face(f)->at_boundary()) && (scratch_data.type_of_pde_rhs == 3));
 
    Assert(!(scratch_data.do_robin && scratch_data.do_K), ExcInternalError());
 
    if (scratch_data.do_robin || scratch_data.do_K ||
        scratch_data.do_T_on_boundary) {
 
      scratch_data.fe_face_values.reinit(cell, f);
 
      if (scratch_data.do_robin) {
        scratch_data.gamma.value_list(
          scratch_data.fe_face_values.get_quadrature_points(),
          scratch_data.fe_face_values.get_normal_vectors(),
          cell->face(f)->boundary_id(),
          cell->material_id(),
          cell->user_index(),
          cell->face(f)->user_index(),
          scratch_data.gamma_list);
 
        scratch_data.robin_rhs.value_list(
          scratch_data.fe_face_values.get_quadrature_points(),
          scratch_data.fe_face_values.get_normal_vectors(),
          cell->face(f)->boundary_id(),
          cell->material_id(),
          cell->user_index(),
          cell->face(f)->user_index(),
          scratch_data.robin_rhs_list);
      }
 
      if (scratch_data.do_K) {
        scratch_data.free_surface_current.value_list(
          scratch_data.fe_face_values.get_quadrature_points(),
          scratch_data.fe_face_values.get_normal_vectors(),
          cell->material_id(),
          cell->user_index(),
          cell->face(f)->user_index(),
          scratch_data.free_surface_current_list);
      }
 
      if (scratch_data.do_T_on_boundary) {
        scratch_data.pde_rhs.value_list(
          scratch_data.fe_face_values.get_quadrature_points(),
          cell->material_id(),
          cell->user_index(),
          scratch_data.pde_rhs_list_face);
      }
 
      for (unsigned int q_index_face = 0;
           q_index_face < scratch_data.n_q_points_face;
           ++q_index_face) {
        for (unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i) {
          if (scratch_data.do_robin) {
            for (unsigned int j = 0; j < scratch_data.dofs_per_cell; ++j) {
              if (dim == 2) {
                // Integral I_a2 in recipe (2).
                copy_data.cell_matrix(i, j) +=
                  scratch_data.gamma_list[q_index_face] * // gamma
                  (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
                     scratch_data.fe_face_values[VE].value(i, q_index_face)[1] -
                   scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
                     scratch_data.fe_face_values[VE].value(i,
                                                           q_index_face)[0]) *
                  (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
                     scratch_data.fe_face_values[VE].value(j, q_index_face)[1] -
                   scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
                     scratch_data.fe_face_values[VE].value(j,
                                                           q_index_face)[0]) *
                  scratch_data.fe_face_values.JxW(q_index_face); // dl
              } else if (dim == 3) {
                // Integral I_a2 in recipe (1).
                copy_data.cell_matrix(i, j) +=
                  scratch_data.gamma_list[q_index_face] * // gamma
                  ((scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
                      scratch_data.fe_face_values[VE].value(i,
                                                            q_index_face)[2] -
                    scratch_data.fe_face_values.normal_vector(q_index_face)[2] *
                      scratch_data.fe_face_values[VE].value(i,
                                                            q_index_face)[1]) *
                     (scratch_data.fe_face_values.normal_vector(
                        q_index_face)[1] *
                        scratch_data.fe_face_values[VE].value(j,
                                                              q_index_face)[2] -
                      scratch_data.fe_face_values.normal_vector(
                        q_index_face)[2] *
                        scratch_data.fe_face_values[VE].value(
                          j, q_index_face)[1]) +
                   (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
                      scratch_data.fe_face_values[VE].value(i,
                                                            q_index_face)[2] -
                    scratch_data.fe_face_values.normal_vector(q_index_face)[2] *
                      scratch_data.fe_face_values[VE].value(i,
                                                            q_index_face)[0]) *
                     (scratch_data.fe_face_values.normal_vector(
                        q_index_face)[0] *
                        scratch_data.fe_face_values[VE].value(j,
                                                              q_index_face)[2] -
                      scratch_data.fe_face_values.normal_vector(
                        q_index_face)[2] *
                        scratch_data.fe_face_values[VE].value(
                          j, q_index_face)[0]) +
                   (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
                      scratch_data.fe_face_values[VE].value(i,
                                                            q_index_face)[1] -
                    scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
                      scratch_data.fe_face_values[VE].value(i,
                                                            q_index_face)[0]) *
                     (scratch_data.fe_face_values.normal_vector(
                        q_index_face)[0] *
                        scratch_data.fe_face_values[VE].value(j,
                                                              q_index_face)[1] -
                      scratch_data.fe_face_values.normal_vector(
                        q_index_face)[1] *
                        scratch_data.fe_face_values[VE].value(
                          j, q_index_face)[0])) *
                  scratch_data.fe_face_values.JxW(
                    q_index_face); // (n x N_i) . (n x N_j) dS
              } else {
                Assert(false, ExcInternalError());
              }
            }
 
            // Integral I_b1 in recipes (1), (2), and (3).
            copy_data.cell_rhs(i) +=
              -scratch_data.robin_rhs_list[q_index_face] *             // Q
              scratch_data.fe_face_values[VE].value(i, q_index_face) * // N_i
              scratch_data.fe_face_values.JxW(q_index_face); // dS (dl in 2D)
          } // if (scratch_data.do_robin)
 
          if (scratch_data.do_K) {
            // Integral I_b2 in recipes (1) and (2).
            copy_data.cell_rhs(i) +=
              scratch_data.free_surface_current_list[q_index_face] *   // K_f
              scratch_data.fe_face_values[VE].value(i, q_index_face) * // N_i
              scratch_data.fe_face_values.JxW(q_index_face); // dS (dl in 2D)
          }
          
          if (scratch_data.do_T_on_boundary) {
            if (dim == 2) {
              // Integral I_b3-2 in recipe (2).
              copy_data.cell_rhs(i) -=
                scratch_data.pde_rhs_list_face[q_index_face][0] * // T
                (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
                   scratch_data.fe_face_values[VE].value(i, q_index_face)[1] -
                 scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
                   scratch_data.fe_face_values[VE].value(i,
                                                         q_index_face)[0] //   S
                 ) *                                           // n x N_i
                scratch_data.fe_face_values.JxW(q_index_face); // dl
            } else if (dim == 3) {
              // Integral I_b3-2 in recipes (1) and (3).
              copy_data.cell_rhs(i) -=
                (scratch_data.pde_rhs_list_face[q_index_face][0] *
                   (scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
                      scratch_data.fe_face_values[VE].value(i,
                                                            q_index_face)[2] -
                    scratch_data.fe_face_values.normal_vector(q_index_face)[2] *
                      scratch_data.fe_face_values[VE].value(i,
                                                            q_index_face)[1]) -
                 scratch_data.pde_rhs_list_face[q_index_face][1] *
                   (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
                      scratch_data.fe_face_values[VE].value(i,
                                                            q_index_face)[2] -
                    scratch_data.fe_face_values.normal_vector(q_index_face)[2] *
                      scratch_data.fe_face_values[VE].value(i,
                                                            q_index_face)[0]) +
                 scratch_data.pde_rhs_list_face[q_index_face][2] *
                   (scratch_data.fe_face_values.normal_vector(q_index_face)[0] *
                      scratch_data.fe_face_values[VE].value(i,
                                                            q_index_face)[1] -
                    scratch_data.fe_face_values.normal_vector(q_index_face)[1] *
                      scratch_data.fe_face_values[VE].value(i,
                                                            q_index_face)[0])) *
                scratch_data.fe_face_values.JxW(
                  q_index_face); // If recipe (1): T.(n x N_i)dS.
            }                    // If recipe (3): J_f.(n x N_i)dS.
              // So, the scratch_data.pde_rhs_list_face must contain the values
              // of T or the values of J_f, depending on the context.
            else {
              Assert(false, ExcInternalError());
            }
          }
        }
      }
    }
  }
  cell->get_dof_indices(copy_data.local_dof_indices);
}
 
template<int dim, int stage>
void
Solver1<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
Solver1<dim, stage>::compute_error_norms()
{
  Weight<dim, stage> weight;
  const Function<dim, double>* mask = &weight;
 
  Constants::QuadratureTableVector<dim> qt(dof_handler.get_fe().degree - 1);
  QGauss<dim> quadrature(qt.enorm());
 
  VectorTools::integrate_difference(MappingQ<dim>(mapping_degree),
                                    dof_handler,
                                    solution,
                                    *exact_solution,
                                    L2_per_cell,
                                    quadrature,
                                    VectorTools::L2_norm,
                                    mask);
 
  L2_norm = VectorTools::compute_global_error(
    triangulation, L2_per_cell, VectorTools::L2_norm);
 
  VectorTools::integrate_difference(MappingQ<dim>(mapping_degree),
                                    dof_handler,
                                    solution,
                                    *exact_solution,
                                    Linfty_per_cell,
                                    QGauss<dim>(1),
                                    VectorTools::Linfty_norm,
                                    mask);
 
  Linfty_norm = Linfty_per_cell.linfty_norm();
}
 
template<int dim, int stage>
void
Solver1<dim, stage>::project_exact_solution_fcn()
{
  Constants::QuadratureTableVector<dim> qt(fe.degree - 1);
 
  AffineConstraints<double> constraints_empty;
  constraints_empty.close();
 
  VectorTools::project(MappingQ<dim>(mapping_degree),
                       dof_handler,
                       constraints_empty,
                       QGauss<dim>(qt.sim()),
                       *exact_solution,
                       projected_exact_solution);
}
 
template<int dim, int stage>
void
Solver1<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, solution, solution_names, interpretation);
 
  if (project_exact_solution && exact_solution) {
    std::vector<std::string> solution_names_ex(dim, "VectorFieldExact");
 
    data_out.add_data_vector(
      dof_handler, projected_exact_solution, solution_names_ex, interpretation);
  }
 
  if (exact_solution) {
    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.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
Solver1<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();
}
 
} // namespace StaticVectorSolver
 
#endif