This tutorial depends on step-9.

Table of contents

Introduction
The commented program

Results
- Possibilities for extensions
The plain program

This program was contributed by Siarhei Uzunbajakau, CEM Books, 2025.

Introduction

Problem definition

This tutorial illustrates application of the FE_Nedelec and FE_RaviartThomas finite elements to problems in electromagnetics. To this end, we will compute the magnetic field induced by a coil with a magnetic core and compare the result to an exact closed-form analytical expression of the magnetic field. We will repeat the comparison on a set of progressively refined meshes and observe that the $L^2$ error norm decreases with the size of mesh cells. We will also vary the degree of the finite elements and observe that the higher-degree finite elements imply faster convergence of the $L^2$ error norms. Experimenting with boundary conditions is suggested as a possibility for extension.

The first figure below illustrates a cross section of the coil. The second figure illustrates the corresponding problem domain, coordinate system, and the relevant notations. The coil consists of two spherical shells. The inner shell, depicted in green, is the magnetic core. The walls of the core are filled with a soft magnetic material of permeability $\mu_1$ . The magnetic material is assumed to be homogeneous, linear, lossless, and isotropic. The permeability of the rest of the space (the ball-shaped space inside the core and the space outside the core) equals to that of the free space, $\mu_0$ . The permeability of the entire space can be expressed as

$\begin{equation}\mu = \left\{ \begin{aligned} & \mu_0 &\text{if } \text{ } & r < a_1\\ & \mu_1 &\text{if } \text{ } & a_1 \le r \le b_1\\ & \mu_0 &\text{if } \text{ } & r > b_1. \end{aligned} \right. \end{equation}$

The outer shell, depicted in blue, contains the current-carrying windings of the coil. The windings are modeled by a prescribed free-current density. The prescribed free-current density fills all the space available in the wall of the shell. In the figure below, however, only three cross sections of the free-current density are shown. The free-current density is the highest on the equator and equals zero at the two poles. It can be described as

$\‍[\vec{J}_f = \left\{ \begin{aligned} & 0 &&\text{if }&& r < a_2\\ & K_0 (-y\hat{i} + x\hat{j})&&\text{if }&& a_2 \le r \le b_2\\ & 0 &&\text{if }&& r > b_2\\ \end{aligned} \right. \‍]$

in the Cartesian coordinate system.

The coil

There is no consensus in the literature on what to call a magnetic field. On this page we adopt the nomenclature from [griffiths2017b]. Consequently, we call a magnetic field the vector quantity $\vec{B}$ such that the Ampere's law in free space in static approximation reads

$\‍[\begin{equation} \vec{\nabla} \times\vec{B} = \mu_0\vec{J}, \end{equation} \‍]$

where $\vec{J}$ is the total current density, free plus bound.

Problem domain and coordinate systems

We adopt the Cartesian coordinate system, the coordinates being $x$ , $y$ , and $z$ , as shown in the figure below. We denote the corresponding unitary basis vectors as $\hat{i}$ , $\hat{j}$ , and $\hat{k}$ . The problem domain consists of a ball centered at the origin and four spherical shells around it. The spheres $\Gamma_{I1}$ and $\Gamma_{I2}$ are interfaces between the magnetic material and free space. Strictly speaking, the magnetic field induced by the coil extends to infinity. Therefore, the problem we would like to solve is unbounded. As the mesh needs to have finite dimensions, we truncate the space with the sphere $\Gamma_{R1}$ . This sphere will represent infinity. We denote the spherical coordinates as $r$ , $\theta$ , and $\phi$ , with $\theta$ and $\phi$ being polar and azimuthal angles, respectively. The figure below illustrates the basis vectors of the spherical coordinate system, $\hat{r}$ , $\hat{\theta}$ , and $\hat{\phi}$ .

The domain

Exact solution

The problem of a permeable spherical shell exposed to a uniform magnetic field is ubiquitous in the literature on electromagnetics. The problem of the magnetic field generated by a spherical coil can be found in the literature as well. Remarkably, the magnetic field inside the coil is uniform. Due to this fact one can combine the solution to these two problems to get the exact closed-form analytical expressions for the magnetic field generated by the coil described above. The exact solution presented below is a combination of the solutions to the problems of permeable shell and spherical coil. The derivation of the equations presented below can be found in [uzunbajakau2024b].

We express the magnetic field induced by the coil as a sum of the magnetic field induced by the free-current density in free space, $\vec{B}_J$ , and the magnetic field induced by the magnetization of the magnetic core, $\vec{B}_{\mu}$ ,

$\‍[\vec{B} = \vec{B}_J +\vec{B}_{\mu}. \‍]$

It is convenient to express the term $\vec{B}_J$ in a spherical coordinate system,

$\‍[\vec{B}_J = \left\{ \begin{aligned} & \frac{1}{2} (b_2^2-a_2^2) \vec{F}_1(\theta) &&\text{if }&& r \le a_2\\ & \frac{1}{2} (b_2^2-r^2) \vec{F}_1(\theta) + \frac{1}{5} (r^5-a_2^5) \vec{F}_2(r, \theta) && \text{if } && a_2 \le r \le b_2\\ & \frac{1}{5} (b_2^5-a_2^5) \vec{F}_2(r, \theta) &&\text{if }&& r \ge b_2,\\ \end{aligned} \right. \‍]$

where

$\‍[\begin{aligned} &\vec{F}_1(\theta) = \dfrac{2}{3}\mu_0 K_0\big(\cos(\theta)\hat{r}-\sin(\theta)\hat{\theta}\big), \\ &\vec{F}_2(r, \theta) = \dfrac{2}{3}\mu_0 K_0 \frac{1}{r^3}\big(\cos(\theta)\hat{r}+\frac{1}{2}\sin(\theta)\hat{\theta}\big). \end{aligned} \‍]$

We, however, need to program in the Cartesian coordinate system. The expressions above can be converted into the Cartesian coordinate system by invoking the following identities:

$\‍[r = \sqrt{x^2+y^2+z^2}, \‍]$

$\‍[\cos(\theta) = \frac{z}{r}, \‍]$

$\‍[\sin(\theta) = \frac{\sqrt{x^2+y^2}}{r}, \‍]$

$\‍[\cos(\phi) = \frac{x}{\sqrt{x^2+y^2}}, \‍]$

$\‍[\sin(\phi) = \frac{y}{\sqrt{x^2+y^2}}, \‍]$

$\‍[\hat{r} = \frac{1}{r} (x\hat{i}+y\hat{j}+z\hat{k}), \‍]$

and

$\‍[\hat{\theta}= \cos(\theta)\cos{\phi}\hat{i} + \cos(\theta)\sin{\phi}\hat{j} - \sin(\theta) \hat{k}. \‍]$

The magnetic field induced by the magnetic core can be expressed in Cartesian coordinates as

$\‍[\vec{B}_{\mu} = - \mu\vec{\nabla} \Psi - \mu_0 H_0 \hat{k}, \‍]$

where

$\‍[H_0 = \frac{1}{3} K_0 (b_2^2-a_2^2), \‍]$

$\begin{equation}\mu = \left\{ \begin{aligned} & \mu_0 &\text{if } \text{ } & r < a_1\\ & \mu_1 &\text{if } \text{ } & a_1 \le r \le b_1\\ & \mu_0 &\text{if } \text{ } & r > b_1, \end{aligned} \right. \end{equation}$

and

$\begin{equation}\vec{\nabla} \Psi = \left\{ \begin{aligned} &\delta_1 \hat{k} &\text{if } \text{ } & r \le a_1\\ &-3\gamma_1\frac{xz}{r^5}\hat{i} -3\gamma_1\frac{yz}{r^5}\hat{j} + \bigg(\beta_1 + \gamma_1\frac{1}{r^3} - 3\gamma_1\frac{z^2}{r^5} \bigg)\hat{k} &\text{if } \text{ } & a_1 \le r \le b_1\\ &-3\alpha_1\frac{xz}{r^5}\hat{i} -3\alpha_1\frac{yz}{r^5}\hat{j} + \bigg(-H_0 + \alpha_1\frac{1}{r^3} - 3\alpha_1\frac{z^2}{r^5} \bigg)\hat{k} &\text{if } \text{ } & r \ge b_1. \end{aligned} \right. \end{equation}$

The following constants were used in the last equation:

$\begin{equation}\begin{aligned} &\Omega = \frac{(\mu_r - 1)}{(\mu_r + 2)} \frac{a_1^3}{b_1^3}, \\ &\gamma_1 = \frac{-3b_1^3 H_0 \Omega}{(2\mu_r+1)-2 (\mu_r-1) \Omega}, \\ &\beta_1 = \frac{(2\mu_r+1) \gamma_1}{(\mu_r-1)a_1^3}, \\ &\alpha_1 = \frac{-b_1^3 H_0+2\mu_r\gamma_1-\mu_r b_1^3\beta_1}{2}, \\ &\delta_1 = \frac{\mu_r a_1^3\beta_1-2\mu_r\gamma_1}{a_1^3}. \end{aligned} \end{equation}$

The relative permeability in the last five equations is computed as

$\‍[\mu_r = \frac{\mu_1}{\mu_0}. \‍]$

Boundary value problems

Magnetic vector potential - A

Ampere's law in static approximation can be written in terms of the auxiliary vector field $\vec{H}$ as

$\‍[\vec{\nabla}\times\vec{H} = \vec{J}_f, \‍]$

where $\vec{J}_f$ is the free current and $\vec{H}$ is given by

$\begin{equation}\begin{array}{ll} \vec{H} = \dfrac{1}{\mu}\vec{B} & \text{(constitutive relation)}. \end{array} \end{equation}$

We introduce the magnetic vector potential, $\vec{A}$ , as

$\begin{equation}\begin{array}{ll} \vec{B} = \vec{\nabla}\times\vec{A} & \text{(definition of magnetic vector potential)}. \end{array} \end{equation}$

By combining the last three equations we arrive at

$\‍[\vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A}\bigg) = \vec{J}_f, \‍]$

which is the curl-curl partial differential equation we would like to solve. There are, however, two subtleties we need to address before solving any curl-curl equations: the gauge and compatibility condition.

Gauge

The Helmholtz decomposition theorem suggests that any vector field studied in electromagnetics, say $\vec{A}$ , can be represented as a sum of a conservative and a solenoidal vector field,

$\‍[\vec{A} = \underbrace{- \vec{\nabla} V}_{\text{conservative}} + \underbrace{\vec{\nabla} \times \vec{W}}_{\text{solenoidal}} = \vec{C} + \vec{S}. \‍]$

The curl of a conservative vector field always equals zero (conservative fields are curl free),

$\‍[\vec{\nabla} \times \vec{C} = -\vec{\nabla} \times \big( \vec{\nabla} V \big) = 0. \‍]$

Consequently, any vector field derived by adding a conservative vector field to a solution to the curl-curl equation is also a solution,

$\‍[\vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\Big(\vec{A} - \vec{\nabla} V'\Big)\bigg) = \vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A}\bigg). \‍]$

We can restrict this freedom by specifying the divergence of the solution. Indeed, the divergence of a solenoidal vector field always equals zero (solenoidal fields are divergence free),

$\‍[\vec{\nabla} \cdot \vec{S} = \vec{\nabla} \cdot \big(\vec{\nabla} \times \vec{W}\big) = 0. \‍]$

Then the divergence of a solution can be expressed as

$\‍[\vec{\nabla} \cdot \vec{A} = \vec{\nabla} \cdot \big(\vec{C} + \vec{S}\big) = \vec{\nabla} \cdot \vec{C}. \‍]$

By specifying the divergence of the solution, we specify the divergence of the conservative part of the solution, and, thus, we specify the conservative part of the solution itself. In short: the magnetic vector potential consists of conservative and solenoidal components; the curl-curl equation fixes the solenoidal component, the gauge fixes the conservative component. For example, the Coulomb gauge,

$\‍[\vec{\nabla} \cdot \vec{A} = 0, \‍]$

specifies that there is no conservative component in $\vec{A}$ . That is, the Coulomb gauge selects a purely solenoidal solution. We are interested in the magnetic field, $\vec{B}$ . In the realm of classical electromagnetics, the magnetic vector potential, $\vec{A}$ , has an instrumental value - we just use it to compute the magnetic field, $\vec{B}$ . The magnetic vector potential is a useful mathematical tool but it has no physical meaning as it cannot be measured. On the contrary, the magnetic field, $\vec{B}$ , can be measured and, as a consequence, has a distinctive physical meaning. From a physical perspective, there are (infinitely) many magnetic vector potentials that lead to the same physically measurable magnetic field and for a physicist, it doesn't matter which of these potentials we choose. But from a mathematical perspective, we need to pick one to make the problem well-posed. A gauge is a way to pick one. Because the magnetic field is computed as

$\‍[\vec{B} = \vec{\nabla}\times\vec{A} = \vec{\nabla}\times\big(\vec{C} + \vec{S} \big) = \vec{\nabla}\times\vec{S}, \‍]$

we are interested in the solenoidal part of the solution. Which conservative component of the solution is brought by the gauge is not important.

A good question is: how to enforce the gauge? The literature suggests that the best approach is not to gauge the magnetic vector potential explicitly. Instead, one should invoke an advanced linear solver (tree-cotree splitting, domain decomposition, geometric multigrid, algebraic multigrid, etc.). For example, in [manges1995a] it is suggested that the tree-cotree splitting can be used us a gauge. The hypre AMS utilizes a more modern approach [hypre1998b]. It can solve positive semidefinite linear systems yielded by the ungauged curl-curl equation. In general, if a solver is able to choose one solution out of many ungauged solutions, we can claim that such solution is implicitly gauged by the linear solver. The linear solver will not communicate to us which gauge has been applied, i.e., the conservative part of the solution will be fixed, but we will not be able to tell what it is, exactly. We will call such gauge an implicit gauge. An implicit gauge is good enough for our purpose.

We will use the CG solver which is the next best choice to the advanced linear solvers mentioned above. The CG solver is designed to handle symmetric positive definite matrices. This means that all eigenvalues of the system matrix are expected to be positive. On the other hand, the curl-curl equation above being augmented with proper boundary conditions will always yield a symmetric positive semidefinite matrix. The last means that most of the eigenvalues of the system matrix will be positive but some of them will be zero. We elevate the zero eigenvalues above zero just a bit by adding a very small gauging term, $\eta^2\vec{A}$ , to the partial differential equation,

$\‍[\vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A}\bigg) + \eta^2\vec{A} = \vec{J}_f, \‍]$

so all eigenvalues become positive. Formally, the system matrix with elevated eigenvalues is symmetric positive definite. We can feed it to the CG solver. The partial differential equation is different now. We can expect a different solution. Moreover, $\eta^2 \vec{A}$ belongs to the $H(\text{curl})$ function space, while the rest of the terms of the partial differential equation belong to the $H(\text{div})$ function space, see the Bossavit's diagram below. They exhibit different behavior on interfaces between dissimilar materials. Strictly speaking, we destroy the compatibility between the two sides of the curl-curl equation by adding the gauging term. However, by choosing the parameter $\eta^2$ to be small enough we can make these horrible crimes insignificant. In practice, we need to set $\eta^2$ to zero. If the CG solver chokes, we will increase $\eta^2$ just a bit. Setting $\eta^2=\dfrac{10^{-6}}{\mu_0}$ can save the day. The $\eta^2$ -trick is an implicit gauge: it helps the CG solver to select one solution out of many, but it is impossible to predict what the conservative part of the solution will look like. The last is not important for the problem we would like to solve because, as mentioned above, the vector potential $\vec{A}$ has no physical reality and all choices of a gauge (including the implicit gauge used here) lead to the same physical fields.

Compatibility condition

To solve the curl-curl equation numerically, we need to feed to the solver the right-hand side of the equation, $\vec{J}_f$ . Whatever we do, we will end up feeding to the solver a discretized version of $\vec{J}_f$ which will be a combination of conservative and solenoidal vector fields even if the problem we would like to solve specifies a closed-form analytical expression of a purely solenoidal $\vec{J}_f$ . The conservative component will be added implicitly by discretization. On the other hand, the left-hand side of the curl-curl equation is purely solenoidal as the curl of any vector field is always solenoidal (assume $\eta^2 = 0$ for a moment). Therefore, the two sides of the equation will be somewhat incompatible: pure solenoidal vector field on the left-hand side and a combination of a conservative and solenoidal vector fields on the right-hand side. As a consequence, we cannot expect that we can get the norm of the residual to zero. In other words, the linear solver cannot converge to a solution that has zero residual. In practice, this means that CG will either not converge at all, or terminate only after an unreasonably large number of iterations. To make the two sides of the equation compatible, we derive the free-current density from a current vector potential, $\vec{T}$ , as

$\‍[\vec{J}_f = \vec{\nabla}\times\vec{T}. \‍]$

Consequently, the curl-curl equation becomes

$\‍[\vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A}\bigg) + \eta^2\vec{A} = \vec{\nabla}\times\vec{T}. \‍]$

This equation contains purely solenoidal vector fields on both sides if $\eta^2=0$ . Now, however, we need to convert $\vec{J}_f$ into $\vec{T}$ . There are a number of ways of doing so. One of them is solving another curl-curl equation as discussed below. For now we assume that $\vec{T}$ is known.

Boundary and interface conditions

Even if the magnetic vector potential is gauged and the compatibility issue is solved, one still needs to set up the boundary conditions to pick up one solution out of many gauged vector potentials that satisfy the curl-curl equation. Moreover, the behavior of the magnetic field on interfaces between dissimilar materials is defined by the Maxwell's equations. Consequently, the behavior of the magnetic vector potential on interfaces is defined by the Maxwell's equations as well. For this reason, we need to add the interface conditions next to the boundary conditions to ensure the correct behavior of the current vector potential on interfaces. The curl-curl partial differential equation augmented with the boundary and interface conditions constitutes the boundary value problem. We compose the following boundary value problem for the problem discussed above:

$\begin{equation}\begin{array}{lrcll} \text{ } & \vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A}\bigg) + \eta^2 \vec{A} = \vec{\nabla}\times\vec{T} & \text{in} & \Omega & \text{(i)},\\ \text{(n)}& \dfrac{1}{\mu}\hat{n}\times\bigg(\vec{\nabla}\times\vec{A}\bigg) + \gamma \hat{n}\times\bigg(\hat{n}\times\vec{A}\bigg) =0 &\text{on} & \Gamma_{R1} & \text{(ii)}, \\ \text{(e)}&\hat{n}\times\vec{A}_{+} = \hat{n}\times\vec{A}_{-}&\text{on}&\Gamma_{I1}\cup\Gamma_{I2}&\text{(iii)},\\ \text{(n)}&\dfrac{1}{\mu}_{+}\hat{n}\times \bigg( \vec{\nabla} \times \vec{A}_{+} \bigg) - \dfrac{1}{\mu}_{-}\hat{n}\times \bigg( \vec{\nabla} \times \vec{A}_{-} \bigg) = 0 & \text{on} & \Gamma_{I1}\cup\Gamma_{I2} & \text{(iv)}, \end{array} \end{equation}$

where

$\‍[\gamma = \frac{1}{\mu r}. \‍]$

In this boundary value problem the vector $\hat{n}$ is the unit vector normal to the surface on which the boundary or interface condition is defined. By convention, $\hat{n}$ always points outside a closed surface. The subscript "+" refers to the space immediately next to the interface in the direction of vector $\hat{n}$ . The subscript "-" refers to the space immediately next to the interface in the direction opposite to $\hat{n}$ .

Equation (i) in the boundary value problem has been discussed above and needs no further comments.

The surface $\Gamma_{R1}$ represents infinity. As most problems in magnetics, the problem we would like to solve is unbounded. That is to say, the magnetic field as well as the magnetic vector potential extend to infinity. A straightforward modeling of an unbounded domain by means of finite elements is impossible as such modeling will require an infinite amount of cells. To overcome this unfortunate predicament, we artificially truncate the unbounded problem domain with the surface $\Gamma_{R1}$ . This truncation will introduce a simulation error. We, however, will try to minimize this error by making the interior of the artificial surface $\Gamma_{R1}$ as spacious as possible. The magnetic field induced by the coil vanishes at infinity, so does the magnetic vector potential. For this reason, the homogeneous Dirichlet,

$\begin{equation}\begin{array}{ll} \hat{n} \times \vec{A} = 0 & \text{(Dirichlet)}, \end{array} \end{equation}$

and homogeneous Neumann,

$\begin{equation}\begin{array}{ll} \dfrac{1}{\mu}\hat{n}\times\bigg(\vec{\nabla}\times\vec{A}\bigg) = 0 & \text{(Neumann)}, \end{array} \end{equation}$

boundary conditions are used on $\Gamma_{R1}$ most often. For this approach to be efficient, one needs to choose the radius of the sphere $\Gamma_{R1}$ large enough so that the tangential component of the magnetic vector potential (Dirichlet) or the tangential component of the magnetic field (Neumann) is negligibly small on $\Gamma_{R1}$ . There is, however, a more efficient approach.

The first term of the multipole expansion of magnetic vector potential is the magnetic vector potential of a magnetic dipole moment unless the source of the magnetic field is explicitly configured as a higher-order multipole. For this reason, the magnetic vector potential induced by the coil looks more like that of a magnetic dipole at a sufficient distance. If a magnetic dipole is placed at the center of a sphere, the magnetic vector potential induced by the dipole on the sphere satisfies the following identity:

$\‍[\hat{n} \times \bigg( \vec{\nabla}\times \vec{A} \bigg) = - \frac{1}{r} \hat{n} \times \bigg( \hat{n} \times \vec{A} \bigg). \‍]$

Note that this is a purely geometric statement. It does not depend on how we choose our coordinate system. Consequently, the identity given by the last equation is true for a magnetic dipole of any orientation as long as it is situated at the center of the sphere. By multiplying the last equation by $\dfrac{1}{\mu}$ and rearranging the terms we arrive into

$\begin{equation}\begin{array}{ll} \dfrac{1}{\mu}\hat{n}\times\bigg(\vec{\nabla}\times\vec{A}\bigg) + \dfrac{1}{\mu r}\hat{n}\times\bigg(\hat{n}\times\vec{A}\bigg) =0 & \text{(Robin)}, \end{array} \end{equation}$

which is, essentially, the Robin boundary condition (ii) in the boundary value problem above. The Robin boundary condition in the current context is also called the first-order asymptotic boundary condition, first-order ABC for short. The Robin boundary condition is expected to be superior to both Dirichlet and Neumann boundary conditions as it allows to reach the same level of simulation error with a sphere $\Gamma_{R1}$ of a smaller radius.

All three boundary conditions, Dirichlet, Neumann, and Robin, guarantee the uniqueness of the curl of the solution, $\vec{\nabla}\times\vec{A}$ . The uniqueness of the solution itself, $\vec{A}$ , can only be guaranteed if one of the boundary conditions is applied in combination with a gauge, see the discussion above. The program uses the Robin boundary condition by default, but can easily be switched to Dirichlet or Robin boundary conditions.

The behavior of the magnetic field on the interfaces between dissimilar materials can be described by

$\begin{equation}\begin{array}{ll} \hat{n}\cdot\vec{B}_+ = \hat{n}\cdot\vec{B}_- ,\\ \hat{n}\times\vec{H}_+ - \hat{n}\times\vec{H}_- = \vec{K}_f, \end{array} \end{equation}$

where $\vec{K}_f$ is a surface free-current density that can be present on an interface. No free currents flow on the surface of the magnetic core of the coil described above. Therefore, we set $\vec{K}_f = 0$ . Then by substituting the constitutive relation for the magnetic field and the definition of the magnetic vector potential, see above, into the last two equations we get

$\begin{equation}\begin{array}{ll} \hat{n}\cdot\bigg(\vec{\nabla}\times\vec{A}_+\bigg) = \hat{n}\cdot\bigg(\vec{\nabla}\times\vec{A}_-\bigg) &\text{(condition A)},\\ \hat{n}\times\bigg(\dfrac{1}{\mu_+}\vec{\nabla}\times\vec{A}_+\bigg) - \hat{n}\times\bigg(\dfrac{1}{\mu_-}\vec{\nabla}\times\vec{A}_-\bigg) = 0. \end{array} \end{equation}$

Next, we simply replace the first condition with another condition and rearrange the terms of the second condition:

$\begin{equation}\begin{array}{ll} \hat{n}\times\vec{A}_+ = \hat{n}\times\vec{A}_- & \text{(condition B)}, \\ \dfrac{1}{\mu_+}\hat{n}\times\bigg(\vec{\nabla}\times\vec{A}_+\bigg) - \dfrac{1}{\mu_-}\hat{n}\times\bigg(\vec{\nabla}\times\vec{A}_-\bigg) = 0, \end{array} \end{equation}$

The last two equations are identical to equations (iii) and (iv) in the boundary value problem above. These are the interface conditions that must be observed on the inner and outer surfaces of the magnetic core, $\Gamma_{I1}\cup\Gamma_{I2}$ .

We have replaced condition A with condition B. This replacement deserves an explanation. The are two reasons for this replacement. The first reason is the following. A magnetic flux through an infinitesimally small closed loop must be infinitesimally small even if an interface between dissimilar magnetic materials passes through the loop. We can guarantee this only if we require that condition B holds. Condition A will not do. The second reason is the uniqueness of the curl of the solution. If we include condition A into the boundary value problem above, the uniqueness of the curl of the solution, $\vec{\nabla}\times\vec{A}$ , cannot be guaranteed. It is possible to guarantee it if condition B is included instead. A good question is: does this replacement disturb the behavior of the magnetic field on the interfaces? Not at all. Consider the following. The normal component of the curl, $\hat{n}\cdot\big(\vec{\nabla}\times\vec{A}\big)$ , is defined by the derivatives of the tangential components of the vector potential. If two vector potentials, $\vec{A}_-$ and $\vec{A}_+$ , have the same tangential components on the interface, i.e., condition B, the derivatives of these tangential components along the interface will be the same as well. Therefore, condition B implicitly implies condition A and, thus, the replacement is justified. We can conclude that condition B is more restrictive than condition A, so the behavior of the magnetic vector potential on the interfaces is observed. Note, that the interface conditions (iii) and (iv) in the boundary value problem above have absolutely no influence on the uniqueness of the solution regardless which boundary condition, Dirichlet, Neumann, or Robin, is used.

It remains to comment on the labels (e) and (n) in the boundary value problem above. They label (e)ssential and (n)atural boundary and interface conditions. How to sort conditions in these two categories and why it is important to do so is discussed below in the section on variational formulations.

Current vector potential - T

To be able to solve the boundary value problem above, we need to convert the free current density, $\vec{J}_f$ , into current vector potential, $\vec{T}$ . We define the current vector potential as

$\‍[\vec{J}_f = \vec{\nabla}\times\vec{T}. \‍]$

By taking the curl of the last equation and rearranging the terms we arrive into

$\‍[\vec{\nabla}\times\bigg(\vec{\nabla}\times\vec{T}\bigg) = \vec{\nabla}\times\vec{J}_f. \‍]$

We will compute the current vector potential by solving this equation. As soon as there are curls of vector fields on both sides of the equation, the compatibility condition is observed. We add to this equation a gauging term,

$\‍[\vec{\nabla}\times\bigg(\vec{\nabla}\times\vec{T}\bigg) + \eta^2 \vec{T} = \vec{\nabla}\times\vec{J}_f. \‍]$

The discussion on the gauging term in the preceding section is valid in the current context as well.

The current vector potential is computed in free space, i.e., $\mu = \mu_0$ everywhere. For this reason, there is no need for interface conditions as there are no interfaces in free space.

We still need to apply a boundary condition to the surface $\Gamma_{R1}$ . It is beneficial to apply the homogeneous Dirichlet boundary condition as it allows to nullify the integral $I_{b3-2}$ in the recipe for calculating the magnetic vector potential, see below. Consequently, we can neglect this integral when programming and this makes the code more efficient.

We compose a boundary value problem by combining the curl-curl equation for the current vector potential with the Dirichlet boundary condition:

$\begin{equation}\begin{array}{lrcll} \text{ } & \vec{\nabla}\times\bigg(\vec{\nabla}\times\vec{T}\bigg) + \eta^2 \vec{T} = \vec{\nabla}\times\vec{J}_f & \text{in} & \Omega & \text{(i)},\\ \text{(e)}& \hat{n}\times\vec{T} = 0 &\text{on} & \Gamma_{R1} & \text{(ii)}. \end{array} \end{equation}$

Variational formulations

Magnetic vector potential

The problem of solving the boundary value problem for the magnetic vector potential can be replaced by the problem of minimizing the following functional:

$\begin{equation}\begin{aligned} &F(\vec{A}) = \iiint_{\Omega}\frac{1}{\mu}\bigg|\vec{\nabla}\times\vec{A}\bigg|^2 dV + \iint_{\Gamma_{R1}} \gamma \bigg|\hat{n}\times\vec{A}\bigg|^2 dS + \eta^2 \iiint_{\Omega}\mid\vec{A}\mid^2 dV -2\iiint_{\Omega} \vec{T} \cdot \bigg( \vec{\nabla} \times\vec{A} \bigg) dV +2\iint_{\Gamma_{R1}} \vec{T} \cdot \bigg(\hat{n}\times\vec{A}\bigg) dS. \end{aligned} \end{equation}$

Any magnetic vector potential that minimizes this functional will satisfy the curl-curl equation (i), the boundary condition (ii) and the interface condition (iv). On the other hand, this functional is invariant to the interface condition (iii). This can be verified as follows: The first variation of the functional at the minimum vanishes,

$\‍[\delta F(\vec{A}) = 0. \‍]$

We can express the first variation and deduce that the identities given by equations (i), (ii), and (iv) in the boundary value problem must be satisfied for the last equation to be true. The interface condition (iii), however, is not needed for the last equation to be satisfied. The boundary and interface conditions that are needed to drive the first variation of the functional to zero are called natural and are labeled by (n) in the boundary value problem. If a condition is not natural, it is essential. We label essential conditions by (e).

The natural conditions are enforced by minimization of the functional. The essential boundary and interface conditions must be enforced by other means. The interface condition (iii) in the boundary value problem above is enforced by choosing the finite elements wisely. The FE_Nedelec finite elements ensure continuity of the tangential component of the vector field they model on the faces of the cells. Therefore, if we construct the mesh such that all interfaces are made of cell faces, i.e., no interface runs through a cell, and choose the FE_Nedelec finite elements to model the magnetic vector potential, the interface condition (iii) will be enforced automatically. The Dirichlet boundary condition is essential. It is enforced by constraining the degrees of freedom in the system of linear equations. This can be done by invoking the deal.II function VectorTools::project_boundary_values_curl_conforming_l2(). The other two boundary conditions, Neumann and Robin, are natural.

Note that the homogeneous Neumann boundary condition is embedded into the first integral of the functional. This integral is present in all functionals related to the curl-curl equation, even in the most minimalistic ones. Therefore, application of no boundary condition is quite impossible: the homogeneous Neumann boundary condition,

$\‍[\frac{1}{\mu}\hat{n}\times\bigg(\vec{\nabla}\times\vec{A}\bigg) = 0, \‍]$

will be applied implicitly by default as a result of the functional minimization.

We can switch between the three boundary conditions as the following.

The functional above as it is corresponds to a boundary value problem with Robin boundary condition, i.e., the boundary value problem given above.
To switch to Neumann boundary condition, we need to discard the second integral from the functional. Such functional will correspond to the boundary value problem above with equation (ii) replaced by the Neumann boundary condition.
To switch to Dirichlet boundary condition, we use the functional with discarded second integral and use VectorTools::project_boundary_values_curl_conforming_l2() to constrain the system of linear equations. Such configuration will correspond to the boundary value problem above with equation (ii) replaced by the Dirichlet boundary condition.

In general, we need to sort all boundary and interface conditions into two categories, natural and essential. The reason for doing so is simple: we need to know which conditions are not taken care of by the functional (the essential conditions), so we can take care of them by other means (choice of finite elements, restricting dofs). Interrogating the first variation of the functional as discussed above is a common method of sorting the conditions. A more detailed discussion on boundary conditions can be found in step-22. (See also video lecture 21.5, video lecture 21.55, video lecture 21.6, video lecture 21.65.)

Next, we need to convert the functional above into a numerical recipe that can be programmed into a computer. To do so, we note that the final result, i.e., the magnetic vector potential computed by the finite element method, is represented as a sum,

$\‍[\vec{A}(\vec{r}) = \sum_{j=0}^{m-1} c_j \vec{N}_j(\vec{r}), \‍]$

where $c_j$ are the degrees of freedom and $\vec{N}_j$ are the vector-valued shape functions of the FE_Nedelec finite elements. We convert the functional into a multivariate function by substituting the last equation into the functional,

$\‍[f(c_0, c_1, ..., c_{m-1}) = F\bigg(\sum_{j=0}^{m-1} c_j \vec{N}_j\bigg). \‍]$

Then we compose the system of linear equations by computing the partial derivatives of the multivariate function:

$\begin{equation}\begin{array}{lr} \dfrac{\partial}{\partial c_0} f(c_0, c_1, ..., c_{m-1})& = 0, \\ \dfrac{\partial}{\partial c_1} f(c_0, c_1, ..., c_{m-1})& = 0, \\ ... & \\ \dfrac{\partial}{\partial c_{m-1}} f(c_0, c_1, ..., c_{m-1})& = 0. \end{array} \end{equation}$

This system of linear equations can be written in matrix form as

$\‍[\boldsymbol{A} \boldsymbol{c} = \boldsymbol{b}, \‍]$

where $\boldsymbol{c}$ is a column vector filled with the degrees of freedom,

$\‍[\boldsymbol{c} = [c_0, c_1, ... c_{m-1}]^T. \‍]$

The system matrix and the right-hand side column vector are computed as

$\begin{equation}\begin{aligned} & A_{ij} = \underbrace{\iiint_{\Omega}\frac{1}{\mu} \bigg(\vec{\nabla}\times\vec{N}_i\bigg) \cdot \bigg(\vec{\nabla}\times\vec{N}_j\bigg) dV}_{I_{a1}} + \underbrace{\iint_{\Gamma_{R1}} \gamma \bigg(\hat{n}\times \vec{N}_i \bigg) \cdot \bigg(\hat{n}\times \vec{N}_j \bigg) dS}_{I_{a2}} + \underbrace{\eta^2\iiint_{\Omega} \vec{N}_i \cdot \vec{N}_j dV}_{I_{a3}} \end{aligned} \end{equation}$

and

$\begin{equation}b_i = \underbrace{\iiint_{\Omega}\vec{T}\cdot\bigg(\vec{\nabla}\times\vec{N}_i\bigg) dV }_{I_{b3-1}} - \underbrace{\iint_{\Gamma_{R1}}\vec{T}\cdot\bigg(\hat{n}\times\vec{N}_i\bigg)dS }_{I_{b3-2} = 0}. \end{equation}$

The last two equations are implemented by the computer code. The integral $I_{b3-2}$ equals zero as the tangential component of the current vector potential, $\hat{n} \times \vec{T}$ , is forced to zero on the boundary $\Gamma_{R1}$ by the Dirichlet boundary condition. We will neglect this integral in the computer code.

Current vector potential

The problem of solving the boundary value problem for the current vector potential can be replaced by the problem of minimizing the following functional:

$\begin{equation}\begin{aligned} &F(\vec{T}) = \iiint_{\Omega}\bigg|\vec{\nabla}\times\vec{T}\bigg|^2 dV + \eta^2 \iiint_{\Omega}\mid\vec{T}\mid^2 dV -2\iiint_{\Omega} \vec{J}_f \cdot \bigg( \vec{\nabla} \times\vec{T} \bigg) dV +2\iint_{\Gamma_{R1}} \vec{J}_f \cdot \bigg(\hat{n}\times\vec{T}\bigg) dS, \end{aligned} \end{equation}$

and constraining the system of linear equations with a help of VectorTools::project_boundary_values_curl_conforming_l2().

This functional is minimized by solving a system of linear equations with the following system matrix and right-hand side column vector:

$\begin{equation}\begin{aligned} & A_{ij} = \underbrace{\iiint_{\Omega} \bigg(\vec{\nabla}\times\vec{N}_i\bigg) \cdot \bigg(\vec{\nabla}\times\vec{N}_j\bigg) dV}_{I_{a1}} + \underbrace{\eta^2\iiint_{\Omega} \vec{N}_i \cdot \vec{N}_j dV}_{I_{a3}}, \end{aligned} \end{equation}$

$\begin{equation}b_i = \underbrace{\iiint_{\Omega}\vec{J}_f\cdot\bigg(\vec{\nabla}\times\vec{N}_i\bigg) dV }_{I_{b3-1}} - \underbrace{\iint_{\Gamma_{R1}}\vec{J}_f\cdot\bigg(\hat{n}\times\vec{N}_i\bigg)dS }_{I_{b3-2}=0}. \end{equation}$

Here again, the integral $I_{b3-2}$ equals zero as the free-current density, $\vec{J}_f$ , equals zero at the boundary $\Gamma_{R1}$ by definition of the problem. We will neglect this integral in the computer code.

Converting potentials into fields

As discussed above, we are interested in the magnetic field induced by the coil. For this reason, we convert the numerically computed vector potential into magnetic field as

$\‍[\vec{B} = \vec{\nabla} \times \vec{A}. \‍]$

The problem of computing this equation can be replaced by the problem of minimizing the following functional:

$\‍[F(\vec{B}) = \iiint_{\Omega} \big| \vec{B} \big|^2 dV - 2 \iiint_{\Omega} \bigg( \vec{\nabla} \times \vec{A} \bigg) \cdot \vec{B} dV. \‍]$

This functional is minimized by solving a system of linear equations with the system matrix

$\‍[A_{ij} = \underbrace{\iiint_{\Omega} \vec{N}_i \cdot \vec{N}_j dV}_{I_a} \‍]$

and right-hand side column vector

$\‍[b_i = \underbrace{\iiint_{\Omega} \bigg( \vec{\nabla} \times \vec{A} \bigg) \cdot \vec{N}_i dS}_{I_b}. \‍]$

This time, however, $\vec{N}_i$ are the shape functions of the FE_RaviartThomas finite elements. The vector field $\vec{A}$ in the last equation is the numerically computed magnetic vector potential, i.e., a linear combination of the shape functions of the FE_Nedelec finite elements.

It is informative to verify the quality of the computed current vector potential, $\vec{T}$ . We can do this by converting it back into the free-current density, $\vec{J}_f$ , and comparing the result to the closed-form analytical expression given by the definition of the problem, see above. The current vector potential relates to the free-current density as

$\‍[\vec{J}_f = \vec{\nabla} \times \vec{T}. \‍]$

This equation is, essentially, the same as the equation for converting $\vec{A}$ into $\vec{B}$ , see the first equation of this section. Consequently, we can simply adapt the equations for the functional, system matrix, and the right-hand side given above by making the following substitutions:

$\begin{equation}\begin{array}{lcr} \vec{A} \rightarrow \vec{T} &\text{ and }& \vec{B} \rightarrow \vec{J}_f. \end{array} \end{equation}$

The two conversions, $\vec{A}$ into $\vec{B}$ and $\vec{T}$ into $\vec{J}_f$ , can be done by the same piece of code. We will call it a projector from $H(\text{curl})$ to $H(\text{div})$ .

Selecting finite elements

The preceding sections leave an impression that we are going to model $\vec{A}$ and $\vec{T}$ by the FE_Nedelec finite elements, while $\vec{B}$ and $\vec{J}_f$ by the FE_RaviartThomas finite elements. It makes sense to discuss how this kind of choice can be made.

The simplest method of assigning a particular type of finite elements to a physical quantity studied in electromagnetics is a contemplation of the Bossavit's diagram. It is presented below. The Bossavit's diagram describes the relations between physical quantities as they are given by Maxwell's equations and constitutive relations. Most equations derived from Maxwell's equations can be captured by this diagram.

The Bossavit's diagram

Let us consider assigning the correct type of finite elements to the magnetic vector potential, $\vec{A}$ . First, we contemplate the Bossavit's diagram and observe that $\vec{A}$ belongs to the $H(\text{curl})$ function space. Second, we look at the table below and conclude that the physical quantities that belong to the $H(\text{curl})$ function space are modeled by the FE_Nedelec finite elements. Therefore, we need to model the magnetic vector potential, $\vec{A}$ , by the FE_Nedelec finite elements. The same assigning procedure can be applied to the rest of the physical quantities, $\vec{B}$ , $\vec{T}$ , and $\vec{J}_f$ .

Funct. space	Finite elements
$H(\text{grad})$	FE_Q
$H(\text{curl})$	FE_Nedelec
$H(\text{div})$	FE_RaviartThomas
	FE_DGQ

Alternatively, one of the types of the finite elements listed in the table can be assigned to a physical quantity by considering the behavior of the physical quantity on the interfaces between dissimilar materials. Let us consider the magnetic field $\vec{B}$ as an example. The normal component of the magnetic field is continuous on interfaces. On the contrary, the tangential component is discontinuous, see the first figure on this page. This type of behavior can be modeled by the FE_RaviartThomas finite elements as they guarantee the continuity of the normal component of the vector field on the faces of mesh cells.

Mesh

The mesh is, essentially, a discretization of the problem domain shown above. The mesh is constructed such that all spherical surfaces of the problem domain are delineated by cell faces. That is, no spherical surface runs through a mesh cell. The figure below illustrates the mesh. The main elements of the mesh are listed below.

The cube with a side of in the middle of the mesh. The purpose of this cube is instrumental. This is the only way to construct a mesh that has spherical interfaces and contains no tetrahedral cells as the cells around the origin must be hexahedral.
The spherical shell nr.1. This shell represents the magnetic core. The permeability of this region is $\mu_1$ . The permeability outside this region is $\mu_0$ .
The spherical shell nr.2. This region contains the free current. There is no free current outside this region.
Sphere nr.1. The error norms are computed inside the region delineated by this surface.
Sphere nr.2. The boundary of the problem domain. Represents infinity. The boundary condition is applied to this surface.

The program uses four meshes created with a help of gmsh. Each mesh has a different degree of refinement. The degree of refinement in each mesh is defined by a number of mesh nodes on the transfinite lines, $r$ . The figure below illustrates a mesh with $r=5$ . The program uses four meshes with the amount of nodes on the transfinite lines in the interval $r = [6, 9]$ .

Note that this tutorial uses only globally refined meshes. There are no non-conforming cells, hanging nodes, and hanging node constrains in this tutorial. We will use constrains only for the purpose of enforcing the Dirichlet boundary condition.

The mesh

Overview of the program

The program runs in a loop. In each iteration of the loop a different parameter $r$ (the mesh refinement parameter, i.e., the amount of nodes on transfinite lines) is assumed. Each iteration consists of four stages:

The current vector potential, $\vec{T}$ , is computed given the closed-form analytical expression for the free-current density, $\vec{J}_f$ , see above. No error norms are computed. The mesh is loaded at this stage. This mesh is simply reused at all other stages.
The current vector potential computed at the preceding stage, $\vec{T}$ , is converted numerically back into free-current density, $\vec{J}_f$ , for verification. The error norm is computed by comparing the numerical result with the closed-form analytical expression. The error norm is saved in a convergence table table_Jf.
The magnetic vector potential, $\vec{A}$ , is computed given the current vector potential computed at the first stage, $\vec{T}$ . No error norms are computed at this stage.
The magnetic vector potential computed at the preceding stage, $\vec{A}$ , is converted into magnetic field, $\vec{B}$ . The norm is computed by comparing the numerical result with the closed-form analytical solution. The error norm is saved in a convergence table table_B.

The first stage is implemented by the code contained in the name space SolverT. The third stage is implemented by the code contained in the name space SolverA. The second and the fourth stages are implemented by the code contained in the name space ProjectorHcurlToHdiv. The name space ExactSolutions contains exact closed-form analytical expressions for $\vec{B}$ and $\vec{J}_f$ .

No error norms are computed at the first and the third stages. Computing error norms at these stages makes no sense as we apply implicit gauges and the conservative portion of the solution is unknown. We can, however, judge the quality of the simulations at these two stages indirectly by evaluating the quality of the simulations at the second and the fourth stages.

The commented program

  #include <deal.II/base/tensor.h>
  #include <deal.II/base/function.h>
  #include <deal.II/base/convergence_table.h>
  #include <deal.II/base/timer.h>
  #include <deal.II/base/types.h>
  #include <deal.II/base/work_stream.h>
  #include <deal.II/base/multithread_info.h>
  
  #include <deal.II/grid/tria.h>
  #include <deal.II/grid/grid_in.h>
  #include <deal.II/grid/grid_out.h>
  #include <deal.II/grid/manifold_lib.h>
  
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/dofs/dof_tools.h>
  
  #include <deal.II/fe/fe_nedelec.h>
  #include <deal.II/fe/fe_raviart_thomas.h>
  #include <deal.II/fe/mapping_q.h>
  #include <deal.II/fe/fe_values.h>
  
  #include <deal.II/lac/vector.h>
  #include <deal.II/lac/sparse_matrix.h>
  #include <deal.II/lac/affine_constraints.h>
  #include <deal.II/lac/sparsity_pattern.h>
  #include <deal.II/lac/dynamic_sparsity_pattern.h>
  #include <deal.II/lac/solver_control.h>
  #include <deal.II/lac/vector_memory.h>
  #include <deal.II/lac/solver_cg.h>
  #include <deal.II/lac/precondition.h>
  
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/numerics/data_component_interpretation.h>
  #include <deal.II/numerics/data_out.h>
  
  #include <string>
  #include <ostream>
  
  using namespace dealii;
  

This enumeration is used to switch between boundary conditions when solving for the magnetic vector potential, $\vec{A}$ . The boundary condition is fixed (Dirichlet) when solving for the current vector potential, $\vec{T}$ .

  enum BoundaryConditionType
  {
    Dirichlet,
    Neumann,
    Robin
  };
  

Settings

The following is the control panel of the program. The scaling of the program can be changed by setting mu_0 = 1.0. Then the computed magnetic field, $\vec{B}$ , will have to be multiplied by a factor of $\mu_0 = 1.25664 \cdot 10^{-6}$ . The free-current density, $\vec{J}_f$ , does not depend on scaling. The setting boundary_condition_type_A can be used to switch between Dirichlet, Neumann, and Robin (ABC) boundary conditions. This has an effect only on the solver that solves for the magnetic vector potential, $\vec{A}$ . The solver that solves for the current vector potential, $\vec{T}$ , applies the Dirichlet boundary condition. This cannot be changed. Recall, that forcing to zero the tangential component of $\vec{T}$ on the boundary allows us to skip the boundary integral $I_{b3-2}$ . The boundary ID is set in the geo files that describe the mesh geometry. This is done by specifying the physical surface, i.e.,

Physical Surface(2) = {115, 479, 665, 297, 847, 1029};

The boundary ID in the geo file must match the boundary ID setting below. If project_exact_solution = true, the program projects the exact solutions for $\vec{J}_f$ and $\vec{B}$ onto $H(\text{div})$ function space and saves the results into corresponding vtu files next to the numerical solutions. If the projected exact solution and the numerical solution look alike, it is a good sign. In case of problems with the convergence of the CG solver the parameter $\eta^2$ must be increased just a bit.

  namespace Settings
  {
    const double permeability_fs = 1.2566370614359172954e-6;
    const double mu_0            = permeability_fs; // Permeability of free space.
  
    const double mu_r = 4; // Relative permeability of the magnetic material.
    const double mu_1 = mu_0 * mu_r; // Permeability of the magnetic material.
    const double d1   = 0.1; // Half-side of the cube in the center of the mesh.
    const double a1   = 0.3; // Inner radius of the magnetic core.
    const double b1   = 0.6; // Outer radius of the magnetic core.
    const double a2   = 0.9; // Inner radius of the free-current region.
    const double b2   = 1.2; // Outer radius of the free-current region.
    const double d2   = 2.0; // Error norms are computed if r < d2.
    const double K0   = 1.0; // Magnitude of the free-current density.
    const double H0   = // Magnitude of H-field inside the coil, core removed.
      (1.0 / 3.0) * K0 * (std::pow(b2, 2) - std::pow(a2, 2));
    const types::material_id material_id_free_space =
      1; // Material ID of the free space.
    const types::material_id material_id_core =
      2; // Material ID of the magnetic core.
    const types::material_id material_id_free_current =
      3; // Material ID of the free-current region.
    const types::boundary_id boundary_id_infinity =
      2; // Boundary ID. Set in the geo file.
  
    const unsigned int mapping_degree = 2; // Mapping degree used in all solvers.
    const unsigned int fe_degree      = 0; // Degree of the finite elements.
    const BoundaryConditionType boundary_condition_type_A = Robin;
  
    const double eta_squared_T = 0.0; // eta^2 when solving for T.
    const double eta_squared_A = 0.0; // eta^2 when solving for A.
  
    const unsigned int n_threads_max = 0; // If >0 limits the number of threads.
    const double eps = 1e-12; // Two doubles are equal if their difference < eps.
    const bool   log_cg_convergence = false; // Save the CG convergence data.
    const bool print_time_tables = false; // Print the time tables on the screen.
    const bool project_exact_solution = false; // Save the exact solution.
  } // namespace Settings
  

Next comes the weight function used to limit the region in which the $L^2$ error norms are computed. The error norms are computed if $r < d2$ .

  class Weight : public Function<3>
  {
  public:
    virtual double value(const Point<3> &p,
                         const unsigned int) const override final
    {
      if (p.norm() > Settings::d2)
        return 0.0;
  
      return 1.0;
    }
  };
  

This class describes a convergence table. The convergence tables are saved on disk in TeX format.

  class MainOutputTable : public ConvergenceTable
  {
  public:
    MainOutputTable() = delete;
  
    MainOutputTable(const unsigned int dimensions)
      : ConvergenceTable()
      , dimensions(dimensions)
    {}
  
    void set_new_order(const std::vector<std::string> &new_order_in)
    {
      new_order = new_order_in;
    }
  
    void append_new_order(const std::string &new_column)
    {
      new_order.push_back(new_column);
    }
  
    void format()
    {
      set_precision("L2", 2);
  
      set_scientific("L2", true);
  
      evaluate_convergence_rates("L2",
                                 "ncells",
                                 ConvergenceTable::reduction_rate_log2,
                                 dimensions);
  
      set_tex_caption("p", "p");
      set_tex_caption("r", "r");
      set_tex_caption("ncells", "nr. cells");
      set_tex_caption("ndofs", "nr. dofs");
      set_tex_caption("L2", "L2 norm");
  
      set_column_order(new_order);
    }
  
    void save(const std::string &fname)
    {
      format();
  
      std::ofstream ofs(fname + ".tex");
      write_tex(ofs);
    }
  
  private:
    const unsigned int       dimensions;
    std::vector<std::string> new_order = {"p", "r", "ncells", "ndofs", "L2"};
  };
  

Equations

This name space contains all closed-form analytical expressions mentioned in the introduction to this tutorial.

  namespace ExactSolutions
  {
  

This function describes the free-current density, $\vec{J}_f$ , inside the current region.

    inline Tensor<1, 3> volume_free_current_density(const Point<3> &p,
                                                    const double    K0)
    {
      return Tensor<1, 3>({-K0 * p[1], K0 * p[0], 0.0});
    }
  

This function computes the magnetic field induced by the free current in absence of the magnetic core, see the expression for $\vec{B}_J$ in the introduction.

    inline Tensor<1, 3> magnetic_field_coil(const Point<3> &p,
                                            const double    K0,
                                            const double    mu_0,
                                            const double    a2,
                                            const double    b2)
    {
      const double r = p.norm();
  
      const double cos_theta = p[2] / r;
      const double sin_theta =
        std::sqrt(std::pow(p[0], 2) + std::pow(p[1], 2)) / r;
  
      const double cos_phi =
        p[0] / std::sqrt(std::pow(p[0], 2) + std::pow(p[1], 2));
      const double sin_phi =
        p[1] / std::sqrt(std::pow(p[0], 2) + std::pow(p[1], 2));
  
      const Tensor<1, 3> r_hat({p[0] / r, p[1] / r, p[2] / r});
      const Tensor<1, 3> theta_hat(
        {cos_theta * cos_phi, cos_theta * sin_phi, -sin_theta});
  
      const Tensor<1, 3> F1 =
        (2.0 / 3.0) * mu_0 * K0 * (cos_theta * r_hat - sin_theta * theta_hat);
      const Tensor<1, 3> F2 = (2.0 / 3.0) * mu_0 * K0 *
                              (cos_theta * r_hat + 0.5 * sin_theta * theta_hat) /
                              std::pow(r, 3);
  
      if (r <= a2)
        {
          return 0.5 * (std::pow(b2, 2) - std::pow(a2, 2)) * F1;
        }
      else if (r >= b2)
        {
          return 0.2 * (std::pow(b2, 5) - std::pow(a2, 5)) * F2;
        }
      else
        {
          return 0.5 * (std::pow(b2, 2) - std::pow(r, 2)) * F1 +
                 0.2 * (std::pow(r, 5) - std::pow(a2, 5)) * F2;
        }
  
      return Tensor<1, 3>({0.404, 0.404, 0.404});
    }
  

This function computes the magnetic field induced by the magnetic core, see the expression for $\vec{B}_{\mu}$ in the introduction.

    inline Tensor<1, 3> magnetic_field_core(const Point<3> &p,
                                            const double    H0,
                                            const double    mur,
                                            const double    mu0,
                                            const double    a1,
                                            const double    b1)
    {
      const double a3 = std::pow(a1, 3);
      const double b3 = std::pow(b1, 3);
  
      const double OMEGA   = ((mur - 1.0) / (mur + 2.0)) * (a3 / b3);
      const double gamma_1 = (-3.0 * b3 * H0 * OMEGA) /
                             ((2.0 * mur + 1.0) - 2.0 * (mur - 1.0) * OMEGA);
      const double beta_1 = ((2.0 * mur + 1.0) * gamma_1) / ((mur - 1.0) * a3);
      const double alpha_1 =
        (-b3 * H0 + 2.0 * mur * gamma_1 - mur * b3 * beta_1) / 2.0;
      const double delta_1 = (mur * a3 * beta_1 - 2.0 * mur * gamma_1) / a3;
  
      const double r  = p.norm();
      const double r3 = std::pow(r, 3);
      const double r5 = std::pow(r, 5);
  
      const double zz = -3.0 * p[2] * p[2] / r5;
      const double xz = -3.0 * p[0] * p[2] / r5;
      const double yz = -3.0 * p[1] * p[2] / r5;
  
      if (r <= a1)
        {
          return -mu0 * (Tensor<1, 3>({0.0, 0.0, delta_1}) +
                         Tensor<1, 3>({0.0, 0.0, H0}));
        }
      else if (r >= b1)
        {
          return -mu0 * (Tensor<1, 3>({alpha_1 * xz,
                                       alpha_1 * yz,
                                       -H0 + alpha_1 / r3 + alpha_1 * zz}) +
                         Tensor<1, 3>({0.0, 0.0, H0}));
        }
      else
        {
          return -mu0 *
                 (mur * Tensor<1, 3>({gamma_1 * xz,
                                      gamma_1 * yz,
                                      beta_1 + gamma_1 / r3 + gamma_1 * zz}) +
                  Tensor<1, 3>({0.0, 0.0, H0}));
        }
  
      return Tensor<1, 3>({0.404, 0.404, 0.404});
    }
  

This class implements the closed-form analytical expression for the free-current density, $\vec{J}_f$ , in the entire domain.

    class FreeCurrentDensity : public Function<3>
    {
    public:
      FreeCurrentDensity()
        : Function<3>(3)
      {}
  
      virtual void
      vector_value_list(const std::vector<Point<3>> &p,
                        std::vector<Vector<double>> &values) const override final
      {
        Assert(values.size() == p.size(),
               ExcDimensionMismatch(values.size(), p.size()));
  
        double r;
  
        for (unsigned int i = 0; i < values.size(); i++)
          {
            r = p[i].norm();
  
            if ((r >= Settings::a2) && (r <= Settings::b2))
              {
                const Tensor<1, 3> Jf =
                  volume_free_current_density(p[i], Settings::K0);
  
                values[i][0] = Jf[0];
                values[i][1] = Jf[1];
                values[i][2] = Jf[2];
              }
            else
              {
                values[i][0] = 0.0;
                values[i][1] = 0.0;
                values[i][2] = 0.0;
              }
          }
      }
    };
  

This class implements the closed-form analytical expression for the magnetic field induced by the coil, $\vec{B} = \vec{B}_J + \vec{B}_{\mu}$ , see the introduction.

    class MagneticField : public Function<3>
    {
    public:
      MagneticField()
        : Function<3>(3)
      {}
  
      virtual void
      vector_value_list(const std::vector<Point<3>> &p,
                        std::vector<Vector<double>> &values) const override final
      {
        Assert(values.size() == p.size(),
               ExcDimensionMismatch(values.size(), p.size()));
  
        Tensor<1, 3> B;
  
        for (unsigned int i = 0; i < values.size(); i++)
          {
            using namespace Settings;
            B = magnetic_field_coil(p[i], K0, mu_0, a2, b2) +
                magnetic_field_core(p[i], H0, mu_r, mu_0, a1, b1);
  
            values[i][0] = B[0];
            values[i][1] = B[1];
            values[i][2] = B[2];
          }
      }
    };
  
  } // namespace ExactSolutions
  

Solver - T

This name space contains all the code related to the computation of the current vector potential, $\vec{T}$ .

namespace SolverT

{

This class describes the free-current density, $\vec{J}_f$ , on the right-hand side of the curl-curl equation, see equation (i) in the boundary value problem for $\vec{T}$ . The free-current density is given as a closed-form analytical expression by the definition of the problem.

    class FreeCurrentDensity
    {
    public:
      void value_list(const std::vector<Point<3>> &p, // Quadrature points.
                      const types::material_id   mid, // Material ID of the cell.
                      std::vector<Tensor<1, 3>> &values) const
      {
        Assert(p.size() == values.size(),
               ExcDimensionMismatch(p.size(), values.size()));
  
        if ((mid == Settings::material_id_free_space) ||
            (mid == Settings::material_id_core))
          for (unsigned int i = 0; i < values.size(); i++)
            values[i] = Tensor<1, 3>({0.0, 0.0, 0.0});
  
        if (mid == Settings::material_id_free_current)
          for (unsigned int i = 0; i < values.size(); i++)
            values[i] =
              ExactSolutions::volume_free_current_density(p[i], Settings::K0);
      }
    };
  

This class implements the solver that minimizes the functional $F(\vec{T})$ , see the introduction. The mesh is loaded in this class. All other solvers use a reference to this mesh.

      class Solver
      {
      public:
        Solver() = delete;
        Solver(const unsigned int p, // Degree of the Nedelec finite elements.
               const unsigned int r, // The mesh refinement parameter.
               const unsigned int mapping_degree,
               const double       eta_squared = 0.0,
               const std::string &fname       = "data");
    
        void make_mesh();  // Loads the mesh. Assigns mat. IDs. Attaches manifold.
        void setup();      // Initializes dofs, vectors, matrices.
        void assemble();   // Assembles the system of linear equations.
        void solve();      // Solves the system of linear equations.
        void save() const; // Saves computed T 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 six functions in the proper order
*                      and measures the execution time for each function. */
    

These three get-functions are used to channel the mesh and the solution to the next solver.

      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:

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.

      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;
      TimerOutput        timer;
  

The program utilizes the WorkStream technology. The step-9 tutorial does a much better job of explaining the workings of WorkStream. Reading the "WorkStream paper", see the glossary, is recommended. The following structures and functions are related to WorkStream.

      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;
        const FEValuesExtractors::Vector ve;
      };
  
      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)
      , timer(std::cout,
              (Settings::print_time_tables) ? TimerOutput::summary :
                                              TimerOutput::never,
              TimerOutput::cpu_and_wall_times_grouped)
    {}
  

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.

      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);
      }
    

This function initializes the dofs, applies the Dirichlet boundary condition, and initializes the vectors and matrices.

    void Solver::setup()
    {
      dof_handler.reinit(triangulation);
      dof_handler.distribute_dofs(fe);
  

The following segment of the code 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.

      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();
  

The rest of the function arranges the dofs in a sparsity pattern and initializes the system matrices and vectors.

      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());
    }
  

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().

    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());
    }
  

The following two constructors initialize scratch data from the input parameters and from another object of the same type, i.e., a copy constructor.

    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)
      , ve(0)
    {}
  
    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)
      , ve(0)
    {}
  

This 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, $A_{ij}$ , and the right-hand side, $b_i$ , by the function Solver::copy_local_to_global().

    void Solver::system_matrix_local(
      const typename DoFHandler<3>::active_cell_iterator &cell,
      AssemblyScratchData                                &scratch_data,
      AssemblyCopyData                                   &copy_data)
    {

First, we reinitialize the matrices and vectors related to the current cell and compute the FE values.

      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);

Second, we compute the free-current density, $\vec{J}_f$ , at the quadrature points.

      scratch_data.Jf.value_list(scratch_data.fe_values.get_quadrature_points(),
                                 cell->material_id(),
                                 scratch_data.Jf_list);
  

Third, we compute the components of the cell matrix and cell right-hand side. The labels of the integrals are the same as in the introduction to this tutorial.

      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[scratch_data.ve].curl(
                       i, q_index) * // curl N_i
                       scratch_data.fe_values[scratch_data.ve].curl(
                         j, q_index)              // curl N_j
                     + scratch_data.eta_squared * // eta^2
                         scratch_data.fe_values[scratch_data.ve].value(
                           i, q_index) * // N_i
                         scratch_data.fe_values[scratch_data.ve].value(
                           j, q_index) // N_j
                     ) *
                    scratch_data.fe_values.JxW(q_index); // dV
                }
              copy_data.cell_rhs(i) += // Integral I_b3-1.
                (scratch_data.Jf_list[q_index] *
                 scratch_data.fe_values[scratch_data.ve].curl(i, q_index)) *
                scratch_data.fe_values.JxW(q_index); // J_f.(curl N_i)dV.
            }
        }
  

Finally, 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.

      cell->get_dof_indices(copy_data.local_dof_indices);
    }
  

This function copies the components of a cell matrix and a cell right-hand side into the system matrix, $A_{i,j}$ , and the system right-hand side, $b_i$ .

    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);
    }
  

This function solves the system of linear equations. If Settings::log_cg_convergence == true, the convergence data is saved into a file. In theory, a CG solver can solve an $m \times m$ system of linear equations in at most $m$ steps. In practice, it can take more steps to converge. The convergence of the algorithm depends on the spectral properties of the system matrix. The best case is if the eigenvalues form a compact cluster away from zero. In our case, however, the eigenvalues are spread in between zero and the maximal eigenvalue. Consequently, we expect a poor convergence and increase the maximal number of iteration steps by a factor of 10, i.e., 10*system_rhs.size(). The stopping condition is

$\‍[ |\boldsymbol{b} - \boldsymbol{A}\boldsymbol{c}| < 10^{-6} |\boldsymbol{b}|. \‍]$

As soon as we use constraints, we must not forget to distribute them.

    void Solver::solve()
    {
      SolverControl control(10 * system_rhs.size(),
                            1.0e-6 * system_rhs.l2_norm(),
                            false,
                            false);
  
      if (Settings::log_cg_convergence)
        control.enable_history_data();
  
      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);
  
      if (Settings::log_cg_convergence)
        {
          const std::vector<double> history_data = control.get_history_data();
  
          std::ofstream ofs(fname + "_cg_convergence.csv");
  
          for (unsigned int i = 1; i < history_data.size(); i++)
            ofs << i << ", " << history_data[i] << "\n";
        }
    }
  

This function saves the computed current vector potential into a vtu file.

    void Solver::save() 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);
    }
  
    void Solver::run()
    {
      {
        TimerOutput::Scope timer_section(timer, "Make mesh");
        make_mesh();
      }
      {
        TimerOutput::Scope timer_section(timer, "Setup");
        setup();
      }
      {
        TimerOutput::Scope timer_section(timer, "Assemble");
        assemble();
      }
      {
        TimerOutput::Scope timer_section(timer, "Solve");
        solve();
      }
      {
        TimerOutput::Scope timer_section(timer, "Save");
        save();
      }
      {
        TimerOutput::Scope timer_section(timer, "Clear");
        clear();
      }
    }
  } // namespace SolverT
  

Solver - A

This name space contains all the code related to the computation of the magnetic vector potential, $\vec{A}$ . The main difference between this solver and the solver for the current vector potential, $\vec{T}$ , is in how the information on the source is fed to respective solvers. The solver for $\vec{T}$ is fed data sampled from the analytical closed-form expression for $\vec{J}_f$ . The solver for $\vec{A}$ is fed a field function, i.e., a numerically computed current vector potential, $\vec{T}$ .

namespace SolverA

{

This class describes the permeability in the entire problem domain. The permeability is given by the definition of the problem, see the introduction.

    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);
      }
    };
  

This class describes the parameter $\gamma$ 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

$\‍[ \gamma = \dfrac{1}{\mu_0 r}. \‍]$

    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());
      }
    };
  

This class implements the solver that minimizes the functional $F(\vec{A})$ . The numerically computed current vector potential, $\vec{T}$ , 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 $\vec{T}$ has been computed. The reference to the mesh is passed via the input parameter triangulation_T.

      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 save() const; // Saves computed T 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;
    

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.

      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;
      TimerOutput        timer;
  

This time we have two dof handlers, dof_handler_T for $\vec{T}$ and dof_handler for $\vec{A}$ . 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.

      using IteratorTuple =
        std::tuple<typename DoFHandler<3>::active_cell_iterator,
                   typename DoFHandler<3>::active_cell_iterator>;
  
      using IteratorPair = SynchronousIterators<IteratorTuple>;
  

The program utilizes the WorkStream technology. The step-9 tutorial does a much better job of explaining the workings of WorkStream. Reading the "WorkStream paper", see the glossary, is recommended. The following structures and functions are related to WorkStream.

      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 FEValuesExtractors::Vector ve;
  
        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)
      , timer(std::cout,
              (Settings::print_time_tables) ? TimerOutput::summary :
                                              TimerOutput::never,
              TimerOutput::cpu_and_wall_times_grouped)
    {}
  

This function initializes the dofs, applies the Dirichlet boundary condition, and initializes the vectors and matrices.

    void Solver::setup()
    {
      dof_handler.reinit(triangulation_T);
      dof_handler.distribute_dofs(fe);
  

The following segment of the code 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 code does the constraining.

      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();
  

The rest of the function arranges the dofs in a sparsity pattern and initializes the system matrix and the system vectors.

      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());
    }
  

Formally, this 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.

    void Solver::assemble()
    {
      WorkStream::run(IteratorPair(IteratorTuple(dof_handler.begin_active(),
                                                 dof_handler_T.begin_active())),
                      IteratorPair(
                        IteratorTuple(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());
    }
  

The following two constructors initialize scratch data from the input parameters and from another object of the same type, i.e., a copy constructor.

    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)
      , ve(0)
      , 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)
      , ve(0)
      , 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)
    {}
  

This 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, $A_{ij}$ , and the right-hand side, $b_i$ , by the function Solver::copy_local_to_global().

    void Solver::system_matrix_local(const IteratorPair  &IP,
                                     AssemblyScratchData &scratch_data,
                                     AssemblyCopyData    &copy_data)
    {

First we reinitialize the matrices and vectors related to the current cell.

      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);

Second, we extract the cells from the pair. We extract them in the correct order, see above.

      auto cell   = std::get<0>(*IP);
      auto cell_T = std::get<1>(*IP);
  

Third, we compute the ordered FE values, the permeability, and the values of the current vector potential, $\vec{T}$ , on the cell.

      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);
  
      scratch_data.fe_values_T[scratch_data.ve].get_function_values(
        scratch_data.dofs_T, scratch_data.T_values);
  

Fourth, we compute the components of the cell matrix and cell right-hand side. The labels of the integrals are the same as in the introduction to this tutorial.

      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[scratch_data.ve].curl(
                       i, q_index) * // curl N_i
                       scratch_data.fe_values[scratch_data.ve].curl(
                         j, q_index)              // curl N_j
                     + scratch_data.eta_squared * // eta^2
                         scratch_data.fe_values[scratch_data.ve].value(
                           i, q_index) * // N_i
                         scratch_data.fe_values[scratch_data.ve].value(
                           j, q_index) // N_j
                     ) *
                    scratch_data.fe_values.JxW(q_index); // dV
                }
              copy_data.cell_rhs(i) += // Integral I_b3-1.
                (scratch_data.T_values[q_index] *
                 scratch_data.fe_values[scratch_data.ve].curl(i, q_index)) *
                scratch_data.fe_values.JxW(q_index); // T.(curl N_i)dV
            }
        }
  

If the Robin boundary condition (first-order ABC) is ordered, we compute an extra integral over the boundary.

        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[scratch_data.ve]
                                       .value(i, q_index_face)) *
                                   cross_product_3d(
                                     scratch_data.fe_face_values.normal_vector(
                                       q_index_face),
                                     scratch_data.fe_face_values[scratch_data.ve]
                                       .value(j,
                                              q_index_face))) // (n x N_i).(n x N_j)
                                  * scratch_data.fe_face_values.JxW(
                                      q_index_face); // dS
                              }                      // for j
                          }                          // for i
                      }                              // for q_index_face
                  } // if (cell->face(f)->at_boundary())
              }     // for f
          }         /* if (scratch_data.boundary_condition_type ==
*                      BoundaryConditionType::Robin) */
    

Finally, 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.

      cell->get_dof_indices(copy_data.local_dof_indices);
    }
  

This function copies the components of a cell matrix and a cell right-hand side into the system matrix, $A_{i,j}$ , and the system right-hand side, $b_i$ .

    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);
    }
  

This function solves the system of linear equations. If Settings::log_cg_convergence == true, the convergence data is saved into a file. In theory, a CG solver can solve an $m \times m$ system of linear equations in at most $m$ steps. In practice, it can take more steps to converge. The convergence of the algorithm depends on the spectral properties of the system matrix. The best case is if the eigenvalues form a compact cluster away from zero. In our case, however, the eigenvalues are spread in between zero and the maximal eigenvalue. Consequently, we expect a poor convergence and increase the maximal number of iteration steps by a factor of 10, i.e., 10*system_rhs.size(). The stopping condition is

$\‍[ |\boldsymbol{b} - \boldsymbol{A}\boldsymbol{c}| < 10^{-6} |\boldsymbol{b}|. \‍]$

As soon as we use constraints, we must not forget to distribute them.

    void Solver::solve()
    {
      SolverControl control(10 * system_rhs.size(),
                            1.0e-6 * system_rhs.l2_norm(),
                            false,
                            false);
  
      if (Settings::log_cg_convergence)
        control.enable_history_data();
  
      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);
  
      if (Settings::log_cg_convergence)
        {
          const std::vector<double> history_data = control.get_history_data();
  
          std::ofstream ofs(fname + "_cg_convergence.csv");
  
          for (unsigned int i = 1; i < history_data.size(); i++)
            ofs << i << ", " << history_data[i] << "\n";
        }
    }
  

This function saves the computed magnetic vector potential into a vtu file.

    void Solver::save() 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);
    }
  
    void Solver::run()
    {
      {
        TimerOutput::Scope timer_section(timer, "Setup");
        setup();
      }
      {
        TimerOutput::Scope timer_section(timer, "Assemble");
        assemble();
      }
      {
        TimerOutput::Scope timer_section(timer, "Solve");
        solve();
      }
      {
        TimerOutput::Scope timer_section(timer, "Save");
        save();
      }
      {
        TimerOutput::Scope timer_section(timer, "Clear");
        clear();
      }
    }
  } // namespace SolverA
  
  

Projector from H(curl) to H(div)

This name space contains all the code related to the conversion of the magnetic vector potential, $\vec{A}$ , into magnetic field, $\vec{B}$ . 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, $\vec{T}$ into the free-current density, $\vec{J}_f$ .

  namespace ProjectorHcurlToHdiv
  {
  

This class implements the solver that minimizes the functional $F(\vec{B})$ or $F(\vec{J}_f)$ , see the introduction. The input vector field, $\vec{A}$ or $\vec{T}$ , 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.

      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 save() const;          // Saves computed T 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;
    

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.

      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;
      TimerOutput       timer;
  

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.

      using IteratorTuple =
        std::tuple<typename DoFHandler<3>::active_cell_iterator,
                   typename DoFHandler<3>::active_cell_iterator>;
  
      using IteratorPair = SynchronousIterators<IteratorTuple>;
  

The program utilizes the WorkStream technology. The step-9 tutorial does a much better job of explaining the workings of WarkStream. Reading the "WorkStream paper", see the glossary, is recommended. The following structures and functions are related to WorkStream.

      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 FEValuesExtractors::Vector ve;
  
        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)
      , timer(std::cout,
              (Settings::print_time_tables) ? TimerOutput::summary :
                                              TimerOutput::never,
              TimerOutput::cpu_and_wall_times_grouped)
    {
      Assert(exact_solution != nullptr,
             ExcMessage("The exact solution is missing."));
    }
  

This function initializes the dofs, vectors and matrices. This time there are no constraints as we do not apply Dirichlet boundary condition.

    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());
    }
  

Formally, this 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.

    void Solver::assemble()
    {
      WorkStream::run(IteratorPair(
                        IteratorTuple(dof_handler_Hdiv.begin_active(),
                                      dof_handler_Hcurl.begin_active())),
                      IteratorPair(IteratorTuple(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());
    }
  

The following two constructors initialize scratch data from the input parameters and from another object of the same type, i.e., a copy constructor.

    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)
      , ve(0)
      , 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)
      , ve(0)
      , dof_hand_Hcurl(scratch_data.dof_hand_Hcurl)
      , dofs_Hcurl(scratch_data.dofs_Hcurl)
    {}
  

This 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, $A_{ij}$ , and the right-hand side, $b_i$ , by the function Solver::copy_local_to_global().

    void Solver::system_matrix_local(const IteratorPair  &IP,
                                     AssemblyScratchData &scratch_data,
                                     AssemblyCopyData    &copy_data)
    {

First we reinitialize the matrices and vectors related to the current cell, update the FE values, and compute the curl of the input vector field.

      copy_data.cell_matrix.reinit(scratch_data.dofs_per_cell,
                                   scratch_data.dofs_per_cell);
  
      copy_data.cell_rhs.reinit(scratch_data.dofs_per_cell);
  
      copy_data.local_dof_indices.resize(scratch_data.dofs_per_cell);
  
      scratch_data.fe_values_Hdiv.reinit(std::get<0>(*IP));
      scratch_data.fe_values_Hcurl.reinit(std::get<1>(*IP));
  

The variable curl_vec_in_Hcurl denotes the curl of the input vector field, $\vec{\nabla} \times \vec{T}$ or $\vec{\nabla} \times \vec{A}$ , depending on the context.

      scratch_data.fe_values_Hcurl[scratch_data.ve].get_function_curls(
        scratch_data.dofs_Hcurl, scratch_data.curl_vec_in_Hcurl);
  

Second, we compute the components of the cell matrix and cell right-hand side. The labels of the integrals are the same as in the introduction to this tutorial.

      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[scratch_data.ve].value(i,
                                                                       q_index) *
                    scratch_data.fe_values_Hdiv[scratch_data.ve].value(j,
                                                                       q_index) *
                    scratch_data.fe_values_Hdiv.JxW(q_index);
                }
  
              copy_data.cell_rhs(i) += // Integral I_b
                scratch_data.curl_vec_in_Hcurl[q_index] *
                scratch_data.fe_values_Hdiv[scratch_data.ve].value(i, q_index) *
                scratch_data.fe_values_Hdiv.JxW(q_index);
            }
        }
  

Finally, 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.

      std::get<0>(*IP)->get_dof_indices(copy_data.local_dof_indices);
    }
  

This function copies the components of a cell matrix and a cell right-hand side into the system matrix, $A_{i,j}$ , and the system right-hand side, $b_i$ .

    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);
    }
  

The following two functions compute the error norms and project the exact solution.

    void Solver::compute_error_norms()
    {
      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()
    {
      AffineConstraints<double> 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);
    }
  

This function solves the system of linear equations. This time we are dealing with a mass matrix. It has good spectral properties. Consequently, we do not use the factor of 10 as in preceding two solvers. The stopping condition is

$\‍[ |\boldsymbol{b} - \boldsymbol{A}\boldsymbol{c}| < 10^{-6} |\boldsymbol{b}|. \‍]$

    void Solver::solve()
    {
      SolverControl control(system_rhs.size(),
                            1.0e-6 * system_rhs.l2_norm(),
                            false,
                            false);
  
      if (Settings::log_cg_convergence)
        control.enable_history_data();
  
      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);
  
      if (Settings::log_cg_convergence)
        {
          const std::vector<double> history_data = control.get_history_data();
  
          std::ofstream ofs(fname + "_cg_convergence.csv");
  
          for (unsigned int i = 1; i < history_data.size(); i++)
            ofs << i << ", " << history_data[i] << "\n";
        }
    }
  

This function saves the computed fields into a vtu file. This time we also save the projected exact solution and the $L^2$ error norm. The exact solution is only saved if the Settings::project_exact_solution = true

    void Solver::save() 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);
    }
  
    void Solver::run()
    {
      {
        TimerOutput::Scope timer_section(timer, "Setup");
        setup();
      }
      {
        TimerOutput::Scope timer_section(timer, "Assemble");
        assemble();
      }
      {
        TimerOutput::Scope timer_section(timer, "Solve");
        solve();
      }
  
      if (exact_solution)
        {
          {
            TimerOutput::Scope timer_section(timer, "Compute error norms");
            compute_error_norms();
          }
  
          if (Settings::project_exact_solution)
            {
              {
                TimerOutput::Scope timer_section(timer, "Project exact solution");
                project_exact_solution_fcn();
              }
            }
        }
      {
        TimerOutput::Scope timer_section(timer, "Save");
        save();
      }
      {
        TimerOutput::Scope timer_section(timer, "Clear");
        clear();
      }
    }
  } // namespace ProjectorHcurlToHdiv
  
  

The main loop

This class contains the main loop of the program.

  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);
  
      table_Jf.clear();
      table_B.clear();
  
      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);
  

Stage 1. Computing $\vec{T}$ .

          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();
  

Stage 2. Computing $\vec{J}_f$ .

          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());
  

Stage 3. Computing $\vec{A}$ .

          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();
  

Stage 4. Computing $\vec{B}$ .

          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());

End stage 4.

          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;
  }

Results

The program generates the following output in the command line interface by default.

Solving for (p = 0): T Jf A B T Jf A B T Jf A B T Jf A B

The program assumes the finite elements of the lowermost degree, $p = 0$ . To change the degree of the finite elements, say $p = 2$ , 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 Notes on visualizing high order output 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, $p = 0, 1, 2$ .

p	r	cells	dofs	$\\|e\\|_{L^2}$	$\alpha_{L^2}$
0	6	4625	13950	1.66e-01	-
0	7	7992	24084	1.38e-01	0.99
0	8	12691	38220	1.19e-01	0.99
0	9	18944	57024	1.04e-01	0.99
1	6	4625	111300	8.12e-04	-
1	7	7992	192240	4.97e-04	2.69
1	8	12691	305172	3.32e-04	2.61
1	9	18944	455424	2.37e-04	2.54
2	6	4625	375300	6.78e-04	-
2	7	7992	648324	3.94e-04	2.97
2	8	12691	1029294	2.49e-04	2.98
2	9	18944	1536192	1.67e-04	2.99

Table 1. Convergence table. Free-current density, $\vec{J}_f$ .

p	r	cells	dofs	$\\|e\\|_{L^2}$	$\alpha_{L^2}$
0	6	4625	13950	8.84e-08	-
0	7	7992	24084	7.36e-08	1.00
0	8	12691	38220	6.30e-08	1.01
0	9	18944	57024	5.51e-08	1.00
1	6	4625	111300	4.41e-09	-
1	7	7992	192240	3.11e-09	1.91
1	8	12691	305172	2.23e-09	2.18
1	9	18944	455424	1.71e-09	1.96
2	6	4625	375300	1.84e-10	-
2	7	7992	648324	1.03e-10	3.21
2	8	12691	1029294	6.08e-11	3.40
2	9	18944	1536192	4.04e-11	3.07

Table 2. Convergence table. Magnetic field, $\vec{B}$ .

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.
$\|e\|_{L^2}$ - the error norm.
$\alpha_{L^2}$ - the order of convergence of the error norm.

If Settings::log_cg_convergence = true, the program saves the convergence data of the CG solver into csv files.

The vector representations of the calculated vector fields, $\vec{J}_f$ and $\vec{B}$ , 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 $p = 2$ and $r = 9$ . Visual inspection of the vector potentials is not very informative as their conservative portions are unknown.

The result - free-current
density

The result - magnetic field
density

Possibilities for extensions

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 $\eta^2$ 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., $p = 1$ and $p = 2$ .

The Robin boundary condition as described above is also called the first-order asymptotic boundary condition (ABC). There exist ABCs of higher orders [gratkowski2010p]. Implement and test the second-order ABC to see if it performs any better. There exist improvised asymptotic boundary conditions, IABCs, [meeker2013a]. Try to implement the first order IABC.

The plain program

/* ------------------------------------------------------------------------
 *
 * SPDX-License-Identifier: LGPL-2.1-or-later
 * Copyright (C) 2024 - 2025 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * Part of the source code is dual licensed under Apache-2.0 WITH
 * LLVM-exception OR LGPL-2.1-or-later. Detailed license information
 * governing the source code and code contributions can be found in
 * LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
 *
 * ------------------------------------------------------------------------
 *
 * This program was contributed by Siarhei Uzunbajakau, www.cembooks.nl, 2025.
 */
 
#include <deal.II/base/tensor.h>
#include <deal.II/base/function.h>
#include <deal.II/base/convergence_table.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/types.h>
#include <deal.II/base/work_stream.h>
#include <deal.II/base/multithread_info.h>
 
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/manifold_lib.h>
 
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
 
#include <deal.II/fe/fe_nedelec.h>
#include <deal.II/fe/fe_raviart_thomas.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/fe/fe_values.h>
 
#include <deal.II/lac/vector.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/vector_memory.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
 
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/data_component_interpretation.h>
#include <deal.II/numerics/data_out.h>
 
#include <string>
#include <ostream>
 
using namespace dealii;
 
enum BoundaryConditionType
{
  Dirichlet,
  Neumann,
  Robin
};
 
 
namespace Settings
{
  const double permeability_fs = 1.2566370614359172954e-6;
  const double mu_0            = permeability_fs; // Permeability of free space.
 
  const double mu_r = 4; // Relative permeability of the magnetic material.
  const double mu_1 = mu_0 * mu_r; // Permeability of the magnetic material.
  const double d1   = 0.1; // Half-side of the cube in the center of the mesh.
  const double a1   = 0.3; // Inner radius of the magnetic core.
  const double b1   = 0.6; // Outer radius of the magnetic core.
  const double a2   = 0.9; // Inner radius of the free-current region.
  const double b2   = 1.2; // Outer radius of the free-current region.
  const double d2   = 2.0; // Error norms are computed if r < d2.
  const double K0   = 1.0; // Magnitude of the free-current density.
  const double H0   = // Magnitude of H-field inside the coil, core removed.
    (1.0 / 3.0) * K0 * (std::pow(b2, 2) - std::pow(a2, 2));
  const types::material_id material_id_free_space =
    1; // Material ID of the free space.
  const types::material_id material_id_core =
    2; // Material ID of the magnetic core.
  const types::material_id material_id_free_current =
    3; // Material ID of the free-current region.
  const types::boundary_id boundary_id_infinity =
    2; // Boundary ID. Set in the geo file.
 
  const unsigned int mapping_degree = 2; // Mapping degree used in all solvers.
  const unsigned int fe_degree      = 0; // Degree of the finite elements.
  const BoundaryConditionType boundary_condition_type_A = Robin;
 
  const double eta_squared_T = 0.0; // eta^2 when solving for T@f$.
  const double eta_squared_A = 0.0; // eta^2 when solving for A@f$.
 
  const unsigned int n_threads_max = 0; // If >0 limits the number of threads.
  const double eps = 1e-12; // Two doubles are equal if their difference < eps.
  const bool   log_cg_convergence = false; // Save the CG convergence data.
  const bool print_time_tables = false; // Print the time tables on the screen.
  const bool project_exact_solution = false; // Save the exact solution.
} // namespace Settings
 
class Weight : public Function<3>
{
public:
  virtual double value(const Point<3> &p,
                       const unsigned int) const override final
  {
    if (p.norm() > Settings::d2)
      return 0.0;
 
    return 1.0;
  }
};
 
class MainOutputTable : public ConvergenceTable
{
public:
  MainOutputTable() = delete;
 
  MainOutputTable(const unsigned int dimensions)
    : ConvergenceTable()
    , dimensions(dimensions)
  {}
 
  void set_new_order(const std::vector<std::string> &new_order_in)
  {
    new_order = new_order_in;
  }
 
  void append_new_order(const std::string &new_column)
  {
    new_order.push_back(new_column);
  }
 
  void format()
  {
    set_precision("L2", 2);
 
    set_scientific("L2", true);
 
    evaluate_convergence_rates("L2",
                               "ncells",
                               ConvergenceTable::reduction_rate_log2,
                               dimensions);
 
    set_tex_caption("p", "p");
    set_tex_caption("r", "r");
    set_tex_caption("ncells", "nr. cells");
    set_tex_caption("ndofs", "nr. dofs");
    set_tex_caption("L2", "L2 norm");
 
    set_column_order(new_order);
  }
 
  void save(const std::string &fname)
  {
    format();
 
    std::ofstream ofs(fname + ".tex");
    write_tex(ofs);
  }
 
private:
  const unsigned int       dimensions;
  std::vector<std::string> new_order = {"p", "r", "ncells", "ndofs", "L2"};
};
 
 
namespace ExactSolutions
{
 
  inline Tensor<1, 3> volume_free_current_density(const Point<3> &p,
                                                  const double    K0)
  {
    return Tensor<1, 3>({-K0 * p[1], K0 * p[0], 0.0});
  }
 
  inline Tensor<1, 3> magnetic_field_coil(const Point<3> &p,
                                          const double    K0,
                                          const double    mu_0,
                                          const double    a2,
                                          const double    b2)
  {
    const double r = p.norm();
 
    const double cos_theta = p[2] / r;
    const double sin_theta =
      std::sqrt(std::pow(p[0], 2) + std::pow(p[1], 2)) / r;
 
    const double cos_phi =
      p[0] / std::sqrt(std::pow(p[0], 2) + std::pow(p[1], 2));
    const double sin_phi =
      p[1] / std::sqrt(std::pow(p[0], 2) + std::pow(p[1], 2));
 
    const Tensor<1, 3> r_hat({p[0] / r, p[1] / r, p[2] / r});
    const Tensor<1, 3> theta_hat(
      {cos_theta * cos_phi, cos_theta * sin_phi, -sin_theta});
 
    const Tensor<1, 3> F1 =
      (2.0 / 3.0) * mu_0 * K0 * (cos_theta * r_hat - sin_theta * theta_hat);
    const Tensor<1, 3> F2 = (2.0 / 3.0) * mu_0 * K0 *
                            (cos_theta * r_hat + 0.5 * sin_theta * theta_hat) /
                            std::pow(r, 3);
 
    if (r <= a2)
      {
        return 0.5 * (std::pow(b2, 2) - std::pow(a2, 2)) * F1;
      }
    else if (r >= b2)
      {
        return 0.2 * (std::pow(b2, 5) - std::pow(a2, 5)) * F2;
      }
    else
      {
        return 0.5 * (std::pow(b2, 2) - std::pow(r, 2)) * F1 +
               0.2 * (std::pow(r, 5) - std::pow(a2, 5)) * F2;
      }
 
    return Tensor<1, 3>({0.404, 0.404, 0.404});
  }
 
  inline Tensor<1, 3> magnetic_field_core(const Point<3> &p,
                                          const double    H0,
                                          const double    mur,
                                          const double    mu0,
                                          const double    a1,
                                          const double    b1)
  {
    const double a3 = std::pow(a1, 3);
    const double b3 = std::pow(b1, 3);
 
    const double OMEGA   = ((mur - 1.0) / (mur + 2.0)) * (a3 / b3);
    const double gamma_1 = (-3.0 * b3 * H0 * OMEGA) /
                           ((2.0 * mur + 1.0) - 2.0 * (mur - 1.0) * OMEGA);
    const double beta_1 = ((2.0 * mur + 1.0) * gamma_1) / ((mur - 1.0) * a3);
    const double alpha_1 =
      (-b3 * H0 + 2.0 * mur * gamma_1 - mur * b3 * beta_1) / 2.0;
    const double delta_1 = (mur * a3 * beta_1 - 2.0 * mur * gamma_1) / a3;
 
    const double r  = p.norm();
    const double r3 = std::pow(r, 3);
    const double r5 = std::pow(r, 5);
 
    const double zz = -3.0 * p[2] * p[2] / r5;
    const double xz = -3.0 * p[0] * p[2] / r5;
    const double yz = -3.0 * p[1] * p[2] / r5;
 
    if (r <= a1)
      {
        return -mu0 * (Tensor<1, 3>({0.0, 0.0, delta_1}) +
                       Tensor<1, 3>({0.0, 0.0, H0}));
      }
    else if (r >= b1)
      {
        return -mu0 * (Tensor<1, 3>({alpha_1 * xz,
                                     alpha_1 * yz,
                                     -H0 + alpha_1 / r3 + alpha_1 * zz}) +
                       Tensor<1, 3>({0.0, 0.0, H0}));
      }
    else
      {
        return -mu0 *
               (mur * Tensor<1, 3>({gamma_1 * xz,
                                    gamma_1 * yz,
                                    beta_1 + gamma_1 / r3 + gamma_1 * zz}) +
                Tensor<1, 3>({0.0, 0.0, H0}));
      }
 
    return Tensor<1, 3>({0.404, 0.404, 0.404});
  }
 
  class FreeCurrentDensity : public Function<3>
  {
  public:
    FreeCurrentDensity()
      : Function<3>(3)
    {}
 
    virtual void
    vector_value_list(const std::vector<Point<3>> &p,
                      std::vector<Vector<double>> &values) const override final
    {
      Assert(values.size() == p.size(),
             ExcDimensionMismatch(values.size(), p.size()));
 
      double r;
 
      for (unsigned int i = 0; i < values.size(); i++)
        {
          r = p[i].norm();
 
          if ((r >= Settings::a2) && (r <= Settings::b2))
            {
              const Tensor<1, 3> Jf =
                volume_free_current_density(p[i], Settings::K0);
 
              values[i][0] = Jf[0];
              values[i][1] = Jf[1];
              values[i][2] = Jf[2];
            }
          else
            {
              values[i][0] = 0.0;
              values[i][1] = 0.0;
              values[i][2] = 0.0;
            }
        }
    }
  };
 
  class MagneticField : public Function<3>
  {
  public:
    MagneticField()
      : Function<3>(3)
    {}
 
    virtual void
    vector_value_list(const std::vector<Point<3>> &p,
                      std::vector<Vector<double>> &values) const override final
    {
      Assert(values.size() == p.size(),
             ExcDimensionMismatch(values.size(), p.size()));
 
      Tensor<1, 3> B;
 
      for (unsigned int i = 0; i < values.size(); i++)
        {
          using namespace Settings;
          B = magnetic_field_coil(p[i], K0, mu_0, a2, b2) +
              magnetic_field_core(p[i], H0, mu_r, mu_0, a1, b1);
 
          values[i][0] = B[0];
          values[i][1] = B[1];
          values[i][2] = B[2];
        }
    }
  };
 
} // namespace ExactSolutions
 
 
namespace SolverT
{
  class FreeCurrentDensity
  {
  public:
    void value_list(const std::vector<Point<3>> &p, // Quadrature points.
                    const types::material_id   mid, // Material ID of the cell.
                    std::vector<Tensor<1, 3>> &values) const
    {
      Assert(p.size() == values.size(),
             ExcDimensionMismatch(p.size(), values.size()));
 
      if ((mid == Settings::material_id_free_space) ||
          (mid == Settings::material_id_core))
        for (unsigned int i = 0; i < values.size(); i++)
          values[i] = Tensor<1, 3>({0.0, 0.0, 0.0});
 
      if (mid == Settings::material_id_free_current)
        for (unsigned int i = 0; i < values.size(); i++)
          values[i] =
            ExactSolutions::volume_free_current_density(p[i], Settings::K0);
    }
  };
 
  class Solver
  {
  public:
    Solver() = delete;
    Solver(const unsigned int p, // Degree of the Nedelec finite elements.
           const unsigned int r, // The mesh refinement parameter.
           const unsigned int mapping_degree,
           const double       eta_squared = 0.0,
           const std::string &fname       = "data");
 
    void make_mesh();  // Loads the mesh. Assigns mat. IDs. Attaches manifold.
    void setup();      // Initializes dofs, vectors, matrices.
    void assemble();   // Assembles the system of linear equations.
    void solve();      // Solves the system of linear equations.
    void save() const; // Saves computed T 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 six functions in the proper order
                   and measures the execution time for each function. */
 
    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:
 
    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;
    TimerOutput        timer;
 
    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;
      const FEValuesExtractors::Vector ve;
    };
 
    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)
    , timer(std::cout,
            (Settings::print_time_tables) ? TimerOutput::summary :
                                            TimerOutput::never,
            TimerOutput::cpu_and_wall_times_grouped)
  {}
 
  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);
  }
 
  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());
  }
 
  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());
  }
 
  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)
    , ve(0)
  {}
 
  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)
    , ve(0)
  {}
 
  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);
 
    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[scratch_data.ve].curl(
                     i, q_index) * // curl N_i
                     scratch_data.fe_values[scratch_data.ve].curl(
                       j, q_index)              // curl N_j
                   + scratch_data.eta_squared * // eta^2
                       scratch_data.fe_values[scratch_data.ve].value(
                         i, q_index) * // N_i
                       scratch_data.fe_values[scratch_data.ve].value(
                         j, q_index) // N_j
                   ) *
                  scratch_data.fe_values.JxW(q_index); // dV
              }
            copy_data.cell_rhs(i) += // Integral I_b3-1.
              (scratch_data.Jf_list[q_index] *
               scratch_data.fe_values[scratch_data.ve].curl(i, q_index)) *
              scratch_data.fe_values.JxW(q_index); // J_f.(curl N_i)dV.
          }
      }
 
    cell->get_dof_indices(copy_data.local_dof_indices);
  }
 
  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);
  }
 
  void Solver::solve()
  {
    SolverControl control(10 * system_rhs.size(),
                          1.0e-6 * system_rhs.l2_norm(),
                          false,
                          false);
 
    if (Settings::log_cg_convergence)
      control.enable_history_data();
 
    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);
 
    if (Settings::log_cg_convergence)
      {
        const std::vector<double> history_data = control.get_history_data();
 
        std::ofstream ofs(fname + "_cg_convergence.csv");
 
        for (unsigned int i = 1; i < history_data.size(); i++)
          ofs << i << ", " << history_data[i] << "\n";
      }
  }
 
  void Solver::save() 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);
  }
 
  void Solver::run()
  {
    {
      TimerOutput::Scope timer_section(timer, "Make mesh");
      make_mesh();
    }
    {
      TimerOutput::Scope timer_section(timer, "Setup");
      setup();
    }
    {
      TimerOutput::Scope timer_section(timer, "Assemble");
      assemble();
    }
    {
      TimerOutput::Scope timer_section(timer, "Solve");
      solve();
    }
    {
      TimerOutput::Scope timer_section(timer, "Save");
      save();
    }
    {
      TimerOutput::Scope timer_section(timer, "Clear");
      clear();
    }
  }
} // namespace SolverT
 
 
namespace SolverA
{
  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);
    }
  };
 
  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());
    }
  };
 
  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 save() const; // Saves computed T 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;
 
    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;
    TimerOutput        timer;
 
    using IteratorTuple =
      std::tuple<typename DoFHandler<3>::active_cell_iterator,
                 typename DoFHandler<3>::active_cell_iterator>;
 
    using IteratorPair = SynchronousIterators<IteratorTuple>;
 
    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 FEValuesExtractors::Vector ve;
 
      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)
    , timer(std::cout,
            (Settings::print_time_tables) ? TimerOutput::summary :
                                            TimerOutput::never,
            TimerOutput::cpu_and_wall_times_grouped)
  {}
 
  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());
  }
 
  void Solver::assemble()
  {
    WorkStream::run(IteratorPair(IteratorTuple(dof_handler.begin_active(),
                                               dof_handler_T.begin_active())),
                    IteratorPair(
                      IteratorTuple(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());
  }
 
  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)
    , ve(0)
    , 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)
    , ve(0)
    , 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)
  {}
 
  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);
 
    auto cell   = std::get<0>(*IP);
    auto 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);
 
    scratch_data.fe_values_T[scratch_data.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[scratch_data.ve].curl(
                     i, q_index) * // curl N_i
                     scratch_data.fe_values[scratch_data.ve].curl(
                       j, q_index)              // curl N_j
                   + scratch_data.eta_squared * // eta^2
                       scratch_data.fe_values[scratch_data.ve].value(
                         i, q_index) * // N_i
                       scratch_data.fe_values[scratch_data.ve].value(
                         j, q_index) // N_j
                   ) *
                  scratch_data.fe_values.JxW(q_index); // dV
              }
            copy_data.cell_rhs(i) += // Integral I_b3-1.
              (scratch_data.T_values[q_index] *
               scratch_data.fe_values[scratch_data.ve].curl(i, q_index)) *
              scratch_data.fe_values.JxW(q_index); // T.(curl N_i)dV
          }
      }
 
    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[scratch_data.ve]
                                   .value(i, q_index_face)) *
                               cross_product_3d(
                                 scratch_data.fe_face_values.normal_vector(
                                   q_index_face),
                                 scratch_data.fe_face_values[scratch_data.ve]
                                   .value(j,
                                          q_index_face))) // (n x N_i).(n x N_j)
                              * scratch_data.fe_face_values.JxW(
                                  q_index_face); // dS
                          }                      // 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);
  }
 
  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);
  }
 
  void Solver::solve()
  {
    SolverControl control(10 * system_rhs.size(),
                          1.0e-6 * system_rhs.l2_norm(),
                          false,
                          false);
 
    if (Settings::log_cg_convergence)
      control.enable_history_data();
 
    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);
 
    if (Settings::log_cg_convergence)
      {
        const std::vector<double> history_data = control.get_history_data();
 
        std::ofstream ofs(fname + "_cg_convergence.csv");
 
        for (unsigned int i = 1; i < history_data.size(); i++)
          ofs << i << ", " << history_data[i] << "\n";
      }
  }
 
  void Solver::save() 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);
  }
 
  void Solver::run()
  {
    {
      TimerOutput::Scope timer_section(timer, "Setup");
      setup();
    }
    {
      TimerOutput::Scope timer_section(timer, "Assemble");
      assemble();
    }
    {
      TimerOutput::Scope timer_section(timer, "Solve");
      solve();
    }
    {
      TimerOutput::Scope timer_section(timer, "Save");
      save();
    }
    {
      TimerOutput::Scope timer_section(timer, "Clear");
      clear();
    }
  }
} // namespace SolverA
 
 
namespace ProjectorHcurlToHdiv
{
 
  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 save() const;          // Saves computed T 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;
 
    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;
    TimerOutput       timer;
 
    using IteratorTuple =
      std::tuple<typename DoFHandler<3>::active_cell_iterator,
                 typename DoFHandler<3>::active_cell_iterator>;
 
    using IteratorPair = SynchronousIterators<IteratorTuple>;
 
    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 FEValuesExtractors::Vector ve;
 
      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)
    , timer(std::cout,
            (Settings::print_time_tables) ? TimerOutput::summary :
                                            TimerOutput::never,
            TimerOutput::cpu_and_wall_times_grouped)
  {
    Assert(exact_solution != nullptr,
           ExcMessage("The exact solution is missing."));
  }
 
  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());
  }
 
  void Solver::assemble()
  {
    WorkStream::run(IteratorPair(
                      IteratorTuple(dof_handler_Hdiv.begin_active(),
                                    dof_handler_Hcurl.begin_active())),
                    IteratorPair(IteratorTuple(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());
  }
 
  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)
    , ve(0)
    , 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)
    , ve(0)
    , dof_hand_Hcurl(scratch_data.dof_hand_Hcurl)
    , dofs_Hcurl(scratch_data.dofs_Hcurl)
  {}
 
  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);
 
    scratch_data.fe_values_Hdiv.reinit(std::get<0>(*IP));
    scratch_data.fe_values_Hcurl.reinit(std::get<1>(*IP));
 
    scratch_data.fe_values_Hcurl[scratch_data.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[scratch_data.ve].value(i,
                                                                     q_index) *
                  scratch_data.fe_values_Hdiv[scratch_data.ve].value(j,
                                                                     q_index) *
                  scratch_data.fe_values_Hdiv.JxW(q_index);
              }
 
            copy_data.cell_rhs(i) += // Integral I_b
              scratch_data.curl_vec_in_Hcurl[q_index] *
              scratch_data.fe_values_Hdiv[scratch_data.ve].value(i, q_index) *
              scratch_data.fe_values_Hdiv.JxW(q_index);
          }
      }
 
    std::get<0>(*IP)->get_dof_indices(copy_data.local_dof_indices);
  }
 
  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);
  }
 
  void Solver::compute_error_norms()
  {
    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()
  {
    AffineConstraints<double> 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);
  }
 
  void Solver::solve()
  {
    SolverControl control(system_rhs.size(),
                          1.0e-6 * system_rhs.l2_norm(),
                          false,
                          false);
 
    if (Settings::log_cg_convergence)
      control.enable_history_data();
 
    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);
 
    if (Settings::log_cg_convergence)
      {
        const std::vector<double> history_data = control.get_history_data();
 
        std::ofstream ofs(fname + "_cg_convergence.csv");
 
        for (unsigned int i = 1; i < history_data.size(); i++)
          ofs << i << ", " << history_data[i] << "\n";
      }
  }
 
  void Solver::save() 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);
  }
 
  void Solver::run()
  {
    {
      TimerOutput::Scope timer_section(timer, "Setup");
      setup();
    }
    {
      TimerOutput::Scope timer_section(timer, "Assemble");
      assemble();
    }
    {
      TimerOutput::Scope timer_section(timer, "Solve");
      solve();
    }
 
    if (exact_solution)
      {
        {
          TimerOutput::Scope timer_section(timer, "Compute error norms");
          compute_error_norms();
        }
 
        if (Settings::project_exact_solution)
          {
            {
              TimerOutput::Scope timer_section(timer, "Project exact solution");
              project_exact_solution_fcn();
            }
          }
      }
    {
      TimerOutput::Scope timer_section(timer, "Save");
      save();
    }
    {
      TimerOutput::Scope timer_section(timer, "Clear");
      clear();
    }
  }
} // namespace ProjectorHcurlToHdiv
 
 
 
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);
 
    table_Jf.clear();
    table_B.clear();
 
    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);
 
 
        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();
 
 
        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());
 
 
        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();
 
 
        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());
 
        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;
}