Return to Revisions

I am physicist from the first education, so it is apparently my field. As it follows from my experience in 3D FEM testing with application to the magnetic field calculation there is a problem with equation $\nabla \times (\nu \nabla \times \vec {A})=\vec {j}$. Therefore we prefer another form of this equation, for example, this one $\nabla \nu \times (\nabla \times \vec {A})+\nu \nabla \times \nabla \times\vec {A} =\vec {j}$ (ungauged form). Then if we have Coulomb gauge $\nabla.\vec {A}$, it automatically turns to $\nabla \nu \times (\nabla \times \vec {A})-\nu \nabla ^2\vec {A} =\vec {j}$ (Coulomb gauge). Now we can compare two form using mesh from Tim Laska answer (thanks to him), and function appro from xzczd answer (thanks to him as well). Let check Coulomb gauge first:

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

Visualization

Now let check ungauged form

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

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