Some problems in magnetostatics can be formulated in terms of the magnetic vector potential, \(\vec{A}\). There are four options:

The magnetic vector potential has the following form (Section 3.3 Magnetic vector potential in three dimensions)
\begin{equation} \vec{A} = A_x(x,y,z) \hat{i} + A_y(x,y,z) \hat{j} + A_z(x,y,z) \hat{k}. \end{equation}
in the three-dimensional Cartesian coordinate system. It cannot be simplified as the three-dimensional problem domain does not exhibit a translation or rotation symmetry. In this case the problem can be described by the three-dimensional static vector boundary value problem shown below on this page.
The magnetic vector potential has the following form (Section 3.4.3 Vector planar problem)
\begin{equation} \vec{A} = A_x(x,y)\hat{i} + A_y(x,y)\hat{j} + 0\hat{k}, \end{equation}
in the three-dimensional Cartesian coordinate system. Such problem exhibits a translation symmetry and the magnetic vector potential curls in planes that are perpendicular to the axis of symmetry. In this case the problem can be described by the two-dimensional static vector boundary value problem shown below on this page.
The magnetic vector potential has the following form (Section 3.4.1 Scalar planar problem)
\begin{equation} \vec{A} = 0 \hat{i} + 0 \hat{j} + A_z(x,y) \hat{k}, \end{equation}
in the three-dimensional Cartesian coordinate system. Such problem exhibits a translation symmetry and the magnetic vector potential runs parallel the axis of symmetry. In this case the problem can be reformulated in terms of the magnitude of the magnetic vector potential \(A\) and described by the two-dimensional static scalar boundary value problem. Note that the specific form of the magnetic vector potential implicitly gauges the vector potential. Therefore, the unique solution in terms of \(A\) is expected.
The magnetic vector potential has the following form (Section 3.4.2 Scalar axisymmetric problem)
\begin{equation} \vec{A} = 0 \hat{r} + A_{\phi}(r,z) \hat{\phi} + 0 \hat{z}, \end{equation}
in the cylindrical coordinate system. Such problem exhibits a rotation symmetry and the magnetic vector potential curls in planes that are perpendicular to the axis of symmetry. In this case the problem can be reformulated in terms of the scaled magnitude of the magnetic vector potential \(A'\) and described by the two-dimensional static scalar boundary value problem. Note that the specific form of the magnetic vector potential implicitly gauges the vector potential. Therefore, the unique solution in terms of \(A'\) is expected.
The magnetic vector potential has the following form (Section 3.4.4 Vector axisymmetric problem)
\begin{equation} \vec{A} = A_r(r,z) \hat{r} + 0 \hat{\phi} + A_z(r,z) \hat{z}, \end{equation}
in the cylindrical coordinate system. Such problem exhibits a rotation symmetry and the magnetic vector potential curls in planes that contain the axis of symmetry. In this case a curl-curl partial differential equation cannot be formulated.

That is, only the first two options lead to boundary value problems based on curl-curl partial differential equations. These two options are considered next.

The first option above leads to the following three-dimensional static vector boundary value problem:

\begin{equation} \begin{array}{lrcll} \text{ } & \vec{\nabla}\times\bigg(\dfrac{1}{\mu} \vec{\nabla}\times\vec{A}\bigg) + \eta^2 \vec{A} = \vec{J}_f & \text{in} & \Omega & \text{(i)},\\ \text{(e)}&\hat{n} \times \vec{A} = \hat{n} \times \vec{G} & \text{on} & \Gamma_{Dn} & \text{(ii)},\\ \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) = \vec{Q} & \text{on} & \Gamma_{Rm} & \text{(iii)}, \\ \text{(e)}&\hat{n}\times\vec{A}_{+} = \hat{n}\times\vec{A}_{-}&\text{on}&\Gamma_{Ik}&\text{(iv)},\\ \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) = \vec{K}_f & \text{on} & \Gamma_{Ik} & \text{(v)}. \end{array} \end{equation}

The following holds for this boundary value problem.

The Robin boundary condition (iii) becomes the Neumann boundary condition for \(\gamma=0\). In the case of the Robin boundary condition

\begin{equation} \begin{array}{rcll} \gamma > 0 & \text{on} & \Gamma_{Rm}. \end{array} \end{equation}
Specification of the Dirichlet boundary condition (ii), or the Robin boundary condition (equation (iii) with \(\gamma\ne 0\)) or the Neumann boundary condition (equation (iii) with \(\gamma = 0\)) will guarantee the uniqueness of the curl of the solution. The curl of the solution equals the magnetic field per definition:
\begin{equation} \vec{B} = \vec{\nabla} \times \vec{A} \end{equation}
The vector potential is not explicitly gauged. The gauging term, \(\eta^2 \vec{A}\), can be regarded as an implicit gauge. That is, it forces the partial differential operator to be positive definite so the conjugate gradient algorithm can find a solution. The conservative part of the solution, however, will remain unknown. The "unknown" means we cannot predict in advance which conservative component of the solution will show up at the output of the conjugate gradient solver. The addition of the gauging term disturbs the initial curl-curl equation, i.e., equation
\begin{equation} \vec{\nabla}\times\bigg(\frac{1}{\mu} \vec{\nabla}\times\vec{A}\bigg) = \vec{J}_f, \end{equation}
and disrupts the structure set by the De Rham complex (the gauging term is in \(H(\text{curl})\) while the other terms are in \(H(\text{div})\)). To switch off the gauging term set \(\eta = 0\). If used, the gauging parameter \( \eta \) must be chosen to be small to make sure that the gauging term is smaller than the other terms of the equation in hope that it will help the CG linear solver find a solution and will not disturb the solution too much.
Presence of the interfaces between the materials, \(\Gamma_{Ik}\), does not affect the uniqueness of the solution in any way. Equations (iv) and (v) describe the interface conditions.
The index "+" refers to the space immediately next to the interface in the direction of vector \(\hat{n}\). The index "-" refers to the space immediately next to the interface in the direction opposite to \(\hat{n}\).
The vector \(\hat{n}\) normal to a boundary always points outside the volume enclosed by the boundary. The same is true for vectors \(\hat{n}\) normal to interfaces.
The volume free-current density, \(\vec{J}_f\), does not include the surface free-current density, \( \vec{K}_f \). Thinking that \(\vec{J}_f = \vec{J}'_f + \vec{K}_f\) is wrong.
The solution to this boundary value problem minimizes the following functional (Chapter 4 Variational formulations):

\begin{equation} \begin{aligned} &F(\vec{A}) = \iiint_{\Omega}\frac{1}{\mu}\bigg|\vec{\nabla}\times\vec{A}\bigg|^2 dV + \sum_m \iint_{\Gamma_{Rm}} \bigg(\gamma\bigg|\hat{n}\times\vec{A}\bigg|^2 + 2\vec{Q}\cdot\vec{A}\bigg) dS + \eta^2 \iiint_{\Omega}\mid\vec{A}\mid^2 dV -2\iiint_{\Omega}\vec{J}_f \cdot \vec{A} dV - \\ & -2\sum_k\iint_{\Gamma_{Ik}} \vec{K}_f \cdot \vec{A} dS. \end{aligned} \end{equation}

The label (e) marks essential boundary or interface conditions. The functional \(F(\vec{A})\) is invariant to these conditions. These conditions must be imposed elsewhere (meaning any part of the system except the functional). For instance, the Dirichlet boundary condition (ii) is imposed by constraining the system of linear equations. The continuity equation (iv) is imposed by the choice of the finite elements. The FE_Nedelec finite elements of deal.II, for instance, do not allow discontinuity of the tangential component on the boundaries of the cells. Therefore, the choice of these elements makes a violation of the continuity condition (iv) utterly impossible.
The label (n) marks natural boundary or interface conditions. The functional \(F(\vec{A})\) implements the boundary or interface conditions that fall into this category. This can be verified by converting the functional back to the boundary value problem.
It is impossible to apply no boundary condition to a particular boundary of the problem domain. Setting no boundary condition on a particular boundary is equivalent to application of the homogeneous Neumann boundary condition,
\[ \frac{1}{\mu}\hat{n}\times\bigg(\vec{\nabla}\times\vec{A}\bigg) = 0 \]
This natural boundary condition is implicitly implied by the first term of the functional above. The first therm of the functional is present even in the most minimalistic problems.
The compatibility condition for the curl-curl equation, i.e., equation (i) above with \( \eta = 0 \), suggests that the volume free-curent density must be derived as a curl of a current vector potential (Section 5.5 Compatibility conditions), i.e.,
\begin{equation} \vec{J}_f = \vec{\nabla} \times \vec{T}. \end{equation}
In this case the functional above must be modified by replacing the integral related to the free-current density with two other integrals:
\begin{equation} \iiint_{\Omega}\vec{J}_f \cdot \vec{A} dV \rightarrow \iiint_{\Omega} \vec{T} \cdot \bigg( \vec{\nabla} \times\vec{A} \bigg) dV - \iint_{\Gamma} \vec{T} \cdot \bigg(\hat{n}\times\vec{A}\bigg) dS. \end{equation}
Normally, the free-current density, \(\vec{J}_f\), is given by the formulation of the problem. In order to satisfy the compatibility condition for the curl-curl equation and derive the free-current density as a curl of the current vector potential, \(\vec{T}\), one needs to convert \(\vec{J}_f\) into \(\vec{T}\). This can be done by solving the following boundary value problem:
\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} = \hat{n} \times \vec{G} & \text{on} & \Gamma_{Dn} & \text{(ii)},\\ \text{(n)}& \hat{n}\times\bigg(\vec{\nabla}\times\vec{T}\bigg) = \vec{Q} & \text{on} & \Gamma_{Rm} & \text{(iii)}. \end{array} \end{equation}
A comparison between this boundary value problem and the three-dimensional boundary value problem formulated in terms of \(\vec{A}\) shown above reveals that these two boundary value problems can be solved by the same computer code. The class template StaticVectorSolver::Solver1 contains such code.

The second option above leads to the following boundary value problem:

\begin{equation} \begin{array}{lrcll} \text{ } & \vec{\nabla} \overset{V}{\times}\bigg(\dfrac{1}{\mu} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg) + \eta^2 \vec{A} = \vec{J}_f & \text{in} & \Omega & \text{(i)},\\ \text{(e)}&\hat{n} \overset{S}{\times} \vec{A} = \hat{n} \overset{S}{\times} \vec{G} & \text{on} & \Gamma_{Dn} & \text{(ii)},\\ \text{(n)}& \dfrac{1}{\mu}\hat{n}\overset{V}{\times}\bigg(\vec{\nabla}\overset{S}{\times}\vec{A}\bigg) + \gamma \hat{n} \overset{V}{\times} \bigg( \hat{n}\overset{S}{\times}\vec{A} \bigg) = \vec{Q} & \text{on} & \Gamma_{Rm} & \text{(iii)}, \\ \text{(e)}&\hat{n}\overset{S}{\times}\vec{A}_{+} = \hat{n}\overset{S}{\times}\vec{A}_{-}&\text{on}&\Gamma_{Ik}&\text{(iv)},\\ \text{(n)}&\dfrac{1}{\mu_{+}}\hat{n}\overset{V}{\times} \bigg( \vec{\nabla}\overset{S}{\times} \vec{A}_{+} \bigg) - \dfrac{1}{\mu_{-}}\hat{n}\overset{V}{\times} \bigg( \vec{\nabla}\overset{S}{\times} \vec{A}_{-} \bigg) = \vec{K}_f & \text{on} & \Gamma_{Ik} & \text{(v)}. \end{array} \end{equation}

The corresponding functional reads

\begin{equation} \begin{aligned} &F(\vec{A}) = \iint_{\Omega}\frac{1}{\mu}\bigg|\vec{\nabla}\overset{S}{\times}\vec{A}\bigg|^2 dS + \sum_m \int_{\Gamma_{Rm}} \bigg(\gamma\bigg|\hat{n}\overset{S}{\times}\vec{A}\bigg|^2 + 2\vec{Q}\cdot\vec{A}\bigg) dl + \eta^2 \iint_{\Omega}\mid\vec{A}\mid^2 dS -2\iint_{\Omega}\vec{J}_f \cdot \vec{A} dS - \\ & -2\sum_k\int_{\Gamma_{Ik}} \vec{K}_f \cdot \vec{A} dl. \end{aligned} \end{equation}

Most of the remarks listed above in the context of the three-dimensional boundary value problem hold for the two-dimensional version of the boundary value problem as well. In two-dimensions the free-current density is derived from a current vector potential as

\begin{equation} \vec{J}_f = \vec{\nabla} \overset{V}{\times} T. \end{equation}

Note that in this case the current vector potential, \( T \), is an out-of-plane vector, i.e., a scalar. If the current vector potential is used, the integral that takes into account the free-current density in the functional must replaced:

\begin{equation} \iint_{\Omega}\vec{J}_f \cdot \vec{A} dS \rightarrow \iint_{\Omega} T \bigg( \vec{\nabla}\overset{S}{\times}\vec{A} \bigg) dS - \int_{\Gamma} T \bigg(\hat{n}\overset{S}{\times} \vec{A}\bigg) dl. \end{equation}

The current vector potential \(T\) can be calculated by solving the following static scalar boundary value problem:

\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big(\vec{\nabla} T \big)= \vec{\nabla}\overset{S}{\times} \vec{J}_f & \text{in} & \Omega & \text{(i)},\\ \text{(e)} & T = G & \text{on} & \Gamma_{Dn} & \text{(ii)},\\ \text{(n)}&\hat{n}\cdot\vec{\nabla} T = Q & \text{on} &\Gamma_{Rm}&\text{(iii)}.\\ \end{array} \end{equation}

The class template StaticScalarSolver::Solver can help to solve this boundary value problem.

The boundary value problems based on the curl-curl partial differential equations discussed above can be solved with a help of the class template StaticVectorSolver::Solver1 or the class template StaticVectorSolver::Solver2.