Most often the problems in electrostatics are formulated in terms of the scalar electric potential, \(\Phi\). The problems in magnetostatics can be roughly sorted in five categories:

The problems that can be formulated in terms of the total scalar magnetic potential, \(\Psi\).
The problems that can be formulated in terms of the reduced scalar magnetic potential, \(\Theta\).
Two- dimensional planar problems that can be formulated in terms of the magnitude of the vector potential, \(A\).
Two- dimensional axisymmetric problems that can be formulated in terms of the scaled magnitude of the vector potential, \(A' = r A\), with \(r\) being the distance to the axis of rotation symmetry.
Three- dimensional problems and two-dimensional planar problems that can be formulated in terms of a vector magnetic potential, \(\vec{A}\).

We can add to this list the current potential, \(T\), which is used for formulating magnetostatic problems in terms of the magnetic vector potential, \(\vec{A}\), in planar two-dimensional problem domains. The problems in electrostatics formulated in term of \(\Phi\) and those problems in magnetostatics that can be formulated in therms of \(\Psi\), \(\Theta\), \(A\), \(A'\), and \(T\) can be described by a single general boundary value problem. It is shown below. This general boundary value problem can be solved with a help of the class template StaticScalarSolver::Solver.

All dielectric and magnetic materials are assumed to be linear and isotropic. Therefore, the permittivity, \(\epsilon\), and permeability, \(\mu\), of all materials are assumed to be real-valued smooth functions of spatial coordinates that are discontinuous on the interfaces between dissimilar materials.

The static scalar boundary value problem (Section 2.1 Electric scalar potential) reads

\begin{equation} \begin{array}{lrcll} \text{ }&- \vec{\nabla} \cdot \big( \epsilon \vec{\nabla} \Phi \big)= \rho_f & \text{in} & \Omega & \text{(i)},\\ \text{(e)} &\Phi = \eta & \text{on} & \Gamma_{Dn} & \text{(ii)},\\ \text{(n)}&\epsilon\hat{n}\cdot\vec{\nabla}\Phi+\gamma\Phi = \sigma&\text{on}&\Gamma_{Rm}&\text{(iii)},\\ \text{(e)} &\Phi_{+} = \Phi_{-} & \text{on} & \Gamma_{Ik} & \text{(iv)},\\ \text{(n)}&\epsilon_{+}\hat{n}\cdot\vec{\nabla}\Phi_{+}-\epsilon_{-}\hat{n}\cdot\vec{\nabla}\Phi_{-}=-\kappa_{f}&\text{on}&\Gamma_{Ik}&\text{(v)}. \end{array} \end{equation}

The other five scalar boundary value problems, i.e., problems formulated in terms of \(\Psi\), \(\Theta\), \(A\), \(A'\), and \(T\), are derived by substitutions listed in the table below.

\(\Phi\)	\(\epsilon\)	\(\rho_f\)	\(\eta\)	\(\sigma\)	\(\gamma\)	\(\kappa_f\)
\(\Psi\)	\(\mu\)	\(\rho_t\)	\(\eta_t\)	\(\sigma_t\)	\(\gamma\)	\(\kappa_t\)	Section 3.1 Total magnetic scalar potential
\(\Theta\)	\(\mu\)	\(\rho_r\)	\(\eta_r\)	\(\sigma_r\)	\(\gamma\)	\(\kappa_r\)	Section 3.2 Reduced magnetic scalar potential
\(A\)	\(\dfrac{1}{\mu}\)	\(J_f\)	\( G \)	\( Q \)	\(\gamma\)	\( K_f\)	Section 3.4 Magnetic vector potential in two dimensions
\(A'\)	\(\dfrac{1}{\mu'}\)	\(J_f\)	\( G \)	\( Q \)	\(\gamma\)	\( K_f\)	Section 3.4 Magnetic vector potential in two dimensions
\(T\)	\( 1\)	\(\vec{\nabla}\overset{S}{\times} \vec{J}_f \)	\( G \)	\( Q \)	0	(iv) and (v) are discarded	Section 5.5 Compatibility conditions

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 (ii) or the Robin (iii) boundary condition on one of the boundaries ensures the uniqueness of the solution.
The Neumann boundary condition, i.e., (iii) with \(\gamma=0\), ensures the uniqueness of the solution with a precision to a constant. However, if a Dirichlet or Robin boundary condition is specified on another boundary, the solution will be unique.
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 (surface in 2D) enclosed by the boundary. The same is true for vectors \(\hat{n}\) normal to interfaces.
The volume free-charge density, \(\rho_f\), does not include the surface free-charge density, \(\kappa_f\). Thinking that \(\rho_f = \rho_f'+\kappa_f\) is wrong.
The solution to this boundary value problem minimizes the following functional (Chapter 4 Variational formulations):
\begin{equation} F(\Phi) = \iiint_{\Omega}\epsilon\mid\vec{\nabla}\Phi\mid^2 dV +\sum_m \iint_{\Gamma_{Rm}}\big(\gamma\Phi^2 - 2 \sigma \Phi \big) dS - 2\sum_k \iint_{\Gamma_{Ik}}\kappa_f \Phi dS - 2\iiint_{\Omega}\rho_f \Phi dV. \end{equation}
The label (e) marks essential boundary or interface conditions. The functional \(F(\Phi)\) 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_Q finite elements of deal.II, for instance, do not allow any discontinuity. 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(\Phi)\) 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,
\[ \epsilon \hat{n}\cdot\vec{\nabla} \Phi = 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 boundary value problem above does not assume any coordinate system. That is, it is formulated in geometric terms. Therefore, the boundary value problem above is correct for thee- dimensional problems as well as for two- dimensional problems. The same cannot be said about the functional \(F(\Phi)\). It is valid only for three- dimensional problems. In general, we can reduce a three- dimensional problem to a two- dimensional problem if the initial three- dimensional problem exhibits a translation or rotation symmetry. A translation symmetry yields a two- dimensional planar domain. To write the functional for the two- dimensional planar domain we just need to remove one dimension of integration (Section 5.6 Numerical recipes):

\begin{equation} F(\Phi) = \iint_{\Omega}\epsilon\mid\vec{\nabla}\Phi\mid^2 dS +\sum_m \int_{\Gamma_{Rm}}\big(\gamma\Phi^2 - 2 \sigma \Phi \big) dl - 2\sum_k \int_{\Gamma_{Ik}} \kappa_f \Phi dl - 2\iint_{\Omega} \rho_f \Phi dS. \end{equation}

A rotation symmetry yields a two- dimensional axisymmetric problem domain. In the case of the axisymmetric problem domain we need to remove one dimension of integration and multiply all integrands by the distance to the axis of rotation symmetry, \(r\), (Section 5.6 Numerical recipes):

\begin{equation} F(\Phi) = \iint_{\Omega}\epsilon\mid\vec{\nabla}\Phi\mid^2 r dS +\sum_m \int_{\Gamma_{Rm}}\big(\gamma\Phi^2 - 2 \sigma \Phi \big)r dl - 2\sum_k \int_{\Gamma_{Ik}}\kappa_f \Phi r dl - 2\iint_{\Omega}\rho_f \Phi r dS. \end{equation}

As mentioned above the same boundary value problem can be used to describe the current vector potential, \(T\), in planar two-dimensional problems, i.e.,

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

The functional in this case reads

\begin{equation} F(T) = \iint_{\Omega}\mid \vec{\nabla} T \mid^2 dS -\sum_m \int_{\Gamma_{Rm}}\big( 2 Q T \big) dl- -2\Bigg[ \iint_{\Omega}\vec{J}_f\cdot\big( \vec{\nabla}\overset{V}{\times} T \big) dS - \oint_{\Gamma}\vec{J}_f \cdot \big( \hat{n} \overset{V}{\times} T \big) dl \Bigg]. \end{equation}

The class template StaticScalarSolver::Solver minimizes one of the four functionals above. Refer to the documentation of this class template for more detail.