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

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["NDSolveFEM"]
(*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?

Best regards Tschibi2000

• I think your equation is not set up correct, you dot multiply a scalar with a vector; also you dirichlet condition is not correct: you can not set u[x,y,z]={0,0,0}; you probably mean something like Thread[Equal[{ux[..].uy[],..},{0,0,0}]] Nov 19 '21 at 15:02
• Thanks you for your input. I tried an alternative approach with bcs = {DirichletCondition[ ux[x, y, z] == 0, {x, y, z} [NotElement] mesh], DirichletCondition[ uy[x, y, z] == 0, {x, y, z} [NotElement] mesh], DirichletCondition[ uz[x, y, z] == 0, {x, y, z} [NotElement] mesh]}; but it didnt work either. Your right the blunder is within the PDE itself I will have a look at i further ... Nov 20 '21 at 8:26
• What is $\mu$ for magnet? Nov 21 '21 at 8:18

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["NDSolveFEM"]
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 = SimplifyPWToUnitStep@
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["NDSolveFEM"]
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 = SimplifyPWToUnitStep@
PiecewiseExpand@If[x^2 + y^2 <= 1 && -.25 <= z <= .25, 1, 0] /.
UnitStep -> appro; nu1 =
SimplifyPWToUnitStep@
PiecewiseExpand@
If[-3.5 <= x <= -3. && -1.5 <= y <= 1.5 && -.25 <= z <= .25,
1/1000, 1] /. UnitStep -> appro; nu2 =
SimplifyPWToUnitStep@
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]}


• Thank you for your help. That was a big step in the right direction :-) Just some questions regarding your approach: Would it be possible to use the "ElementMarker" mechanism to define the magnetization of the permanent magnet? Could other elements with varying permeability be also defined and added? Nov 22 '21 at 7:39
• My plan is to extend my program to that point that I can adapt to more complex geometries combining permanent magnets and metal parts. Nov 22 '21 at 7:49
• @Tschibi Yes, it is possible to use "ElementMarker", but your function M = Piecewise[{{{1, 0, 0}, ElementMarker == 2}}, {0, 0, 0}] is not working. This problem discussed on mathematica.stackexchange.com/questions/230282/… Nov 22 '21 at 14:42
• @Tschibi Could you show geometry of your magnet with all parts and permeability? Nov 23 '21 at 3:44
• Yes of course. What would be the best way to share? Via STL File? Nov 24 '21 at 6:11