Boundary Conditions in Magnetostatics - Calculating surface current density

I'm working on a finite element simulation of the magnetostatic magnetic vector potential in 3D, with the following geometry: .

The internal structure is a hollow cored-rectangular current loop. The outer box is the spacial domain, set to be approximately 10x the characteristic length scale of the problem. As a result I expect magnetic flux to be wholly contained within this box. My equation to be solved is the vector Poisson equation in 3D: $$\nabla^2\vec A=-\mu\vec J$$

My question then is with regards to the boundary conditions upon the faces. For faces 11-16 (the outer box) I choose Dirichlet zero conditions which correspond to the flux-enclosed approximation mentioned.

Faces 1-10, i.e the faces of my coil I'm unsure of. Assuming A is non-divergent, I expect the magnetic vector potential to be continuous across boundary but I also expect the normal derivative of A to be discontinuous by the surface current density K (D.J. Griffiths (1999),"Introduction to Electrodynamics", Upper Saddle River, NJ, Prentice Hall). This corresponds to a tangential discontinuity in the magnetic flux due to some K. However I have noticed as I work through the literature that many papers seem to ignore this component and simply require A to be continuous e.g:

N. Demerdash, T. Nehl and F. Fouad, "Finite element formulation and analysis of three dimensional magnetic field problems," in IEEE Transactions on Magnetics, vol. 16, no. 5, pp. 1092-1094, September 1980. doi: 10.1109/TMAG.1980.1060817

N. A. Demerdash, F. A. Fouad, T. W. Nehl and O. A. Mohammed, "Three Dimensional Finite Element Vector Potential Formulation of Magnetic Fields in Electrical Apparatus," in IEEE Transactions on Power Apparatus and Systems, vol. PAS-100, no. 8, pp. 4104-4111, Aug. 1981. doi: 10.1109/TPAS.1981.317005

The part I am unsure of is how to find the surface current density (assuming this is important). I know for my rectangular-cross sectioned coil with winding height z, inner length L_i, outer length L_o with N turns carrying a current I, the average current density through a winding section is: $$j_{0} = \frac{IN}{z(L_{o} - L_{i})/2}$$.

I also know the surface current density is also likely to vary with my faces in accordance with the geometries.

In general, the boundary conditions for the vector potential can be derived from the defining equation \begin{align} \mathrm{curl}\mathbf{H} &= \mathbf{J} \tag{1} \label{1}\\ \mathrm{div}\mathbf{B} &= 0 \tag{2}\label{2} \\ \mathbf{B}&=\mu\mathbf{H} \tag{3} \label{3}\\ \mathbf{B}&=\mathrm{curl}\mathbf{A} \tag{4} \label{4}\\ \mathrm{div}\mathbf{A} &= 0 \tag{5}\label{5} \\ \end{align} Here $$\eqref{1},\eqref{2}$$ are the static Maxwell's equations, $$\eqref{3}$$ assumes that the material has no hysteresis; Equation $$\eqref{4}$$ solves $$\eqref{2}$$,and finally $$\eqref{5}$$ is the Coulomb gauge that fixes the unavoidable arbitrariness of the vector potential in $$\eqref{4}$$.
The boundary conditions between two materials indexed by $$1$$ and $$2$$ now can be derived by using the vector potential in equations $$\eqref{1}, \eqref{3}$$ that is: \begin{align} \mathrm{curl}(\mu^{-1}\mathbf{\mathrm{curl}\mathbf{A}}) &= \mathbf{J} \tag{6} \label{6}\\ \end{align} From $$\mathrm{div}\mathbf{A} = 0$$ and the vector potential having a finite $$\mathrm{curl}$$ everywhere, the vector potential must also be continuous everywhere, even at a discontinuous interface, too, that is at the interface
$$\mathbf{A}_1=\mathbf{A}_2 \tag{7}\label{7}$$
Next using $$\eqref{6}$$ in regions where the free current density is zero, for example, inside and on the surface of the ferrite, and just outside of it, you get that the tangential component of the $$H$$ field is continuous, that is $$H_t^1=H_t^2$$ or equivalently $$\mathbf{n}\times \mathbf{H}_1=\mathbf{n}\times \mathbf{H}_2$$ where $$\mathbf{n}$$ is the local normal at the discontinuity. When written in terms of the vector potential this gives $$\frac{1}{\mu_1}\mathbf{n}\times\mathbf{\mathrm{curl}\mathbf{A}}_1=\frac{1}{\mu_2}\mathbf{n}\times\mathbf{\mathrm{curl}\mathbf{A}}_2 \tag{8}\label{8}$$
And finally the gauge itself is continuous. So within the magnetic material we have $$\mathrm{\mathbf{div}}\mathrm{\mathbf{grad}}\mathbf{A}=0 \tag{9}\label{9}$$ from which it follows that at the interface $$\mathbf{n}\cdot\mathrm{grad}(\mathbf{n}\cdot\mathbf{A}_1) = \mathbf{n}\cdot\mathrm{grad}(\mathbf{n}\cdot\mathbf{A}_2) \tag{10}\label{10}$$
Summarizing: the boundary conditions for the vector potential at a material interface with discontinuous permeability are $$\eqref{7}, \eqref{8}, \eqref{10}$$