step-97.h

Go to the documentation of this file.

 
 *   
 * @endcode
 * 
 * The following three get-functions are used to channel the mesh and the
 * solution to the next solver.
 * 
 * @code
 *       const Triangulation<3> &get_tria() const
 *       {
 *         return triangulation;
 *       }
 *       const DoFHandler<3> &get_dof_handler() const
 *       {
 *         return dof_handler;
 *       }
 *       const Vector<double> &get_solution() const
 *       {
 *         return solution;
 *       }
 *   
 *     private:
 * @endcode
 * 
 * The following data members are typical for all deal.II simulations:
 * triangulation, finite elements, dof handlers, etc. The constraints
 * are used to enforce the Dirichlet boundary conditions. The names of the
 * data members are self-explanatory.
 * 
 
 * 
 * 
 * @code
 *       Triangulation<3> triangulation;
 *   
 *       const FE_Nedelec<3> fe;
 *       DoFHandler<3>       dof_handler;
 *       Vector<double>      solution;
 *   
 *       SparseMatrix<double> system_matrix;
 *       Vector<double>       system_rhs;
 *   
 *       AffineConstraints<double> constraints;
 *       SparsityPattern           sparsity_pattern;
 *   
 *       SphericalManifold<3> sphere;
 *   
 *       const unsigned int refinement_parameter;
 *       const unsigned int mapping_degree;
 *       const double       eta_squared;
 *       const std::string  fname;
 *   
 * @endcode
 * 
 * The program utilizes the WorkStream technology. The @ref step_9 "step-9" tutorial
 * does a much better job of explaining the workings of WorkStream.
 * Reading the @ref workstream_paper "WorkStream paper" is recommended.
 * The following structures and functions are related to WorkStream.
 * 
 * @code
 *       struct AssemblyScratchData
 *       {
 *         AssemblyScratchData(const FiniteElement<3> &fe,
 *                             const double            eta_squared,
 *                             const unsigned int      mapping_degree);
 *   
 *         AssemblyScratchData(const AssemblyScratchData &scratch_data);
 *   
 *         const FreeCurrentDensity Jf;
 *   
 *         MappingQ<3> mapping;
 *         FEValues<3> fe_values;
 *   
 *         const unsigned int dofs_per_cell;
 *         const unsigned int n_q_points;
 *   
 *         std::vector<Tensor<1, 3>> Jf_list;
 *   
 *         const double eta_squared;
 *       };
 *   
 *       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<3>::active_cell_iterator &cell,
 *         AssemblyScratchData                                &scratch_data,
 *         AssemblyCopyData                                   &copy_data);
 *   
 *       void copy_local_to_global(const AssemblyCopyData &copy_data);
 *     };
 *   
 *     Solver::Solver(const unsigned int p,
 *                    const unsigned int r,
 *                    const unsigned int mapping_degree,
 *                    const double       eta_squared,
 *                    const std::string &fname)
 *       : fe(p)
 *       , refinement_parameter(r)
 *       , mapping_degree(mapping_degree)
 *       , eta_squared(eta_squared)
 *       , fname(fname)
 *     {}
 *   
 * @endcode
 * 
 * The following function loads the mesh, assigns material IDs to all cells,
 * and attaches the spherical manifold to the mesh. The material IDs are
 * assigned on the basis of the distance from the center of a cell to the
 * origin. The spherical manifold is attached to a face if all vertices of
 * the face are at the same distance from the origin provided the cell is
 * outside the cube in the center of the mesh, see mesh description in the
 * introduction.
 * 
 * @code
 *     void Solver::make_mesh()
 *     {
 *       GridIn<3> gridin;
 *   
 *       gridin.attach_triangulation(triangulation);
 *       std::ifstream ifs("sphere_r" + std::to_string(refinement_parameter) +
 *                         ".msh");
 *       gridin.read_msh(ifs);
 *   
 *       triangulation.reset_all_manifolds();
 *   
 *       for (auto cell : triangulation.active_cell_iterators())
 *         {
 *           cell->set_material_id(
 *             Settings::material_id_free_space); // The cell is in free space.
 *   
 *           if ((cell->center().norm() > Settings::a1) &&
 *               (cell->center().norm() < Settings::b1))
 *             cell->set_material_id(
 *               Settings::material_id_core); // The cell is inside the core.
 *   
 *           if ((cell->center().norm() > Settings::a2) &&
 *               (cell->center().norm() < Settings::b2))
 *             cell->set_material_id(
 *               Settings::material_id_free_current); /* The cell is inside the Jf
 *                                                       region. */
 *   
 *           for (unsigned int f = 0; f < cell->n_faces(); f++)
 *             {
 *               double dif_norm = 0.0;
 *               for (unsigned int v = 1; v < cell->face(f)->n_vertices(); v++)
 *                 dif_norm += std::abs(cell->face(f)->vertex(0).norm() -
 *                                      cell->face(f)->vertex(v).norm());
 *   
 *               if ((dif_norm < Settings::eps) &&
 *                   (cell->center().norm() > Settings::d1))
 *                 cell->face(f)->set_all_manifold_ids(1);
 *             }
 *         }
 *   
 *       triangulation.set_manifold(1, sphere);
 *     }
 *   
 * @endcode
 * 
 * The following function initializes the dofs, applies the Dirichlet
 * boundary condition, and initializes the vectors and matrices. The first
 * two lines of the code initialize the dof handler and distribute the dofs.
 * The segment of the code in between `constraints.clear()` and
 * `constraints.close()` applies the homogeneous Dirichlet boundary condition.
 * As discussed in the introduction, the Dirichlet boundary condition is an
 * essential condition and must be enforced by constraining the system matrix.
 * This segment of the code does the constraining. The rest of the function
 * arranges the dofs in a sparsity pattern and initializes the system matrices
 * and vectors.
 * 
 * @code
 *     void Solver::setup()
 *     {
 *       dof_handler.reinit(triangulation);
 *       dof_handler.distribute_dofs(fe);
 *   
 *       constraints.clear();
 *   
 *       DoFTools::make_hanging_node_constraints(dof_handler, constraints);
 *   
 *       VectorTools::project_boundary_values_curl_conforming_l2(
 *         dof_handler,
 *         0, // The first vector component.
 *         Functions::ZeroFunction<3>(3),
 *         Settings::boundary_id_infinity,
 *         constraints,
 *         MappingQ<3>(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());
 *     }
 *   
 * @endcode
 * 
 * Formally, the following function assembles the system of linear equations.
 * In reality, however, it just spells all the magic words to get the
 * WorkStream going. The interesting part, i.e., the actual assembling of the
 * system matrix and the right-hand side, happens below in the function
 * `Solver::system_matrix_local()`.
 * 
 * @code
 *     void Solver::assemble()
 *     {
 *       WorkStream::run(dof_handler.begin_active(),
 *                       dof_handler.end(),
 *                       *this,
 *                       &Solver::system_matrix_local,
 *                       &Solver::copy_local_to_global,
 *                       AssemblyScratchData(fe, eta_squared, mapping_degree),
 *                       AssemblyCopyData());
 *     }
 *   
 * @endcode
 * 
 * The following two constructors initialize scratch data from the input
 * parameters and from another object of the same type, i.e., a copy
 * constructor.
 * 
 * @code
 *     Solver::AssemblyScratchData::AssemblyScratchData(
 *       const FiniteElement<3> &fe,
 *       const double            eta_squared,
 *       const unsigned int      mapping_degree)
 *       : Jf()
 *       , mapping(mapping_degree)
 *       , fe_values(mapping,
 *                   fe,
 *                   QGauss<3>(fe.degree + 2),
 *                   update_gradients | update_values | update_quadrature_points |
 *                     update_JxW_values)
 *       , dofs_per_cell(fe_values.dofs_per_cell)
 *       , n_q_points(fe_values.get_quadrature().size())
 *       , Jf_list(n_q_points, Tensor<1, 3>())
 *       , eta_squared(eta_squared)
 *     {}
 *   
 *     Solver::AssemblyScratchData::AssemblyScratchData(
 *       const AssemblyScratchData &scratch_data)
 *       : Jf()
 *       , mapping(scratch_data.mapping.get_degree())
 *       , 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)
 *       , dofs_per_cell(fe_values.dofs_per_cell)
 *       , n_q_points(fe_values.get_quadrature().size())
 *       , Jf_list(n_q_points, Tensor<1, 3>())
 *       , eta_squared(scratch_data.eta_squared)
 *     {}
 *   
 * @endcode
 * 
 * The following function assembles a fraction of the system matrix and the
 * system right-hand side related to a single cell. These fractions are
 * `copy_data.cell_matrix` and `copy_data.cell_rhs`. They are copied into
 * the system matrix, @f$A_{ij}@f$, and the right-hand side, @f$b_i@f$, by the
 * function `Solver::copy_local_to_global()`.
 *   
 
 * 
 * In the first four statements of the function we reinitialize the matrices
 * and vectors related to the current cell and compute the finite element
 * values. Next, we compute the free-current density, @f$\vec{J}_f@f$, at the
 * quadrature points by calling `scratch_data.Jf.value_list`. After that we
 * declare a vector values extractor, `ve`, and compute the
 * [components](@ref Step97_Numerical_Recipe_T) of the cell matrix and cell
 * right-hand side in the nested `for` loops. The labels of the integrals are
 * the same as in the introduction to this tutorial. In the last line of the
 * function we query the dof indices on the current cell and store them in the
 * copy data structure, so we know to which locations of the system matrix and
 * right-hand side the components of the cell matrix and cell right-hand side
 * must be copied.
 * 
 * @code
 *     void Solver::system_matrix_local(
 *       const typename DoFHandler<3>::active_cell_iterator &cell,
 *       AssemblyScratchData                                &scratch_data,
 *       AssemblyCopyData                                   &copy_data)
 *     {
 *       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.Jf.value_list(scratch_data.fe_values.get_quadrature_points(),
 *                                  cell->material_id(),
 *                                  scratch_data.Jf_list);
 *   
 *       const FEValuesExtractors::Vector ve(0);
 *   
 *       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.
 *                     (scratch_data.fe_values[ve].curl(i,
 *                                                      q_index) * // curl phi_i(x_q)
 *                        scratch_data.fe_values[ve].curl(j,
 *                                                        q_index) // curl phi_j(x_q)
 *                      +
 *                      scratch_data.eta_squared * // eta^2
 *                        scratch_data.fe_values[ve].value(i,
 *                                                         q_index) *  // phi_i(x_q)
 *                        scratch_data.fe_values[ve].value(j, q_index) // phi_j(x_q)
 *                      ) *
 *                     scratch_data.fe_values.JxW(q_index); // dx
 *                 }
 *               copy_data.cell_rhs(i) += // Integral I_b3-1.
 *                 (scratch_data.Jf_list[q_index] *
 *                  scratch_data.fe_values[ve].curl(i, q_index)) *
 *                 scratch_data.fe_values.JxW(
 *                   q_index); // J_f(x_q).(curl phi_i(x_q))dx.
 *             }
 *         }
 *   
 *       cell->get_dof_indices(copy_data.local_dof_indices);
 *     }
 *   
 * @endcode
 * 
 * The following function copies the components of a cell matrix and a cell
 * right-hand side into the system matrix, @f$A_{ij}@f$, and the system right-hand
 * side, @f$b_i@f$.
 * 
 * @code
 *     void Solver::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);
 *     }
 *   
 * @endcode
 * 
 * The following function solves the system of linear equations. In theory,
 * the CG solver can solve an @f$m \times m@f$ system of linear equations in at
 * most @f$m@f$ steps. Accordingly, we set the maximum number of iteration steps
 * to `system_rhs.size()`. The stopping condition is
 * \f[\|\boldsymbol{b} - \boldsymbol{A}\boldsymbol{c}\| < 10^{-6}
 * \|\boldsymbol{b}\|.\f] As soon as we use constraints, we must not forget to
 * distribute them.
 * 
 * @code
 *     void Solver::solve()
 *     {
 *       SolverControl control(system_rhs.size(), 1.0e-6 * system_rhs.l2_norm());
 *   
 *       GrowingVectorMemory<Vector<double>> memory;
 *       SolverCG<Vector<double>>            cg(control, memory);
 *   
 *       PreconditionSSOR<SparseMatrix<double>> preconditioner;
 *       preconditioner.initialize(system_matrix, 1.2);
 *   
 *       cg.solve(system_matrix, solution, system_rhs, preconditioner);
 *   
 *       constraints.distribute(solution);
 *     }
 *   
 * @endcode
 * 
 * The following function saves the computed current vector potential into a
 * `.vtu` file.
 * 
 * @code
 *     void Solver::output_results() const
 *     {
 *       const std::vector<std::string> solution_names(3, "VectorField");
 *       const std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *         interpretation(3,
 *                        DataComponentInterpretation::component_is_part_of_vector);
 *   
 *       DataOut<3> data_out;
 *   
 *       data_out.add_data_vector(dof_handler,
 *                                solution,
 *                                solution_names,
 *                                interpretation);
 *   
 *       DataOutBase::VtkFlags flags;
 *       flags.write_higher_order_cells = true;
 *       data_out.set_flags(flags);
 *   
 *       const MappingQ<3> mapping(mapping_degree);
 *   
 *       data_out.build_patches(mapping,
 *                              fe.degree + 2,
 *                              DataOut<3>::CurvedCellRegion::curved_inner_cells);
 *   
 *       std::ofstream ofs(fname + ".vtu");
 *       data_out.write_vtu(ofs);
 *     }
 *   
 * @endcode
 * 
 * The `run` function below aggregates all the computation steps in the
 * correct order.
 * 
 * @code
 *     void Solver::run()
 *     {
 *       make_mesh();
 *       setup();
 *       assemble();
 *       solve();
 *       output_results();
 *       clear();
 *     }
 *   } // namespace SolverT
 *   
 * @endcode
 * 
 * 
 * <a name="step_97-SolverA"></a> 
 * <h3>Solver - A</h3>
 * 
 
 * 
 * The following namespace contains all the code related to the computation of
 * the magnetic vector potential, @f$\vec{A}@f$. The main difference between this
 * solver and the solver for the current vector potential, @f$\vec{T}@f$, is in how
 * the information on the source is fed to respective solvers. The solver for
 * @f$\vec{T}@f$ is fed data sampled from the analytical closed-form expression for
 * @f$\vec{J}_f@f$. The solver for @f$\vec{A}@f$ is fed a field function, i.e., a
 * numerically computed current vector potential, @f$\vec{T}@f$.
 * 
 * @code
 *   namespace SolverA
 *   {
 * @endcode
 * 
 * The following class describes the permeability in the entire problem
 * domain. The permeability is [given](@ref Step97_Equation_MU) by the
 * definition of the problem, see the introduction.
 * 
 * @code
 *     class Permeability
 *     {
 *     public:
 *       void value_list(const types::material_id mid,
 *                       std::vector<double>     &values) const
 *       {
 *         if ((mid == Settings::material_id_free_space) ||
 *             (mid == Settings::material_id_free_current))
 *           std::fill(values.begin(), values.end(), Settings::mu_0);
 *   
 *         if (mid == Settings::material_id_core)
 *           std::fill(values.begin(), values.end(), Settings::mu_1);
 *       }
 *     };
 *   
 * @endcode
 * 
 * The following class describes the parameter @f$\gamma@f$ in the Robin boundary
 * condition. As soon as it is evaluated on the boundary, the permeability
 * equals to that of free space. Therefore, we evaluate the parameter gamma as
 * \f[
 * \gamma = \dfrac{1}{\mu_0 r}.
 * \f]
 * 
 * @code
 *     class Gamma
 *     {
 *     public:
 *       void value_list(const std::vector<Point<3>> &r,
 *                       std::vector<double>         &values) const
 *       {
 *         Assert(r.size() == values.size(),
 *                ExcDimensionMismatch(r.size(), values.size()));
 *   
 *         for (unsigned int i = 0; i < values.size(); i++)
 *           values[i] = 1.0 / (Settings::mu_0 * r[i].norm());
 *       }
 *     };
 *   
 * @endcode
 * 
 * The following class implements the solver that minimizes the
 * [functional](@ref Step97_Functional_A) @f$F(\vec{A})@f$. The numerically
 * computed current vector potential, @f$\vec{T}@f$, is fed to this solver by
 * means of the input parameters `dof_handler_T` and `solution_T`. Moreover,
 * this solver reuses the mesh on which @f$\vec{T}@f$ has been computed. The
 * reference to the mesh is passed via the input parameter `triangulation_T`.
 * 
 * @code
 *     class Solver
 *     {
 *     public:
 *       Solver() = delete;
 *       Solver(const unsigned int      p, // Degree of the Nedelec finite elements.
 *              const unsigned int      mapping_degree,
 *              const Triangulation<3> &triangulation_T,
 *              const DoFHandler<3>    &dof_handler_T,
 *              const Vector<double>   &solution_T,
 *              const double            eta_squared = 0.0,
 *              const std::string      &fname       = "data");
 *   
 *       void setup();                // Initializes dofs, vectors, matrices.
 *       void assemble();             // Assembles the system of linear equations.
 *       void solve();                // Solves the system of linear equations.
 *       void output_results() const; // Saves computed A into a vtu file.
 *       void clear()                 // Clears the memory for the next solver.
 *       {
 *         system_matrix.clear();
 *         system_rhs.reinit(0);
 *       }
 *       void run(); /* Executes the last five functions in the proper order
 *                      and measures the execution time for each function. */
 *   
 *       const DoFHandler<3> &get_dof_handler() const
 *       {
 *         return dof_handler;
 *       }
 *       const Vector<double> &get_solution() const
 *       {
 *         return solution;
 *       }
 *   
 *     private:
 *       const Triangulation<3> &triangulation_T;
 *       const DoFHandler<3>    &dof_handler_T;
 *       const Vector<double>   &solution_T;
 *   
 * @endcode
 * 
 * The following data members are typical for all deal.II simulations:
 * triangulation, finite elements, dof handlers, etc. The constraints
 * are used to enforce the Dirichlet boundary conditions. The names of the
 * data members are self-explanatory.
 * 
 * @code
 *       const FE_Nedelec<3> fe;
 *       DoFHandler<3>       dof_handler;
 *       Vector<double>      solution;
 *   
 *       SparseMatrix<double> system_matrix;
 *       Vector<double>       system_rhs;
 *   
 *       AffineConstraints<double> constraints;
 *       SparsityPattern           sparsity_pattern;
 *   
 *       const unsigned int mapping_degree;
 *       const double       eta_squared;
 *       const std::string  fname;
 *   
 * @endcode
 * 
 * This time we have two dof handlers, `dof_handler_T` for @f$\vec{T}@f$ and
 * `dof_handler` for @f$\vec{A}@f$. The WorkStream needs to walk through
 * the two dof handlers synchronously. For this purpose we will pair two
 * active cell iterators (one from `dof_handler_T`, another from
 * `dof_handler`). For that we need the `IteratorPair` type.
 * 
 * @code
 *       using IteratorTuple =
 *         std::tuple<typename DoFHandler<3>::active_cell_iterator,
 *                    typename DoFHandler<3>::active_cell_iterator>;
 *   
 *       using IteratorPair = SynchronousIterators<IteratorTuple>;
 *   
 * @endcode
 * 
 * The program utilizes the WorkStream technology. The @ref step_9 "step-9" tutorial
 * does a much better job of explaining the workings of WorkStream.
 * Reading the @ref workstream_paper "WorkStream paper" is recommended.
 * The following structures and functions are related to WorkStream.
 * 
 * @code
 *       struct AssemblyScratchData
 *       {
 *         AssemblyScratchData(const FiniteElement<3>     &fe,
 *                             const DoFHandler<3>        &dof_hand_T,
 *                             const Vector<double>       &dofs_T,
 *                             const unsigned int          mapping_degree,
 *                             const double                eta_squared,
 *                             const BoundaryConditionType boundary_condition_type);
 *   
 *         AssemblyScratchData(const AssemblyScratchData &scratch_data);
 *   
 *         const Permeability permeability;
 *         const Gamma        gamma;
 *   
 *         MappingQ<3>     mapping;
 *         FEValues<3>     fe_values;
 *         FEFaceValues<3> fe_face_values;
 *   
 *         FEValues<3> fe_values_T;
 *   
 *         const unsigned int dofs_per_cell;
 *         const unsigned int n_q_points;
 *         const unsigned int n_q_points_face;
 *   
 *         std::vector<double>       permeability_list;
 *         std::vector<double>       gamma_list;
 *         std::vector<Tensor<1, 3>> T_values;
 *   
 *         const DoFHandler<3>  &dof_hand_T;
 *         const Vector<double> &dofs_T;
 *   
 *         const double                eta_squared;
 *         const BoundaryConditionType boundary_condition_type;
 *       };
 *   
 *       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);
 *     };
 *   
 *     Solver::Solver(const unsigned int      p,
 *                    const unsigned int      mapping_degree,
 *                    const Triangulation<3> &triangulation_T,
 *                    const DoFHandler<3>    &dof_handler_T,
 *                    const Vector<double>   &solution_T,
 *                    const double            eta_squared,
 *                    const std::string      &fname)
 *       : triangulation_T(triangulation_T)
 *       , dof_handler_T(dof_handler_T)
 *       , solution_T(solution_T)
 *       , fe(p)
 *       , mapping_degree(mapping_degree)
 *       , eta_squared(eta_squared)
 *       , fname(fname)
 *     {}
 *   
 * @endcode
 * 
 * The following function initializes the dofs, applies the Dirichlet
 * boundary condition, and initializes the vectors and matrices. The first
 * two lines of the code initialize the dof handler and distribute the dofs.
 * The segment of the code in between `constraints.clear()` and
 * `constraints.close()` applies the homogeneous Dirichlet boundary condition.
 * As discussed in the introduction, the Dirichlet boundary condition is an
 * essential condition and must be enforced by constraining the system matrix.
 * This segment of the code does the constraining. The program can be
 * [switched](@ref Step97_TXT_BCSwitch) between the Dirichlet, Neumann, and
 * Robin boundary conditions. For this reason, we use the `if` statement to
 * make sure that the system matrix is constrained only if the Dirichlet
 * boundary condition is chosen by the user. The rest of the function arranges
 * the dofs in a sparsity pattern and initializes the system matrices and
 * vectors.
 * 
 * @code
 *     void Solver::setup()
 *     {
 *       dof_handler.reinit(triangulation_T);
 *       dof_handler.distribute_dofs(fe);
 *   
 *       constraints.clear();
 *   
 *       DoFTools::make_hanging_node_constraints(dof_handler, constraints);
 *   
 *       if (Settings::boundary_condition_type_A == Dirichlet)
 *         VectorTools::project_boundary_values_curl_conforming_l2(
 *           dof_handler,
 *           0, // The first vector component.
 *           Functions::ZeroFunction<3>(3),
 *           Settings::boundary_id_infinity,
 *           constraints,
 *           MappingQ<3>(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());
 *     }
 *   
 * @endcode
 * 
 * Formally, the following function assembles the system of linear equations.
 * In reality, however, it just spells all the magic words to get the
 * WorkStream going. The interesting part, i.e., the actual assembling of the
 * system matrix and the right-hand side happens below in the
 * `Solver::system_matrix_local` function. Note that this time the first two
 * input parameters to `WorkStream::run` are pairs of iterators, not
 * iterators themselves as per usual. Note also the order in which we package
 * the iterators: first the iterator of `dof_handler`, then the iterator of
 * the `dof_handler_T`. We will extract them in the same order.
 * 
 * @code
 *     void Solver::assemble()
 *     {
 *       WorkStream::run(IteratorPair({dof_handler.begin_active(),
 *                                     dof_handler_T.begin_active()}),
 *                       IteratorPair({dof_handler.end(), dof_handler_T.end()}),
 *                       *this,
 *                       &Solver::system_matrix_local,
 *                       &Solver::copy_local_to_global,
 *                       AssemblyScratchData(fe,
 *                                           dof_handler_T,
 *                                           solution_T,
 *                                           mapping_degree,
 *                                           eta_squared,
 *                                           Settings::boundary_condition_type_A),
 *                       AssemblyCopyData());
 *     }
 *   
 * @endcode
 * 
 * The following two constructors initialize scratch data from the input
 * parameters and from another object of the same type, i.e., a copy
 * constructor.
 * 
 * @code
 *     Solver::AssemblyScratchData::AssemblyScratchData(
 *       const FiniteElement<3>     &fe,
 *       const DoFHandler<3>        &dof_hand_T,
 *       const Vector<double>       &dofs_T,
 *       const unsigned int          mapping_degree,
 *       const double                eta_squared,
 *       const BoundaryConditionType boundary_condition_type)
 *       : permeability()
 *       , gamma()
 *       , mapping(mapping_degree)
 *       , fe_values(mapping,
 *                   fe,
 *                   QGauss<3>(fe.degree + 2),
 *                   update_gradients | update_values | update_JxW_values)
 *       , fe_face_values(mapping,
 *                        fe,
 *                        QGauss<2>(fe.degree + 2),
 *                        update_values | update_normal_vectors |
 *                          update_quadrature_points | update_JxW_values)
 *       , fe_values_T(mapping,
 *                     dof_hand_T.get_fe(),
 *                     QGauss<3>(fe.degree + 2),
 *                     update_values | update_gradients)
 *       , 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())
 *       , permeability_list(n_q_points)
 *       , gamma_list(n_q_points_face)
 *       , T_values(n_q_points)
 *       , dof_hand_T(dof_hand_T)
 *       , dofs_T(dofs_T)
 *       , eta_squared(eta_squared)
 *       , boundary_condition_type(boundary_condition_type)
 *     {}
 *   
 *     Solver::AssemblyScratchData::AssemblyScratchData(
 *       const AssemblyScratchData &scratch_data)
 *       : permeability()
 *       , gamma()
 *       , mapping(scratch_data.mapping.get_degree())
 *       , fe_values(mapping,
 *                   scratch_data.fe_values.get_fe(),
 *                   scratch_data.fe_values.get_quadrature(),
 *                   update_gradients | update_values | 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)
 *       , fe_values_T(mapping,
 *                     scratch_data.fe_values_T.get_fe(),
 *                     scratch_data.fe_values_T.get_quadrature(),
 *                     update_values | update_gradients)
 *       , 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())
 *       , permeability_list(n_q_points)
 *       , gamma_list(n_q_points_face)
 *       , T_values(n_q_points)
 *       , dof_hand_T(scratch_data.dof_hand_T)
 *       , dofs_T(scratch_data.dofs_T)
 *       , eta_squared(scratch_data.eta_squared)
 *       , boundary_condition_type(scratch_data.boundary_condition_type)
 *     {}
 *   
 * @endcode
 * 
 * The following function assembles a fraction of the system matrix and the
 * system right-hand side related to a single cell. These fractions are
 * `copy_data.cell_matrix` and `copy_data.cell_rhs`. They are copied into
 * to the system matrix, @f$A_{ij}@f$, and the right-hand side, @f$b_i@f$, by the
 * function `Solver::copy_local_to_global()`.
 *   
 
 * 
 * In the first three statements of the function we reinitialize the matrices
 * and vectors related to the current cell. Next, we extract the cells from
 * the cell pair. We extract them in the correct order, see above. After that
 * we compute the finite element values for both types of the finite elements,
 * compute the permeability, declare the vector values extractor, `ve`, and
 * compute current vector potential, @f$\vec{T}@f$, at quadrature points. Next,
 * we compute the three [volume integrals](@ref Step97_Numerical_Recipe_A),
 * @f$I_{a1}@f$, @f$I_{a3}@f$, and @f$I_{b3-1}@f$ in the three nested `for` loops. The
 * labels of the integrals are the same as in the introduction to this
 * tutorial. The program can be [switched](@ref Step97_TXT_BCSwitch) between
 * the Dirichlet, Neumann, and Robin boundary conditions. For this reason, we
 * compute the surface integral @f$I_{a2}@f$ only if the Robin boundary condition
 * is chosen by the user. For this reason, we use the `if` statement to make
 * sure that the surface integral @f$I_{a2}@f$ is added to the components of the
 * system matrix only if the Robin boundary condition is chosen by the user.
 * In all other cases the integral @f$I_{a2}@f$ is omitted. Omitting the @f$I_{a2}@f$
 * integral is as good as setting @f$\gamma = 0@f$ in the
 * [boundary value problem](@ref Step97_BVP_A). In the last line of the
 * function we query the dof indices on the current cell and store them in
 * the copy data structure, so we know to which locations of the system
 * matrix and right-hand side the components of the cell matrix and cell
 * right-hand side must be copied.
 * 
 * @code
 *     void Solver::system_matrix_local(const IteratorPair  &IP,
 *                                      AssemblyScratchData &scratch_data,
 *                                      AssemblyCopyData    &copy_data)
 *     {
 *       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);
 *   
 *       const typename DoFHandler<3>::active_cell_iterator cell = std::get<0>(*IP);
 *       const typename DoFHandler<3>::active_cell_iterator cell_T =
 *         std::get<1>(*IP);
 *   
 *       scratch_data.fe_values.reinit(cell);
 *       scratch_data.fe_values_T.reinit(cell_T);
 *   
 *       scratch_data.permeability.value_list(cell->material_id(),
 *                                            scratch_data.permeability_list);
 *   
 *       const FEValuesExtractors::Vector ve(0);
 *   
 *       scratch_data.fe_values_T[ve].get_function_values(scratch_data.dofs_T,
 *                                                        scratch_data.T_values);
 *   
 *       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.
 *                     (1.0 / scratch_data.permeability_list[q_index]) * // 1 / mu
 *                     (scratch_data.fe_values[ve].curl(i,
 *                                                      q_index) * // curl phi_i(x_q)
 *                        scratch_data.fe_values[ve].curl(j,
 *                                                        q_index) // curl phi_j(x_q)
 *                      +
 *                      scratch_data.eta_squared * // eta^2
 *                        scratch_data.fe_values[ve].value(i,
 *                                                         q_index) *  // phi_i(x_q)
 *                        scratch_data.fe_values[ve].value(j, q_index) // phi_j(x_q)
 *                      ) *
 *                     scratch_data.fe_values.JxW(q_index); // dx
 *                 }
 *               copy_data.cell_rhs(i) += // Integral I_b3-1.
 *                 (scratch_data.T_values[q_index] *
 *                  scratch_data.fe_values[ve].curl(i, q_index)) *
 *                 scratch_data.fe_values.JxW(q_index); // T(x_q).(curl phi_i(x_q))dx
 *             }
 *         }
 *   
 *       if (scratch_data.boundary_condition_type == BoundaryConditionType::Robin)
 *         {
 *           for (unsigned int f = 0; f < cell->n_faces(); ++f)
 *             {
 *               if (cell->face(f)->at_boundary())
 *                 {
 *                   scratch_data.fe_face_values.reinit(cell, f);
 *   
 *                   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)
 *                         {
 *                           scratch_data.gamma.value_list(
 *                             scratch_data.fe_face_values.get_quadrature_points(),
 *                             scratch_data.gamma_list);
 *   
 *                           for (unsigned int j = 0; j < scratch_data.dofs_per_cell;
 *                                ++j)
 *                             {
 *                               copy_data.cell_matrix(i, j) += // Integral I_a2.
 *                                 scratch_data.gamma_list[q_index_face] * // gamma
 *                                 (cross_product_3d(scratch_data.fe_face_values
 *                                                     .normal_vector(q_index_face),
 *                                                   scratch_data.fe_face_values[ve]
 *                                                     .value(i, q_index_face)) *
 *                                  cross_product_3d(scratch_data.fe_face_values
 *                                                     .normal_vector(q_index_face),
 *                                                   scratch_data.fe_face_values[ve]
 *                                                     .value(j, q_index_face)))
 *                                 /* (n x phi_i(x_q)).(n x phi_j(x_q)) */
 *                                 * scratch_data.fe_face_values.JxW(
 *                                     q_index_face); // dx
 *                             }                      // for j
 *                         }                          // for i
 *                     }                              // for q_index_face
 *                 } // if (cell->face(f)->at_boundary())
 *             }     // for f
 *         }         /* if (scratch_data.boundary_condition_type ==
 *                      BoundaryConditionType::Robin) */
 *   
 *       cell->get_dof_indices(copy_data.local_dof_indices);
 *     }
 *   
 * @endcode
 * 
 * This function copies the components of a cell matrix and a cell right-hand
 * side into the system matrix, @f$A_{ij}@f$, and the system right-hand side,
 * @f$b_i@f$.
 * 
 * @code
 *     void Solver::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);
 *     }
 *   
 * @endcode
 * 
 * The following function solves the system of linear equations. In theory,
 * the CG solver can solve an @f$m \times m@f$ system of linear equations in at
 * most @f$m@f$ steps. Accordingly, we set the maximum number of iteration steps
 * to `system_rhs.size()`. The stopping condition is
 * \f[\|\boldsymbol{b} - \boldsymbol{A}\boldsymbol{c}\| < 10^{-6}
 * \|\boldsymbol{b}\|.\f] As soon as we use constraints, we must not forget to
 * distribute them.
 * 
 * @code
 *     void Solver::solve()
 *     {
 *       SolverControl control(system_rhs.size(), 1.0e-6 * system_rhs.l2_norm());
 *   
 *       GrowingVectorMemory<Vector<double>> memory;
 *       SolverCG<Vector<double>>            cg(control, memory);
 *   
 *       PreconditionSSOR<SparseMatrix<double>> preconditioner;
 *       preconditioner.initialize(system_matrix, 1.2);
 *   
 *       cg.solve(system_matrix, solution, system_rhs, preconditioner);
 *   
 *       constraints.distribute(solution);
 *     }
 *   
 * @endcode
 * 
 * The following function saves the computed magnetic vector potential into a
 * `.vtu` file.
 * 
 * @code
 *     void Solver::output_results() const
 *     {
 *       std::vector<std::string> solution_names(3, "VectorField");
 *       std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *         interpretation(3,
 *                        DataComponentInterpretation::component_is_part_of_vector);
 *   
 *       DataOut<3> data_out;
 *   
 *       data_out.add_data_vector(dof_handler,
 *                                solution,
 *                                solution_names,
 *                                interpretation);
 *   
 *       DataOutBase::VtkFlags flags;
 *       flags.write_higher_order_cells = true;
 *       data_out.set_flags(flags);
 *   
 *       const MappingQ<3> mapping(mapping_degree);
 *   
 *       data_out.build_patches(mapping,
 *                              fe.degree + 2,
 *                              DataOut<3>::CurvedCellRegion::curved_inner_cells);
 *   
 *       std::ofstream ofs(fname + ".vtu");
 *       data_out.write_vtu(ofs);
 *     }
 *   
 * @endcode
 * 
 * The `run` function below aggregates all the computation steps in the
 * correct order.
 * 
 * @code
 *     void Solver::run()
 *     {
 *       setup();
 *       assemble();
 *       solve();
 *       output_results();
 *       clear();
 *     }
 *   } // namespace SolverA
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_97-ProjectorfromHcurltoHdiv"></a> 
 * <h3>Projector from H(curl) to H(div)</h3>
 * The following namespace contains all the code related to the conversion of
 * the magnetic vector potential, @f$\vec{A}@f$, into magnetic field, @f$\vec{B}@f$.
 * The magnetic vector potential is modeled by the FE_Nedelec finite elements,
 * while the magnetic field is modeled by the FE_RaviartThomas finite elements.
 * This code is also used for converting the current vector potential,
 * @f$\vec{T}@f$, into the free-current density, @f$\vec{J}_f@f$.
 * 
 * @code
 *   namespace ProjectorHcurlToHdiv
 *   {
 *   
 * @endcode
 * 
 * This class implements the solver that minimizes the
 * [functional](@ref Step97_Functional_B) @f$F(\vec{B})@f$ or @f$F(\vec{J}_f)@f$, see
 * the introduction. The input vector field, @f$\vec{A}@f$ or @f$\vec{T}@f$, is fed to
 * the solver by means of the input parameters `dof_handler_Hcurl` and
 * `solution_Hcurl`. Moreover, this solver reuses the mesh on which the input
 * vector field has been computed. The reference to the mesh is passed via the
 * input parameter `triangulation_Hcurl`. There are no constraints this time
 * around as we are not going to apply the Dirichlet boundary condition.
 * 
 * @code
 *     class Solver
 *     {
 *     public:
 *       Solver() = delete;
 *       Solver(const unsigned int p, // Degree of the Raviart-Thomas finite elements
 *              const unsigned int mapping_degree,
 *              const Triangulation<3> &triangulation_Hcurl,
 *              const DoFHandler<3>    &dof_handler_Hcurl,
 *              const Vector<double>   &solution_Hcurl,
 *              const std::string      &fname          = "data",
 *              const Function<3>      *exact_solution = nullptr);
 *   
 *       double get_L2_norm()
 *       {
 *         return L2_norm;
 *       };
 *   
 *       unsigned int get_n_cells() const
 *       {
 *         return triangulation_Hcurl.n_active_cells();
 *       }
 *   
 *       types::global_dof_index get_n_dofs() const
 *       {
 *         return dof_handler_Hdiv.n_dofs();
 *       }
 *   
 *       void setup();                // Initializes dofs, vectors, matrices.
 *       void assemble();             // Assembles the system of linear equations.
 *       void solve();                // Solves the system of linear equations.
 *       void output_results() const; // Saves computed Jf or B into a vtu file.
 *       void compute_error_norms();  // Computes L^2 error norm.
 *       void project_exact_solution_fcn(); // Projects exact solution.
 *       void clear()
 *       {
 *         system_matrix.clear();
 *         system_rhs.reinit(0);
 *       }
 *       void run(); /* Executes the last seven functions in the proper order
 *                      and measures the execution time for each function. */
 *   
 *     private:
 *       const Triangulation<3> &triangulation_Hcurl;
 *       const DoFHandler<3>    &dof_handler_Hcurl;
 *       const Vector<double>   &solution_Hcurl;
 *   
 * @endcode
 * 
 * The following data members are typical for all deal.II simulations:
 * triangulation, finite elements, dof handlers, etc. The constraints
 * are used to enforce the Dirichlet boundary conditions. The names of the
 * data members are self-explanatory.
 * 
 * @code
 *       const FE_RaviartThomas<3> fe_Hdiv;
 *       DoFHandler<3>             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<3> *exact_solution;
 *   
 *       const unsigned int mapping_degree;
 *   
 *       Vector<double> L2_per_cell;
 *       double         L2_norm;
 *   
 *       const std::string fname;
 *   
 * @endcode
 * 
 * This time we have two dof handlers, `dof_handler_Hcurl` for the input
 * vector field and `dof_handler_Hdiv` for the output vector field. The
 * WorkStream needs to walk through the two dof handlers synchronously.
 * For this purpose we will pair two active cells iterators (one from
 * `dof_handler_Hcurl`, another from `dof_handler_Hdiv`) to be walked
 * through synchronously. For that we need the `IteratorPair` type.
 * 
 * @code
 *       using IteratorTuple =
 *         std::tuple<typename DoFHandler<3>::active_cell_iterator,
 *                    typename DoFHandler<3>::active_cell_iterator>;
 *   
 *       using IteratorPair = SynchronousIterators<IteratorTuple>;
 *   
 * @endcode
 * 
 * The program utilizes the WorkStream technology. The @ref step_9 "step-9" tutorial
 * does a much better job of explaining the workings of WarkStream.
 * Reading the @ref workstream_paper "WorkStream paper" is recommended.
 * The following structures and functions are related to WorkStream.
 * 
 * @code
 *       struct AssemblyScratchData
 *       {
 *         AssemblyScratchData(const FiniteElement<3> &fe,
 *                             const DoFHandler<3>    &dof_handr_Hcurl,
 *                             const Vector<double>   &dofs_Hcurl,
 *                             const unsigned int      mapping_degree);
 *   
 *         AssemblyScratchData(const AssemblyScratchData &scratch_data);
 *   
 *         MappingQ<3> mapping;
 *         FEValues<3> fe_values_Hdiv;
 *         FEValues<3> fe_values_Hcurl;
 *   
 *         const unsigned int dofs_per_cell;
 *         const unsigned int n_q_points;
 *   
 *         std::vector<Tensor<1, 3>> curl_vec_in_Hcurl;
 *   
 *         const DoFHandler<3>  &dof_hand_Hcurl;
 *         const Vector<double> &dofs_Hcurl;
 *       };
 *   
 *       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);
 *     };
 *   
 *     Solver::Solver(const unsigned int      p,
 *                    const unsigned int      mapping_degree,
 *                    const Triangulation<3> &triangulation_Hcurl,
 *                    const DoFHandler<3>    &dof_handler_Hcurl,
 *                    const Vector<double>   &solution_Hcurl,
 *                    const std::string      &fname,
 *                    const Function<3>      *exact_solution)
 *       : triangulation_Hcurl(triangulation_Hcurl)
 *       , dof_handler_Hcurl(dof_handler_Hcurl)
 *       , solution_Hcurl(solution_Hcurl)
 *       , fe_Hdiv(p)
 *       , exact_solution(exact_solution)
 *       , mapping_degree(mapping_degree)
 *       , fname(fname)
 *     {
 *       Assert(exact_solution != nullptr,
 *              ExcMessage("The exact solution is missing."));
 *     }
 *   
 * @endcode
 * 
 * The following function initializes the dofs, vectors and matrices. This
 * time there are no constraints as we do not apply Dirichlet boundary
 * condition.
 * 
 * @code
 *     void Solver::setup()
 *     {
 *       constraints.close();
 *   
 *       dof_handler_Hdiv.reinit(triangulation_Hcurl);
 *       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 (Settings::project_exact_solution && exact_solution)
 *         projected_exact_solution.reinit(dof_handler_Hdiv.n_dofs());
 *   
 *       if (exact_solution)
 *         L2_per_cell.reinit(triangulation_Hcurl.n_active_cells());
 *     }
 *   
 * @endcode
 * 
 * Formally, the following function assembles the system of linear equations.
 * In reality, however, it just spells all the magic words to get the
 * WorkStream going. The interesting part, i.e., the actual assembling of the
 * system matrix and the right-hand side happens below in the
 * `Solver::system_matrix_local` function. Note that this time the first two
 * input parameters to `WorkStream::run` are pairs of iterators, not
 * iterators themselves as per usual. Note also the order in which we package
 * the iterators: first the iterator of `dof_handler_Hdiv`, then the iterator
 * of the `dof_handler_Hcurl`. We will extract them in the same order.
 * 
 * @code
 *     void Solver::assemble()
 *     {
 *       WorkStream::run(IteratorPair({dof_handler_Hdiv.begin_active(),
 *                                     dof_handler_Hcurl.begin_active()}),
 *                       IteratorPair(
 *                         {dof_handler_Hdiv.end(), dof_handler_Hcurl.end()}),
 *                       *this,
 *                       &Solver::system_matrix_local,
 *                       &Solver::copy_local_to_global,
 *                       AssemblyScratchData(fe_Hdiv,
 *                                           dof_handler_Hcurl,
 *                                           solution_Hcurl,
 *                                           mapping_degree),
 *                       AssemblyCopyData());
 *     }
 *   
 * @endcode
 * 
 * The following two constructors initialize scratch data from the input
 * parameters and from another object of the same type, i.e., a copy
 * constructor.
 * 
 * @code
 *     Solver::AssemblyScratchData::AssemblyScratchData(
 *       const FiniteElement<3> &fe,
 *       const DoFHandler<3>    &dof_hand_Hcurl,
 *       const Vector<double>   &dofs_Hcurl,
 *       const unsigned int      mapping_degree)
 *       : mapping(mapping_degree)
 *       , fe_values_Hdiv(mapping,
 *                        fe,
 *                        QGauss<3>(fe.degree + 2),
 *                        update_values | update_JxW_values)
 *       , fe_values_Hcurl(mapping,
 *                         dof_hand_Hcurl.get_fe(),
 *                         QGauss<3>(fe.degree + 2),
 *                         update_gradients)
 *       , dofs_per_cell(fe_values_Hdiv.dofs_per_cell)
 *       , n_q_points(fe_values_Hdiv.get_quadrature().size())
 *       , curl_vec_in_Hcurl(n_q_points)
 *       , dof_hand_Hcurl(dof_hand_Hcurl)
 *       , dofs_Hcurl(dofs_Hcurl)
 *     {}
 *   
 *     Solver::AssemblyScratchData::AssemblyScratchData(
 *       const AssemblyScratchData &scratch_data)
 *       : mapping(scratch_data.mapping.get_degree())
 *       , fe_values_Hdiv(mapping,
 *                        scratch_data.fe_values_Hdiv.get_fe(),
 *                        scratch_data.fe_values_Hdiv.get_quadrature(),
 *                        update_values | update_JxW_values)
 *       , fe_values_Hcurl(mapping,
 *                         scratch_data.fe_values_Hcurl.get_fe(),
 *                         scratch_data.fe_values_Hcurl.get_quadrature(),
 *                         update_gradients)
 *       , dofs_per_cell(fe_values_Hdiv.dofs_per_cell)
 *       , n_q_points(fe_values_Hdiv.get_quadrature().size())
 *       , curl_vec_in_Hcurl(scratch_data.n_q_points)
 *       , dof_hand_Hcurl(scratch_data.dof_hand_Hcurl)
 *       , dofs_Hcurl(scratch_data.dofs_Hcurl)
 *     {}
 *   
 * @endcode
 * 
 * The following function assembles a fraction of the system matrix and the
 * system right-hand side related to a single cell. These fractions are
 * `copy_data.cell_matrix` and `copy_data.cell_rhs`. They are copied into
 * to the system matrix, @f$A_{ij}@f$, and the right-hand side, @f$b_i@f$, by the
 * function `Solver::copy_local_to_global()`.
 *   
 
 * 
 * First, we reinitialize the matrices and vectors related to the current
 * cell, update the finite element values, and compute the curl of the input
 * vector field at quadrature points. The variable
 * `scratch_data.curl_vec_in_Hcurl` denotes the curl of the input vector
 * field, @f$\vec{\nabla} \times \vec{T}@f$ or @f$\vec{\nabla} \times \vec{A}@f$,
 * depending on the context. Second, we compute the
 * [components](@ref Step97_Numerical_Recipe_B) of the cell matrix and cell
 * right-hand side in the three nested `for` loops. The labels of the
 * integrals are the same as in the introduction to this tutorial. Third, we
 * query the dof indices on the current cell and store them in the copy data
 * structure, so we know to which locations of the system matrix and
 * right-hand side the components of the cell matrix and cell right-hand side
 * must be copied.
 * 
 * @code
 *     void Solver::system_matrix_local(const IteratorPair  &IP,
 *                                      AssemblyScratchData &scratch_data,
 *                                      AssemblyCopyData    &copy_data)
 *     {
 *       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);
 *   
 *       const typename DoFHandler<3>::active_cell_iterator cell_Hdiv =
 *         std::get<0>(*IP);
 *       const typename DoFHandler<3>::active_cell_iterator cell_Hcurl =
 *         std::get<1>(*IP);
 *   
 *       scratch_data.fe_values_Hdiv.reinit(cell_Hdiv);
 *       scratch_data.fe_values_Hcurl.reinit(cell_Hcurl);
 *   
 *       const FEValuesExtractors::Vector ve(0);
 *   
 *       scratch_data.fe_values_Hcurl[ve].get_function_curls(
 *         scratch_data.dofs_Hcurl, scratch_data.curl_vec_in_Hcurl);
 *   
 *       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_a
 *                     scratch_data.fe_values_Hdiv[ve].value(i,
 *                                                           q_index) * // phi_i(x_q)
 *                     scratch_data.fe_values_Hdiv[ve].value(j,
 *                                                           q_index) * // phi_j(x_q)
 *                     scratch_data.fe_values_Hdiv.JxW(q_index);        // dx
 *                 }
 *   
 *               copy_data.cell_rhs(i) += // Integral I_b
 *                 scratch_data
 *                   .curl_vec_in_Hcurl[q_index] * /* curl A(x_q) OR curl T(x_q),
 *                                                    depending on the context. */
 *                 scratch_data.fe_values_Hdiv[ve].value(i, q_index) * // phi_i(x_q)
 *                 scratch_data.fe_values_Hdiv.JxW(q_index);           // dx
 *             }
 *         }
 *   
 *       cell_Hdiv->get_dof_indices(copy_data.local_dof_indices);
 *     }
 *   
 * @endcode
 * 
 * The following function copies the components of a cell matrix and a cell
 * right-hand side into the system matrix, @f$A_{ij}@f$, and the system right-hand
 * side, @f$b_i@f$.
 * 
 * @code
 *     void Solver::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);
 *     }
 *   
 * @endcode
 * 
 * The following two functions compute the error norms and project the exact
 * solution.
 * 
 * @code
 *     void Solver::compute_error_norms()
 *     {
 *       if (exact_solution)
 *         {
 *           const Weight               weight;
 *           const Function<3, double> *mask = &weight;
 *   
 *           VectorTools::integrate_difference(MappingQ<3>(mapping_degree),
 *                                             dof_handler_Hdiv,
 *                                             solution_Hdiv,
 *                                             *exact_solution,
 *                                             L2_per_cell,
 *                                             QGauss<3>(fe_Hdiv.degree + 4),
 *                                             VectorTools::L2_norm,
 *                                             mask);
 *   
 *           L2_norm = VectorTools::compute_global_error(triangulation_Hcurl,
 *                                                       L2_per_cell,
 *                                                       VectorTools::L2_norm);
 *         }
 *     }
 *   
 *     void Solver::project_exact_solution_fcn()
 *     {
 *       if (Settings::project_exact_solution && exact_solution)
 *         {
 *           AffineConstraints<double> constraints_empty;
 *   
 *           constraints_empty.clear();
 *   
 *           DoFTools::make_hanging_node_constraints(dof_handler_Hdiv,
 *                                                   constraints_empty);
 *   
 *           constraints_empty.close();
 *   
 *           VectorTools::project(MappingQ<3>(mapping_degree),
 *                                dof_handler_Hdiv,
 *                                constraints_empty,
 *                                QGauss<3>(fe_Hdiv.degree + 2),
 *                                *exact_solution,
 *                                projected_exact_solution);
 *         }
 *     }
 *   
 * @endcode
 * 
 * The following function solves the system of linear equations. In theory,
 * the CG solver can solve an @f$m \times m@f$ system of linear equations in at
 * most @f$m@f$ steps. Accordingly, we set the maximum number of iteration steps
 * to `system_rhs.size()`. The stopping condition is
 * \f[\|\boldsymbol{b} - \boldsymbol{A}\boldsymbol{c}\| < 10^{-6}
 * \|\boldsymbol{b}\|.\f] This time the constraints are empty as we do not
 * use the Dirichlet boundary condition. Consequently, we do not have to
 * distribute the constraints.
 * 
 * @code
 *     void Solver::solve()
 *     {
 *       SolverControl control(system_rhs.size(), 1.0e-6 * system_rhs.l2_norm());
 *   
 *       GrowingVectorMemory<Vector<double>> memory;
 *       SolverCG<Vector<double>>            cg(control, memory);
 *   
 *       PreconditionSSOR<SparseMatrix<double>> preconditioner;
 *       preconditioner.initialize(system_matrix, 1.2);
 *   
 *       cg.solve(system_matrix, solution_Hdiv, system_rhs, preconditioner);
 *     }
 *   
 * @endcode
 * 
 * The following function saves the computed fields into a `.vtu` file.
 * This time we also save the projected exact solution and the @f$L^2@f$ error
 * norm. The exact solution is only saved if
 * `Settings::project_exact_solution = true`
 * 
 * @code
 *     void Solver::output_results() const
 *     {
 *       std::vector<std::string> solution_names(3, "VectorField");
 *       std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *         interpretation(3,
 *                        DataComponentInterpretation::component_is_part_of_vector);
 *   
 *       DataOut<3> data_out;
 *   
 *       data_out.add_data_vector(dof_handler_Hdiv,
 *                                solution_Hdiv,
 *                                solution_names,
 *                                interpretation);
 *   
 *       if (Settings::project_exact_solution)
 *         {
 *           std::vector<std::string> solution_names_ex(3, "VectorFieldExact");
 *   
 *           data_out.add_data_vector(dof_handler_Hdiv,
 *                                    projected_exact_solution,
 *                                    solution_names_ex,
 *                                    interpretation);
 *         }
 *   
 *       if (exact_solution)
 *         {
 *           data_out.add_data_vector(L2_per_cell, "L2norm");
 *         }
 *   
 *       DataOutBase::VtkFlags flags;
 *       flags.write_higher_order_cells = true;
 *       data_out.set_flags(flags);
 *   
 *       const MappingQ<3> mapping(mapping_degree);
 *   
 *       data_out.build_patches(mapping,
 *                              fe_Hdiv.degree + 2,
 *                              DataOut<3>::CurvedCellRegion::curved_inner_cells);
 *   
 *       std::ofstream ofs(fname + ".vtu");
 *       data_out.write_vtu(ofs);
 *     }
 *   
 * @endcode
 * 
 * The `run` function below aggregates all the computation steps in the
 * correct order.
 * 
 * @code
 *     void Solver::run()
 *     {
 *       setup();
 *       assemble();
 *       solve();
 *       compute_error_norms();
 *       project_exact_solution_fcn();
 *       output_results();
 *       clear();
 *     }
 *   } // namespace ProjectorHcurlToHdiv
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_97-Themainloop"></a> 
 * <h3>The main loop</h3>
 * 
 
 * 
 * The following class contains the main loop of the program.
 * 
 * @code
 *   class MagneticProblem
 *   {
 *   public:
 *     void run()
 *     {
 *       if (Settings::n_threads_max)
 *         MultithreadInfo::set_thread_limit(Settings::n_threads_max);
 *   
 *       MainOutputTable table_Jf(3);
 *       MainOutputTable table_B(3);
 *   
 *       std::cout << "Solving for (p = " << Settings::fe_degree
 *                 << "): " << std::flush;
 *   
 *       for (unsigned int r = 6; r < 10; r++) // Mesh refinement parameter.
 *         {
 *           table_Jf.add_value("r", r);
 *           table_Jf.add_value("p", Settings::fe_degree);
 *   
 *           table_B.add_value("r", r);
 *           table_B.add_value("p", Settings::fe_degree);
 *   
 * @endcode
 * 
 * Stage 1. Computing @f$\vec{T}@f$.
 * 
 
 * 
 * 
 * @code
 *           std::cout << "T " << std::flush;
 *   
 *           SolverT::Solver T(Settings::fe_degree,
 *                             r,
 *                             Settings::mapping_degree,
 *                             Settings::eta_squared_T,
 *                             "T_p" + std::to_string(Settings::fe_degree) + "_r" +
 *                               std::to_string(r));
 *   
 *           T.run();
 *   
 * @endcode
 * 
 * Stage 2. Computing @f$\vec{J}_f@f$.
 * 
 
 * 
 * 
 * @code
 *           std::cout << "Jf " << std::flush;
 *   
 *           ExactSolutions::FreeCurrentDensity Jf_exact;
 *   
 *           ProjectorHcurlToHdiv::Solver Jf(Settings::fe_degree,
 *                                           Settings::mapping_degree,
 *                                           T.get_tria(),
 *                                           T.get_dof_handler(),
 *                                           T.get_solution(),
 *                                           "Jf_p" +
 *                                             std::to_string(Settings::fe_degree) +
 *                                             "_r" + std::to_string(r),
 *                                           &Jf_exact);
 *   
 *           Jf.run();
 *   
 *           table_Jf.add_value("ndofs", Jf.get_n_dofs());
 *           table_Jf.add_value("ncells", Jf.get_n_cells());
 *           table_Jf.add_value("L2", Jf.get_L2_norm());
 *   
 * @endcode
 * 
 * Stage 3. Computing @f$\vec{A}@f$.
 * 
 
 * 
 * 
 * @code
 *           std::cout << "A " << std::flush;
 *   
 *           SolverA::Solver A(Settings::fe_degree,
 *                             Settings::mapping_degree,
 *                             T.get_tria(),
 *                             T.get_dof_handler(),
 *                             T.get_solution(),
 *                             Settings::eta_squared_A,
 *                             "A_p" + std::to_string(Settings::fe_degree) + "_r" +
 *                               std::to_string(r));
 *   
 *           A.run();
 *   
 * @endcode
 * 
 * Stage 4. Computing @f$\vec{B}@f$.
 * 
 
 * 
 * 
 * @code
 *           std::cout << "B " << std::flush;
 *   
 *           ExactSolutions::MagneticField B_exact;
 *   
 *           ProjectorHcurlToHdiv::Solver B(Settings::fe_degree,
 *                                          Settings::mapping_degree,
 *                                          T.get_tria(),
 *                                          A.get_dof_handler(),
 *                                          A.get_solution(),
 *                                          "B_p" +
 *                                            std::to_string(Settings::fe_degree) +
 *                                            "_r" + std::to_string(r),
 *                                          &B_exact);
 *           B.run();
 *   
 *           table_B.add_value("ndofs", B.get_n_dofs());
 *           table_B.add_value("ncells", B.get_n_cells());
 *           table_B.add_value("L2", B.get_L2_norm());
 * @endcode
 * 
 * End stage 4.
 * 
 * @code
 *         }
 *   
 *       table_Jf.save("table_Jf_p" + std::to_string(Settings::fe_degree));
 *       table_B.save("table_B_p" + std::to_string(Settings::fe_degree));
 *       std::cout << std::endl;
 *     }
 *   };
 *   
 *   int main()
 *   {
 *     try
 *       {
 *         MagneticProblem problem;
 *         problem.run();
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *   
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *   
 *     return 0;
 *   }
 * @endcode
<a name="step_97-Results"></a><h1>Results</h1>
 
 
The program generates the following output in the command line interface by
default.
 
@code
Solving for (p = 0): T Jf A B T Jf A B T Jf A B T Jf A B
@endcode
 
The program assumes the finite elements of the lowermost degree, @f$p = 0@f$.
To change the degree of the finite elements, say @f$p = 2@f$, one needs to change
the setting `Settings::fe_degree = 2` and rebuild the program.
 
The program also dumps a number of files in the current directory. In the default
configuration these files are:
- `.vtu` files. They contain the computed vector fields. Recall that the spherical
  manifold is attached to many cell faces. Consequently, these cell faces are
  curved. They look more like patches of a sphere. Furthermore, the shape
  functions are mapped from the reference cell to the real mesh cells by the
  second-order mapping to accommodate the cells with curved faces. For these
  reasons, one needs to use a visualization software that can deal with curved
  faces and the higher-order mapping. A fresh version of ParaView is recommended.
  Visit will not do. The <a href="https://github.com/dealii/dealii/wiki/Notes-on-visualizing-high-order-output">
  Notes on visualizing high order output</a> provide more information on this topic.
- `.tex` files. These files contain the convergence tables.
 
The following provides examples of the convergence tables simulated with the
default settings for three different degrees of the finite elements,
@f$p = 0, 1, 2@f$.
 
<table>
<caption>Convergence table @f$\vec{J}_f@f$.</caption>
    <tr>
        <th>p</td>
        <th>r</td>
        <th>cells</td>
        <th>dofs</td>
        <th>@f$\|e\|_{L^2}@f$</td>
        <th>@f$\alpha_{L^2}@f$</td>
    </tr>
    <tr>
        <td>0</td>
        <td>6</td>
        <td>4625</td>
        <td>13950</td>
        <td>1.66e-01</td>
        <td>-</td>
    </tr>
    <tr>
        <td>0</td>
        <td>7</td>
        <td>7992</td>
        <td>24084</td>
        <td>1.38e-01</td>
        <td>0.99</td>
    </tr>
    <tr>
        <td>0</td>
        <td>8</td>
        <td>12691</td>
        <td>38220</td>
        <td>1.19e-01</td>
        <td>0.99</td>
    </tr>
    <tr>
        <td>0</td>
        <td>9</td>
        <td>18944</td>
        <td>57024</td>
        <td>1.04e-01</td>
        <td>0.99</td>
    </tr>
    <tr>
        <td>1</td>
        <td>6</td>
        <td>4625</td>
        <td>111300</td>
        <td>8.12e-04</td>
        <td>-</td>
    </tr>
    <tr>
        <td>1</td>
        <td>7</td>
        <td>7992</td>
        <td>192240</td>
        <td>4.97e-04</td>
        <td>2.69</td>
    </tr>
    <tr>
        <td>1</td>
        <td>8</td>
        <td>12691</td>
        <td>305172</td>
        <td>3.32e-04</td>
        <td>2.61</td>
    </tr>
    <tr>
        <td>1</td>
        <td>9</td>
        <td>18944</td>
        <td>455424</td>
        <td>2.37e-04</td>
        <td>2.54</td>
    </tr>
    <tr>
        <td>2</td>
        <td>6</td>
        <td>4625</td>
        <td>375300</td>
        <td>6.78e-04</td>
        <td>-</td>
    </tr>
    <tr>
        <td>2</td>
        <td>7</td>
        <td>7992</td>
        <td>648324</td>
        <td>3.94e-04</td>
        <td>2.97</td>
    </tr>
    <tr>
        <td>2</td>
        <td>8</td>
        <td>12691</td>
        <td>1029294</td>
        <td>2.49e-04</td>
        <td>2.98</td>
    </tr>
    <tr>
        <td>2</td>
        <td>9</td>
        <td>18944</td>
        <td>1536192</td>
        <td>1.67e-04</td>
        <td>2.99</td>
    </tr>
</table>
<br>
<table>
<caption>Convergence table @f$\vec{B}@f$.</caption>
    <tr>
        <th>p</td>
        <th>r</td>
        <th>cells</td>
        <th>dofs</td>
        <th>@f$\|e\|_{L^2}@f$</td>
        <th>@f$\alpha_{L^2}@f$</td>
    </tr>
    <tr>
        <td>0</td>
        <td>6</td>
        <td>4625</td>
        <td>13950</td>
        <td>8.84e-08</td>
        <td>-</td>
    </tr>
    <tr>
        <td>0</td>
        <td>7</td>
        <td>7992</td>
        <td>24084</td>
        <td>7.36e-08</td>
        <td>1.00</td>
    </tr>
    <tr>
        <td>0</td>
        <td>8</td>
        <td>12691</td>
        <td>38220</td>
        <td>6.30e-08</td>
        <td>1.01</td>
    </tr>
    <tr>
        <td>0</td>
        <td>9</td>
        <td>18944</td>
        <td>57024</td>
        <td>5.51e-08</td>
        <td>1.00</td>
    </tr>
    <tr>
        <td>1</td>
        <td>6</td>
        <td>4625</td>
        <td>111300</td>
        <td>4.41e-09</td>
        <td>-</td>
    </tr>
    <tr>
        <td>1</td>
        <td>7</td>
        <td>7992</td>
        <td>192240</td>
        <td>3.11e-09</td>
        <td>1.91</td>
    </tr>
    <tr>
        <td>1</td>
        <td>8</td>
        <td>12691</td>
        <td>305172</td>
        <td>2.23e-09</td>
        <td>2.18</td>
    </tr>
    <tr>
        <td>1</td>
        <td>9</td>
        <td>18944</td>
        <td>455424</td>
        <td>1.71e-09</td>
        <td>1.96</td>
    </tr>
    <tr>
        <td>2</td>
        <td>6</td>
        <td>4625</td>
        <td>375300</td>
        <td>1.84e-10</td>
        <td>-</td>
    </tr>
    <tr>
        <td>2</td>
        <td>7</td>
        <td>7992</td>
        <td>648324</td>
        <td>1.03e-10</td>
        <td>3.21</td>
    </tr>
    <tr>
        <td>2</td>
        <td>8</td>
        <td>12691</td>
        <td>1029294</td>
        <td>6.08e-11</td>
        <td>3.40</td>
    </tr>
    <tr>
        <td>2</td>
        <td>9</td>
        <td>18944</td>
        <td>1536192</td>
        <td>4.04e-11</td>
        <td>3.07</td>
    </tr>
</table>
 
The following notations were used in the headers of the tables:
 
- p - the degree of the finite elements.
 
- r - the mesh refinement parameter, i.e., the number of nodes on the transfinite
lines.
 
- cells - the total amount of active cells.
 
- dofs - the amount of degrees of freedom.
 
-@f$\|e\|_{L^2}@f$ - the @f$L^2@f$ error norm.
 
-@f$\alpha_{L^2}@f$ - the order of convergence of the @f$L^2@f$ error norm.
 
The vector representations of the calculated vector fields, @f$\vec{J}_f@f$ and
@f$\vec{B}@f$, are illustrated above by the first figure on this page. The figures
below illustrate slices of the magnitudes of these fields. The figures below
were simulated with @f$p = 2@f$ and @f$r = 9@f$. Visual inspection of the vector
potentials is not very informative as their conservative portions are unknown.
 
@htmlonly
<p align="center">
  <img src="https://dealii.org/images/steps/developer/step-97-Jf.svg" alt="The result - free-current
 density" height="531">
</p>
@endhtmlonly
 
@htmlonly
<p align="center">
  <img src="https://dealii.org/images/steps/developer/step-97-B.svg" alt="The result - magnetic
  field" height="531">
</p>
@endhtmlonly
 
@anchor Step97_PossibilitiesForExtensions
<a name="step_97-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
Repeat the simulations for the three types of the boundary conditions,
Dirichlet, Neumann, and Robin. The Robin boundary condition is supposed to be
superior to the other two. Look at the simulated data to see that this is indeed
the case. You can save the projected exact solution next to the simulated
solutions into the `.vtu` files, just set `Settings::project_exact_solution = true`.
ParaView has "Plot Over Line" filter. You can use this filter to visualize the
difference between the exact solution and a solution simulated with a particular
boundary condition. You can also draw conclusions by observing the convergence
tables. Keep in mind the @f$\eta^2@f$ parameter. Increase it if the CG solver chokes
while you are experimenting. Note that the benefits offered by the Robin
boundary condition are observed the best when higher-order finite elements are
used, i.e., @f$p = 1@f$ and @f$p = 2@f$.
 
The Robin boundary condition as described above is also called the first-order
asymptotic boundary condition (ABC). There exist ABCs of higher orders
@cite gratkowski2010p. Implement and test the second-order ABC to see if it
performs any better. There exist improvised asymptotic boundary conditions, IABCs,
@cite meeker2013a. Try to implement the first order IABC.
 *
 *
<a name="step_97-PlainProg"></a>
<h1> The plain program</h1>
@include "step-97.cc"
*/