Calculating the 3D Vector field of a permanent magnet with FEM

Question

trying to calculate the magnetic field of a permanent magnet using the finite element method I struggle setting up the equation system properly. In the end I want to be able to calculate the magnetic fields for any kind of shapes but I started with the easy system described bellow:

Needs["NDSolve`FEM`"]
(*Generate Geometry*)
pts = {{0, 0, 0}, {0, 1, 0}, {1, 1, 0}, {1, 0, 0}, {0, 0, 1}, {0, 1,

    1}, {1, 1, 1}, {1, 0, 1}};
pts1 = ScalingTransform[{3, 3, 3}][pts];
pts2 = TranslationTransform[{1, 1, 1}][pts];
hex = {{2, 3, 4, 1}, {1, 4, 8, 5}, {4, 3, 7, 8}, {3, 2, 6, 7}, {2, 1, 
    5, 6}, {5, 8, 7, 6}};
bmesh = BoundaryMeshRegion[Join[pts1, pts2], Polygon[hex], 
   Polygon[hex + 8], MeshCellStyle -> Opacity[0.3]];
HighlightMesh[bmesh, Labeled[0, "Index"]]
bmesh = bmesh = ToBoundaryMesh[bmesh]

To generate the following output:

According to the documentation (Wolfram Documentation) I continued defining markers for the different parts of the region

(*Generate Markers*)
cuboidCoordinates = {{1.5, 1.5, 1.5}};
airCoordinates = {0.5, 0.5, 0.5};
markerColors = {Red, Blue};
markerCoordinates = {{cuboidCoordinates}, {airCoordinates}};
Show[bmesh["Wireframe"],
Graphics3D[
MapThread[{PointSize[0.02], #1, Point /@ #2} &, {markerColors, 

    markerCoordinates}]

]
 ]

Before I started the meshing process:

markerSpecification = {{airCoordinates, 1 }, {cuboidCoordinates[[1]], 
   2}}
mesh = ToElementMesh[bmesh, "RegionMarker" -> markerSpecification]
mesh["Wireframe"]

When I now start to set up the differential equation I always run into errors:

I use:

\[Mu]Air = 4 \[Pi]*10^-7;
u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]};
M = Piecewise[{{{1, 0, 0}, ElementMarker == 2}}, {0, 0, 0}];
pde = Inactive[Curl][
    Inactive[Dot][1/[Mu]Air, Inactive[Curl][u, {x, y, z}]] - M, {x, 
     y, z}] == {0, 0, 0};
bcs = DirichletCondition[ 
   u[x, y, z] == { 0, 0, 0}, {x, y, z} [NotElement] mesh];
usol = NDSolveValue[{bcs, pde}, u, {x, y, z} [Element]  mesh]

Where M should be the magnetization of the permanent magnet and get the following error message:

I compared my solution to the Wolfram documentation as well as the answers in the following Question 3D FEM Vector Potential and unfortunately I can not see my mistake. So my questions are:

1. How can I make this code work?
2. Is my approach suitable for this type of problems?

**UPDATE: For a more advances approach to the computation of magnetic field, please have a look at the question follow up question over here **

Thanks in advance.

Best regards Tschibi2000

Answer 1

We can solve this problem using Coulomb gauge $\nabla.\vec{A}=0$ and continue approximation of magnetization instead of piecewise function from my answer here as follows

Needs["NDSolve`FEM`"]
 mesh = ToElementMesh[Cuboid[{0, 0, 0}, {3, 3, 3}], 
  MaxCellMeasure -> .001]
mesh["Wireframe"] 
u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]}; appro = 
 With[{k = 2. 10^4}, ArcTan[k #]/Pi + 1/2 &];
mx = Simplify`PWToUnitStep@
   PiecewiseExpand@
    If[1 <= x <= 2 && 1 <= y <= 2 && 1 <= z <= 2, 1, 0] /. 
  UnitStep -> appro; bmx[x_, y_, z_] := Curl[{mx, 0, 0}, {x, y, z}]
pde = Inactivate[Laplacian[u, {x, y, z}], Laplacian];
bcs = DirichletCondition[{ux[x, y, z] == 0, uy[x, y, z] == 0, 
    uz[x, y, z] == 0}, True];
{Ax, Ay, Az} = 
 NDSolveValue[{bcs, 
   Table[Activate[pde][[i]] == -bmx[x, y, z][[i]], {i, 3}]}, {ux, uy, 
   uz}, {x, y, z} [Element] mesh]
B = Evaluate[
   Curl[{Ax[x, y, z], Ay[x, y, z], Az[x, y, z]}, {x, y, z}]];

Visualization of vector potential and magnetic field

{VectorPlot3D[{Ax[x, y, z], Ay[x, y, z], 
   Az[x, y, z]}, {x, y, z} \[Element] mesh, 
  VectorStyle -> Arrowheads[0.01], VectorPoints -> Fine], 
 StreamPlot[{Ay[1.5, y, z], Az[1.5, y, z]}, {y, 0, 3}, {z, 0, 3}]}

{VectorPlot3D[B, {x, y, z} \[Element] mesh, 
  VectorStyle -> Arrowheads[0.01], VectorPoints -> 10], 
 StreamPlot[{B[[1]], B[[3]]} /. y -> 1.5, {x, .5, 2.5}, {z, .50, 2.5},
   VectorPoints -> Fine]}

With this method we can also solve more complex geometry including cylinder magnet and 2 iron bars:

Needs["NDSolve`FEM`"]
bmesh = Region[
  RegionUnion[Cylinder[{{0, 0, -0.25}, {0, 0, .25}}, 1.0], 
   Cuboid[{3.0, -1.5, -0.25}, {3.5, 1.5, .25}], 
   Cuboid[{-3.5, -1.5, 0}, {-3.0, 1.5, .5}]]]; Region[bmesh]
mesh = ToElementMesh[Cuboid[{-6, -3, -2}, {6, 3, 2}], 
  MaxCellMeasure -> .005]
m = mesh["Wireframe"]
Show[m, bmesh] 

Properties of magnet and iron bar

u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]}; appro = 
 With[{k = 2. 10^4}, ArcTan[k #]/Pi + 1/2 &];
mz = Simplify`PWToUnitStep@
   PiecewiseExpand@If[x^2 + y^2 <= 1 && -.25 <= z <= .25, 1, 0] /. 
  UnitStep -> appro; nu1 = 
 Simplify`PWToUnitStep@
   PiecewiseExpand@
    If[-3.5 <= x <= -3. && -1.5 <= y <= 1.5 && -.25 <= z <= .25, 
     1/1000, 1] /. UnitStep -> appro; nu2 = 
 Simplify`PWToUnitStep@
   PiecewiseExpand@
    If[3. <= x <= 3.5 && -1.5 <= y <= 1.5 && -.25 <= z <= .25, 1/1000,
      1] /. UnitStep -> appro; 
bmz[x_, y_, z_] := Curl[{0, 0, mz}, {x, y, z}]

PDEs, boundary condition and solution

pde = Inactivate[(nu1 + nu2) Laplacian[u, {x, y, z}] - 
   Cross[Grad[nu1 + nu2, {x, y, z}], Curl[u, {x, y, z}]], 
  Laplacian | Grad | Curl]; bcs = 
 DirichletCondition[{ux[x, y, z] == 0, uy[x, y, z] == 0, 
   uz[x, y, z] == 0}, True];
{Ax, Ay, Az} = 
 NDSolveValue[{bcs, 
   Table[Activate[pde][[i]] == -bmz[x, y, z][[i]], {i, 3}]}, {ux, uy, 
   uz}, {x, y, z} [Element] mesh]

Visualization of magnetic field

B = Evaluate[Curl[{Ax[x, y, z], Ay[x, y, z], Az[x, y, z]}, {x, y, z}]];
{VectorPlot3D[B, {x, y, z} [Element] mesh, 
  VectorStyle -> Arrowheads[0.01], VectorPoints -> 10], 
 StreamPlot[{B[[2]], B[[3]]} /. x -> 0, {y, -2.5, 2.5}, {z, -1.5, 
   1.5}, VectorPoints -> Fine]}