# Calculating the 3D magnetic vector field of a permanent magnet, with shape given by STL File

following my earlier question Link, I am trying to calculate the magnetic field of a permanent magnet. The shape of the magnet was generated with FreeCAD and is provided to Mathematica by an .STL file. The file can be downloaded here.

To test if NDSolve or import has a general problem with the .STL File I followed the documentation for solid mechanics and derived a displacement under load.

Needs["NDSolveFEM"]

(*Import STL File*)
region =
Import[FileNameJoin[{NotebookDirectory[], "Magnet.stl"}], {"STL", "BoundaryMeshRegion"}];

(*Set variables and parameters*)
vars = {{u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z}};
pars = <|"Material" -> Entity["Element", "Titanium"]|>;

(*Set Boundary conditions*)
Subscript[\[CapitalGamma], force] = SolidBoundaryLoadValue[x == 10, vars,pars, <|"Force" -> {0, 0, Quantity[-1000, "Newtons"]}|>];
Subscript[\[CapitalGamma], wall] = SolidFixedCondition[x == -10, vars, pars];

(*Set PDE components*)
op = SolidMechanicsPDEComponent[vars, pars];

(*Solve for displacement*)
regionDisplacement =
NDSolveValue[{op == Subscript[\[CapitalGamma], force], Subscript[\[CapitalGamma], wall]}, {u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z} \[Element] region];

(*Plot results*)
VectorDisplacementPlot3D[regionDisplacement, {x, y, z} \[Element] region]


While this is working:

I struggle to adapt this solution for the vector field of a permanent magnet. Following the solution that is linked as well as the solution here I imported the STL as a region and set up the magnetization as a an If condition containing RegionMember.

(*Import Magnet*)
magnet = Import[FileNameJoin[{NotebookDirectory[], "Magnet.stl"}], {"STL","BoundaryMeshRegion"}];

(*Setting up mesh*)
mesh = ToElementMesh[Cuboid[{-20, -10, -10}, {20, 10, 10}], MaxCellMeasure -> 1]

(*Setting up variables*)
u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]};

(*Setting up magnetization via approximation*)
appro = With[{k = 2. 10^4}, ArcTan[k #]/Pi + 1/2 &];

mx = SimplifyPWToUnitStep@PiecewiseExpand[If[RegionMember[magnet, {x, y, z}], 1, 0], Reals] /. UnitStep -> appro;
bmx[x_, y_, z_] := Curl[{mx, 0, 0}, {x, y, z}]

(*Setting up PDE and boundary conditions*)
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];

(*Solve and Plot System*)
{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}]];
VectorPlot3D[{Ax[x, y, z], Ay[x, y, z], Az[x, y, z]}, {x, y, z} \[Element] mesh, VectorStyle -> Arrowheads[0.01], VectorPoints -> Fine]


This unfortunately, this gives me back an empty solution:

Anyhow, If I replace the STL region by a dummy Cuboid region Cuboid[{-10, -4, -4}, {10, 4, 4}] it works.

(*Import Magnet*)
magnet = Import[FileNameJoin[{NotebookDirectory[], "Magnet.stl"}], {"STL","BoundaryMeshRegion"}];

(*Setting up mesh*)
mesh = ToElementMesh[Cuboid[{-20, -10, -10}, {20, 10, 10}], MaxCellMeasure -> 1]

(*Setting up variables*)
u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]};

(*Setting up magnetization via approximation*)
appro = With[{k = 2. 10^4}, ArcTan[k #]/Pi + 1/2 &];

mx = SimplifyPWToUnitStep@PiecewiseExpand[If[RegionMember[Cuboid[{-10, -4, -4}, {10, 4, 4}], {x, y, z}], 1, 0], Reals] /. UnitStep -> appro;
bmx[x_, y_, z_] := Curl[{mx, 0, 0}, {x, y, z}]

(*Setting up PDE and boundary conditions*)
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];

(*Solve and Plot System*)
{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}]];
VectorPlot3D[{Ax[x, y, z], Ay[x, y, z], Az[x, y, z]}, {x, y, z} \[Element] mesh, VectorStyle -> Arrowheads[0.01], VectorPoints -> Fine


I also tried to use RegionMarkers as an alternative to the RegionMember but in the end i always get the empty solution. In my humble oppinion the problem lies within PiecewiseExpand. While for easy geometries as cuboid it is able to change RegionMember into a piecewise function it fails when given the complex STL geometry. Without the approximation I was not able to derive a magnetic field, but the way the approximation is implemented I doubt it can be adapted for complex geometries.

Could anyone provide help with this topic? How can I derive Magnetic fields for complex shapes and magnetization.

P.S.: I use Mathematica 13.1 on Linux as well as Windows. It fails in both cases. You can download an example Notebook-File over here