Return to Answer

It looks reasonable but pay attention how we play with k and with Tanh[](Coulomb gauge) and ArcTan[] (ungauged form).

It looks reasonable but pay attention how we play with k and with Tanh[](Coulomb gauge) and ArcTan[] (ungauged form). For reference we can compare 3 solution for problem of coil magnetic field first considered by 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. This solution explained without code on https://physics.stackexchange.com/questions/513834/current-density-in-a-3d-loop-discretising-a-model/515657#515657 We have to calculate the vector potential and magnetic field of a rectangular coil with a current of 20A. The number of turns = 861. Inner cross section is 10.42cm×10.42cm, outer cross section is 15.24cm×15.24cm, coil height is 8.89cm. Here we show code for Closed Form Solution Algorithm (CFSA), BEM (Integral) and Mathematica FEM. CFSA code:

h = 0.0889; L1 = 0.1042; L2 = 0.1524; n = 861 (*16AWG wire*); J0 = \
20(*Amper*); j0 = 20*n/(h*(L2 - L1)/2); mu0 = 4 Pi 10^-7; b0 = j0 mu0;
bx[a_, b_, x_, y_, z_] := 
 z/(Sqrt[(-a + x)^2 + (-b + y)^2 + 
    z^2] (-b + y + Sqrt[(-a + x)^2 + (-b + y)^2 + z^2])) - z/(
  Sqrt[(a + x)^2 + (-b + y)^2 + 
    z^2] (-b + y + Sqrt[(a + x)^2 + (-b + y)^2 + z^2])) - z/(
  Sqrt[(-a + x)^2 + (b + y)^2 + 
    z^2] (b + y + Sqrt[(-a + x)^2 + (b + y)^2 + z^2])) + z/(
  Sqrt[(a + x)^2 + (b + y)^2 + 
    z^2] (b + y + Sqrt[(a + x)^2 + (b + y)^2 + z^2]))
by[a_, b_, x_, y_, z_] := 
 z/(Sqrt[(-a + x)^2 + (-b + y)^2 + 
    z^2] (-a + x + Sqrt[(-a + x)^2 + (-b + y)^2 + z^2])) - z/(
  Sqrt[(a + x)^2 + (-b + y)^2 + 
    z^2] (a + x + Sqrt[(a + x)^2 + (-b + y)^2 + z^2])) - z/(
  Sqrt[(-a + x)^2 + (b + y)^2 + 
    z^2] (-a + x + Sqrt[(-a + x)^2 + (b + y)^2 + z^2])) + z/(
  Sqrt[(a + x)^2 + (b + y)^2 + 
    z^2] (a + x + Sqrt[(a + x)^2 + (b + y)^2 + z^2]))
bz[a_, b_, x_, y_, 
  z_] := -((-b + y)/(
   Sqrt[(-a + x)^2 + (-b + y)^2 + 
     z^2] (-a + x + Sqrt[(-a + x)^2 + (-b + y)^2 + z^2]))) - (-a + 
   x)/(Sqrt[(-a + x)^2 + (-b + y)^2 + 
    z^2] (-b + y + Sqrt[(-a + x)^2 + (-b + y)^2 + z^2])) + (-b + y)/(
  Sqrt[(a + x)^2 + (-b + y)^2 + 
    z^2] (a + x + Sqrt[(a + x)^2 + (-b + y)^2 + z^2])) + (a + x)/(
  Sqrt[(a + x)^2 + (-b + y)^2 + 
    z^2] (-b + y + Sqrt[(a + x)^2 + (-b + y)^2 + z^2])) + (b + y)/(
  Sqrt[(-a + x)^2 + (b + y)^2 + 
    z^2] (-a + x + Sqrt[(-a + x)^2 + (b + y)^2 + z^2])) + (-a + x)/(
  Sqrt[(-a + x)^2 + (b + y)^2 + 
    z^2] (b + y + Sqrt[(-a + x)^2 + (b + y)^2 + z^2])) - (b + y)/(
  Sqrt[(a + x)^2 + (b + y)^2 + 
    z^2] (a + x + Sqrt[(a + x)^2 + (b + y)^2 + z^2])) - (a + x)/(
  Sqrt[(a + x)^2 + (b + y)^2 + 
    z^2] (b + y + Sqrt[(a + x)^2 + (b + y)^2 + z^2]))
da = (L2 - L1)/15/2;
dh = h/26/2; a = b = L1/2;
Bz[x_, y_, z_] := 
 Sum[bz[a + da (i - 1), b + da (i - 1), x, y, z + dh j], {i, 1, 
    16}, {j, -26, 26, 1}] + 
  Sum[bz[a, b, x, y, z + dh j], {j, -6, 6, 
    1}];

Code for BEM (Integral)

reg = RegionDifference[
   ImplicitRegion[-L2/2 <= x <= L2/2 && -L2/2 <= y <= L2/2 && -h/2 <= 
      z <= h/2, {x, y, z}], 
   ImplicitRegion[-L1/2 <= x <= L1/2 && -L1/2 <= y <= L1/2 && -h/2 <= 
      z <= h/2, {x, y, z}]];

j[x_, y_, z_] := Boole[{x, y, z} \[Element] reg]
jx[x_, y_, z_] := If[-y <= x <= y || y <= -x <= -y, Sign[y], 0]

jy[x_, y_, z_] := -jx[y, x, z]


Bx1[X_?NumericQ, Y_?NumericQ, Z_?NumericQ] := 
 b0/(4 Pi) NIntegrate[
    j[x, y, z] jy[x, y, 
      z] (Z - z)/(Sqrt[(x - X)^2 + (y - Y)^2 + (z - Z)^2])^3, {x, y, 
      z} \[Element] reg] // Quiet
By1[X_?NumericQ, Y_?NumericQ, 
  Z_?NumericQ] := -b0/(4 Pi) NIntegrate[
    j[x, y, z] jx[x, y, 
      z] (Z - z)/(Sqrt[(x - X)^2 + (y - Y)^2 + (z - Z)^2])^3, {x, y, 
      z} \[Element] reg] // Quiet
Bz1[X_?NumericQ, Y_?NumericQ, Z_?NumericQ] := 
 b0/(4 Pi) NIntegrate[
    j[x, y, z] (jx[x, y, z] (Y - y) - 
        jy[x, y, 
          z] (X - x))/(Sqrt[(x - X)^2 + (y - Y)^2 + (z - Z)^2])^3, {x,
       y, z} \[Element] reg] // Quiet

Code for FEM

eq1 = {Laplacian[A1[x, y, z], {x, y, z}] == -j[x, y, z] jx[x, y, z], 
   Laplacian[A2[x, y, z], {x, y, z}] == -j[x, y, z] jy[x, y, z]};
{Ax1, Ay1} = 
  NDSolveValue[{eq1, 
    DirichletCondition[{A1[x, y, z] == 0, A2[x, y, z] == 0}, 
     True]}, {A1, A2}, {x, y, z} \[Element] 
    ImplicitRegion[-2 L2 <= x <= 2 L2 && -2 L2 <= y <= 
       2 L2 && -2 L2 <= z <= 2 L2, {x, y, z}]];
B = Evaluate[Curl[{Ax1[x, y, z], Ay1[x, y, z], 0}, {x, y, z}]];

Now we compute and visualize data

lst1 = Table[{z1, -b0 B[[3]] /. {x -> 0, y -> 0, 
      z -> z1}}, {z1, -.3, .3, .01}];
lst2 = Table[{z1, Bz[0, 0, z1] mu0 20/(4 Pi)}, {z1, -.3, .3, .01}];
lst3 = Table[{z1, -Bz1[0, 0, z1]}, {z1, -.3, .3, .01}];

{Region[reg], 
 Show[ListLinePlot[lst2, PlotStyle -> Orange, Frame -> True, 
   Axes -> False], 
  ListPlot[{lst1, lst2, lst3}, Frame -> True, 
   FrameLabel -> {"z", "\!\(\*SubscriptBox[\(B\), \(z\)]\)"}, 
   PlotLegends -> {"FEM", "CFSA", "Integral"}]]}

u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]}; appro = 
 With[{k = 2. 10^4}, Tanh[k #]/2 + 1/2 &];
\[Nu]1 = Simplify`PWToUnitStep@
    PiecewiseExpand@If[x <= a && y <= a && z <= a, 1/(\[Mu]r), 1] /. 
   UnitStep -> 
    appro;(*permeability depending on iron cube in mesh*)\
\[CapitalGamma]d = {DirichletCondition[ux[x, y, z] == 0, y == 0], 
   DirichletCondition[ux[x, y, z] == -fluxDensity*b/2, y == b], 
   DirichletCondition[uy[x, y, z] == 0, x == 0], 
   DirichletCondition[uy[x, y, z] == fluxDensity*b/2, x == b], 
   DirichletCondition[uz[x, y, z] == 0, 
    y == b || y == 0 || x == 0 || x == b || z == 0 || z == b]};
\[CapitalGamma]n = {0, 0, 0};
op7 = Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] - \[Nu]1*
   Laplacian[u, {x, y, z}];(*Coulomb gauged*){mvpAx, mvpAy, mvpAz} = 
 NDSolveValue[{op7 == {0, 0, 0}, \[CapitalGamma]d}, {ux, uy, 
   uz}, {x, y, z} \[Element] mesh]; 

u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]}; appro = 
 With[{k = 1. 10^4}, Tanh[k #]/2 + 1/2 &];
\[Nu]1 = Simplify`PWToUnitStep@
    PiecewiseExpand@If[x <= a && y <= a && z <= a, 1/(\[Mu]r), 1] /. 
   UnitStep -> 
    appro;(*permeability depending on iron cube in mesh*)\
\[CapitalGamma]d = {DirichletCondition[ux[x, y, z] == 0, y == 0], 
   DirichletCondition[ux[x, y, z] == -fluxDensity*b/2, y == b], 
   DirichletCondition[uy[x, y, z] == 0, x == 0], 
   DirichletCondition[uy[x, y, z] == fluxDensity*b/2, x == b], 
   DirichletCondition[uz[x, y, z] == 0, 
    y == b || y == 0 || x == 0 || x == b || z == 0 || z == b]};
\[CapitalGamma]n = {0, 0, 0};
op7 = Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] - \[Nu]1*
   Laplacian[u, {x, y, z}];(*Coulomb gauged*){mvpAx, mvpAy, mvpAz} = 
 NDSolveValue[{op7 == {0, 0, 0}, \[CapitalGamma]d}, {ux, uy, 
   uz}, {x, y, z} \[Element] mesh]; 

u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]}; appro = 
 With[{k = 2. 10^4}, Tanh[k #]/2 + 1/2 &];
\[Nu]1 = Simplify`PWToUnitStep@
    PiecewiseExpand@If[x <= a && y <= a && z <= a, 1/(\[Mu]r), 1] /. 
   UnitStep -> 
    appro;(*permeability depending on iron cube in mesh*)\
\[CapitalGamma]d = {DirichletCondition[ux[x, y, z] == 0, y == 0], 
   DirichletCondition[ux[x, y, z] == -fluxDensity*b/2, y == b], 
   DirichletCondition[uy[x, y, z] == 0, x == 0], 
   DirichletCondition[uy[x, y, z] == fluxDensity*b/2, x == b], 
   DirichletCondition[uz[x, y, z] == 0, 
    y == b || y == 0 || x == 0 || x == b || z == 0 || z == b]};;
\[CapitalGamma]n = {0, 0, 0};
op7 = Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] - \[Nu]1*
   Laplacian[u, {x, y, z}];(*Coulomb gauged*){mvpAx, mvpAy, mvpAz} = 
 NDSolveValue[{op7 == {0, 0, 0}, \[CapitalGamma]d}, {ux, uy, 
   uz}, {x, y, z} \[Element] mesh]; 

u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]}; appro = 
 With[{k = 2. 10^4}, ArcTan[k #]/Pi + 1/2 &];
\[Nu]1 = Simplify`PWToUnitStep@
    PiecewiseExpand@If[x <= a && y <= a && z <= a, 1/(\[Mu]r), 1] /. 
   UnitStep -> 
    appro;(*permeability depending on iron cube in mesh*)\
\[CapitalGamma]d = {DirichletCondition[ux[x, y, z] == 0, y == 0], 
   DirichletCondition[ux[x, y, z] == -fluxDensity*b/2, y == b], 
   DirichletCondition[uy[x, y, z] == 0, x == 0], 
   DirichletCondition[uy[x, y, z] == fluxDensity*b/2, x == b], 
   DirichletCondition[uz[x, y, z] == 0, 
    y == b || y == 0 || x == 0 || x == b || z == 0 || z == b]};;
\[CapitalGamma]n = {0, 0, 0};
op7 = Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] - \[Nu]1*
   Laplacian[u, {x, y, z}]; op8 = 
 Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] + \[Nu]1*
   Curl[Curl[u, {x, y, z}], {x, y, z}];(*Coulomb gauged*){mvpAx, 
  mvpAy, mvpAz} = 
 NDSolveValue[{op8 == {0, 0, 0}, \[CapitalGamma]d}, {ux, uy, 
   uz}, {x, y, z} \[Element] mesh]; 

u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]}; appro = 
 With[{k = 2. 10^4}, Tanh[k #]/2 + 1/2 &];
\[Nu]1 = Simplify`PWToUnitStep@
    PiecewiseExpand@If[x <= a && y <= a && z <= a, 1/(\[Mu]r), 1] /. 
   UnitStep -> 
    appro;(*permeability depending on iron cube in mesh*)\
\[CapitalGamma]d = {DirichletCondition[ux[x, y, z] == 0, y == 0], 
   DirichletCondition[ux[x, y, z] == -fluxDensity*b/2, y == b], 
   DirichletCondition[uy[x, y, z] == 0, x == 0], 
   DirichletCondition[uy[x, y, z] == fluxDensity*b/2, x == b], 
   DirichletCondition[uz[x, y, z] == 0, 
    y == b || y == 0 || x == 0 || x == b || z == 0 || z == b]};
\[CapitalGamma]n = {0, 0, 0};
op7 = Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] - \[Nu]1*
   Laplacian[u, {x, y, z}];(*Coulomb gauged*){mvpAx, mvpAy, mvpAz} = 
 NDSolveValue[{op7 == {0, 0, 0}, \[CapitalGamma]d}, {ux, uy, 
   uz}, {x, y, z} \[Element] mesh]; 

u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]}; appro = 
 With[{k = 2. 10^4}, ArcTan[k #]/Pi + 1/2 &];
\[Nu]1 = Simplify`PWToUnitStep@
    PiecewiseExpand@If[x <= a && y <= a && z <= a, 1/(\[Mu]r), 1] /. 
   UnitStep -> 
    appro;(*permeability depending on iron cube in mesh*)\
\[CapitalGamma]d = {DirichletCondition[ux[x, y, z] == 0, y == 0], 
   DirichletCondition[ux[x, y, z] == -fluxDensity*b/2, y == b], 
   DirichletCondition[uy[x, y, z] == 0, x == 0], 
   DirichletCondition[uy[x, y, z] == fluxDensity*b/2, x == b], 
   DirichletCondition[uz[x, y, z] == 0, 
    y == b || y == 0 || x == 0 || x == b || z == 0 || z == b]};
\[CapitalGamma]n = {0, 0, 0};
op7 = Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] - \[Nu]1*
   Laplacian[u, {x, y, z}]; op8 = 
 Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] + \[Nu]1*
   Curl[Curl[u, {x, y, z}], {x, y, z}];(*Coulomb gauged*){mvpAx, 
  mvpAy, mvpAz} = 
 NDSolveValue[{op8 == {0, 0, 0}, \[CapitalGamma]d}, {ux, uy, 
   uz}, {x, y, z} \[Element] mesh];