# 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

Update: I tried another approach by generating the mesh with OpencascadeLink. Unfortunately the problem persists.

Needs["NDSolveFEM"]
Needs["OpenCascadeLink"]

FileNameJoin[{NotebookDirectory[], "Magnet.stl"}]];

air = OpenCascadeShape[Cuboid[{-20, -10, -10}, {20, 10, 10}]];

(*Create Boundary Mesh*)

(*Visualize Surfaces*)
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = {Opacity[0.75], ColorData["BrightBands"][#]} & /@ temp;
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]
mesh = ToElementMesh[bmesh, "RegionMarker" -> {{{0, 0, 0}, 1}}]
mesh["MeshElementMarkerUnion"]
parts = Map[
mesh["Wireframe"[ElementMarker == #[[1]],
"MeshElement" -> "MeshElements",
"ElementMeshDirective" ->
Directive[EdgeForm[], FaceForm[#[[2]]]]]] &, {{0, Gray}, {1,
Pink}}]
Rasterize[Show[parts, PlotRange -> {All, {0, 0}, All}]]

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

(*Setting up magnetization via Element Marker*)

mx = If[ElementMarker == 1, 1, 0]
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]


Update 2: An interesting fact about the solution with the dummie region is, that I many implemented shapes will work fine, even if I combine them via RegionUnion[], but once I want to use a Polyhedron[] the solution fails. After combining a lattice of spheres successfully with

Needs["NDSolveFEM"]
pos = Table[{i, j, k}, {i, -1, 1, 2}, {j, -1, 1, 2}, {k, -1, 1, 2}] //
Flatten[#, 2] & // Union;
magnets = Ball /@ pos;
magRegion = RegionUnion@magnets;

Graphics3D[magnets]
mesh = ToElementMesh[Cuboid[{-5, -5, -5}, {5, 5, 5}],
MaxCellMeasure -> 0.1]

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[RegionMember[magRegion, {x, y, z}], 1, 0],
Reals] /. 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}]];

{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[0, y, z], Az[0, y, z]}, {y, -5, 5}, {z, -5, 5}]}

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


and successfully checking, that Tetrahedron[] could be used. I was intrigued by the idea by converting the STL data to individual Polyhedrons or Tetrahedrons. Anyway it failed as soon as a Polyhedron[] is used, or two Tetrahedron[] will be combined to a Polyhedron[] by RegionUnion[].

Update 4: Thanks to the great answer of User21. I was able to identify my problem with "ElementMarker"-mechanism in NDSolve[]. The demonstrated solution is also adaptable to include µ-Values. Anyway there is one fallback from physics viewpoint at the moment. In order to simulate the behavior of a isotropic permanent magnet we need to define the magnetization vector function.

In the solution provided over here the magnetization Vector is equal to {1,0,0} and we see in the solution, that the derived magnetic field will show the expected behavior. That's because the mx Term with the dummy cuboid will evaluate to a function of mx[x,y,z] and will therefor be handled properly by Curl[]

mx = SimplifyPWToUnitStep@
PiecewiseExpand[
If[RegionMember[Cuboid[{-10, -4, -4}, {10, 4, 4}], {x, y, z}], 1,
0], Reals] /. UnitStep -> appro


If we use ElementMarker we can not simply omit the influence of the Curl within the PDE and therefor have to find a way to implement it. My first approach was to see the ElementMarker function as a function of x,y,z derive the Curl by hand and insert it. My first (obviously not very sophisticated) approach

m1CurlTest[x_, y_, z_] := Evaluate[
{If[ElementMarker == 1, -my1 + mz1, 0],
If[ElementMarker == 1, -mz1 + mx1, 0],
If[ElementMarker == 1, -mx1 + my1, 0]
}
]


was solvable by NDSolve[] but obviously wrong, because a magnetization vector of {1,1,1} will lead to a magnetization of 0.

My overall plan is to extend this solution to a more generalized toolbox to model permanent magnets in combination with metal parts in static conditions. While I could do this with other, more engineering like programs I hope to be able to use the optimization abilities of Mathematica in my favor for some types of problems.

• OK, you created the the CAD in FreeCAD, right? FreeCAD also uses OpenCascade. Try and see if FreeCAD allows to export a *.brep file format. That is the native file format for OpenCascasde. We can then try to import that brep file in Mathematica's OpenCascade and do the region operation there. ElementMarker is definitely the way to go. Upload the brep file and I can help you. Sep 3, 2022 at 16:50
• Another question: Do you want the DirichletCondition be valid on all (also the internal) boundaries of the object. True will do that, but perhaps you only want the bc at the outer faces of the cube? See this Sep 3, 2022 at 17:14
• This 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; Curl[{mx, 0, 0}, {x, y, z}] give a different result than this: mx = SimplifyPWToUnitStep@ PiecewiseExpand[If[ElementMarker == 1, 1, 0], Reals] /. UnitStep -> appro; Curl[{mx, 0, 0}, {x, y, z}] Sep 3, 2022 at 17:24
• Thank you very much for your help. I uploaded the .brep file over here: wolframcloud.com/obj/b34298f8-3087-47fb-ac16-05151d272471 Unfortunately I get an error if I want to use this type of file instead of .stl file. ToElementMesh[] says the mesh quality is below 0. Sep 4, 2022 at 11:55
• About the DirichletCondition you are right it should be only applied to the outer faces of the qube. As the working first working solution did not contain internal boundaries, the True value was suitable, but I should replace it with something like: x == -20 && x == 20 && y == -10 && y == 10 && z == -10 && z == 10 Sep 4, 2022 at 12:08

With the help of user21 I was able to find a solution for my problem. Anyway there is still room for improvement as it is still a little hacky in my humble opinion.

The problem that user21 pointed out rightly was the Curl[] - term that takes care of the magnetization within the PDE. It evaluates (correctly) to 0 within my first approach, as the components mx, my and mz are no function of x, y and z (et least not from the perspective of Curl[]). Calculating the derivatives of a constant with respect to x,yz will therefor always lead to 0.

The solution from the former post, bypasses this problem by not using ElementMarker and defining the components of the magnetization in the following way (example for x component of magnetization):

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



Like This mx will derived to a function of the three coordinates aka mx[x,y,z], which will have a non 0 solution when Curl[] is applied.

So my main goal was to find a way, feeding the different magnetization values for the whole grid in a form of {mx[x,y,z],my[x,y,z],mz[x,y,z]} to the PDE. To do this I used EvaluateOnElementMesh[]:

(*Setting up geometry*)
Needs["NDSolveFEM"]
Needs["OpenCascadeLink"]

(*Magnet*)

(*Surrounding Air*)
cubeShape = OpenCascadeShape[Cuboid[{-20, -10, -10}, {20, 10, 10}]];

(*Stitch the shapes together in one mesh*)

(*Set up parameters*)

regions = {"Magnet", "Air"};
regionCoords = <|"Magnet" -> {0, 0, 0}, "Air" -> {-19, -9, -9}|>;
regionMarker = <|"Magnet" -> 1, "Air" -> 2|>;
regionsColour = <|"Magnet" -> Orange, "Air" -> Blue|>;
magnetization = {1, 0, 0};

(*Plot Mesh with markers*)
Show[Graphics3D[{regionsColour[[#]], PointSize[0.02], Point[regionCoords[[#]]]} & /@ regions], bmesh["Wireframe"], Boxed -> False]

(*Generate Mesh*)
mesh = ToElementMesh[bmesh, "RegionMarker" -> ({regionCoords[#], regionMarker[#]} & /@ regions)]

(*Generate Interpolated functions for Magnetization over mesh*)
{magnetInterpolX, magnetInterpolY, magnetInterpolZ} =
EvaluateOnElementMesh[{x, y, z}, If[ElementMarker == 1, #, 0] & /@ magnetization, mesh]

(*Setting up PDE and boundary conditions*)
u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]};
pde = Laplacian[u, {x, y, z}];
bcs = DirichletCondition[{ux[x, y, z] == 0, uy[x, y, z] == 0, uz[x, y, z] == 0}, x == -20 || x == 20 || y == -10 || y == 10 || z == -10 || z == 10];

{Ax, Ay, Az} = NDSolveValue[{bcs, pde == -Curl[{magnetInterpolX[x, y, z], magnetInterpolY[x, y, z], magnetInterpolZ[x, y, z]}, {x, y, z}]}, {ux, uy, uz}, {x, y, z} \[Element] mesh]

(*Plotting the Results*)

B[x_, y_, z_] =
Curl[{Ax[x, y, z], Ay[x, y, z], Az[x, y, z]}, {x, y, z}];
VectorPlot3D[B[x, y, z], {x, y, z} \[Element] mesh, VectorStyle -> Arrowheads[0.01], VectorPoints -> 10, BoxRatios -> {2, 1, 1}]
SliceVectorPlot3D[B[x, y, z], {x, y, z} \[Element] mesh, VectorStyle -> Arrowheads[0.01], VectorPoints -> 10, BoxRatios -> {2, 1, 1}]
VectorPlot[{B[x, y, 0][[1]], B[x, y, 0][[2]]}, {x, -20, 20}, {y, -10, 10}, VectorPoints -> Fine, AspectRatio -> Automatic]
VectorPlot[{B[x, 0, z][[1]], B[x, 0, z][[3]]}, {x, -20, 20}, {z, -10, 10}, VectorPoints -> Fine, AspectRatio -> Automatic]
VectorPlot[{B[0, y, z][[2]], B[0, y, z][[3]]}, {y, -10, 10}, {z, -10, 10}, VectorPoints -> Fine, AspectRatio -> Automatic]


Anyway this solution has some drawbacks.

• The generation of the interpolated functions magnetInterpolX, magnetInterpolY and magnetInterpolZ takes more time and resources.
• Looking at the ElementMarker-mechanism and other software on the market my solution is definitely a hack. There must be a more elegant way that makes use of ElementMarker for that type of problems instead of hacking around it like I did.

I will continue with implementing a term for the permeability and see if the results are compareable to 3rd Party software. I will post updates to this answer once I did this.

@Everyone who it may concern: If you have any idea how to improve this solution please feel free to take part here. Possible improvements could include:

• another handling of ElementMarker within the build in functions
• changing the PDE bypass the Curl[] term

UPDATE 1: Added code for faster generation of magnetInterpolX, magnetInterpolY and magnetInterpolZ according to suggestion of User21 in the comments

UPDATE 2: I Was able to continue with the task and did the following things: I implemented permeability and set up a second solution using the scalar magnetic potential.

Solution 1 (Vector Potential):

Needs["NDSolveFEM"]

(*Define Regions*)
diskRegion = Cylinder[{{0, 0, 3.097}, {0, 0, 3.097 + 6}}, 3];
airRegion = Cuboid[{-10, -10, -10}, {10, 10, 10}];
dummyRegion = Cuboid[{-5, -5, -2}, {5, 5, 2}];

(*Region and Mesh Properties*)
regionNames = {"Disk", "Air", "Dummy"};
regions = <|"Disk" -> diskRegion, "Air" -> airRegion, "Dummy" -> dummyRegion|>;

regionCoords = <|"Disk" -> RegionCentroid[regions["Disk"]], "Air" -> {9, -9, -9}, "Dummy" -> RegionCentroid[regions["Dummy"]]|>;

regionMarker = <|"Disk" -> 1, "Air" -> 2, "Dummy" -> 3|>;
regionsColour = <|"Disk" -> Red, "Air" -> Blue, "Dummy" -> Green|>;
regionsMagnetization = <|"Disk" -> {0, 0, -1.32}, "Air" -> {0, 0, 0},
"Dummy" -> {0, 0, 0}|>;
regionsMuR = <|"Disk" -> 1.1, "Air" -> 1, "Dummy" -> 1|>;
regionsResolution = <|"Disk" -> 0.1, "Air" -> 0.5, "Dummy" -> 0.05|>;
maxCellMeasure = 1;
maxBoundaryCellmeasure = 0.05;

mu0 = 1; (*Give Magnetization in Tesla Remanent Flux instead of \
Ampere/Meter*)

(*Definition of Mesh dependent µ*)
muDef[x_, y_, z_] :=
Piecewise[{regionsMuR[#]*mu0, ElementMarker == regionMarker[#]} & /@
regionNames, mu0]

(*Visualize System*)
Show[Graphics3D[{regionsColour[[#]], PointSize[0.02], Point[regionCoords[[#]]]} & /@ regionNames], bmesh["Wireframe"],
Boxed -> False]
mesh = ToElementMesh[bmesh, "RegionMarker" -> ({regionCoords[#], regionMarker[#], regionsResolution[#]} & /@ regionNames), MaxCellMeasure -> maxCellMeasure, "MaxBoundaryCellMeasure " ->maxBoundaryCellmeasure ]

(*Function for Remanent Flux of Magnet dependent on Mesh*)
{magnetInterpolX, magnetInterpolY, magnetInterpolZ} = EvaluateOnElementMesh[{x, y, z}, {(Piecewise[{regionsMagnetization[#][[1]],         ElementMarker == regionMarker[#]} & /@ regionNames, {0, 0, 0}]),
(Piecewise[{regionsMagnetization[#][[2]], ElementMarker == regionMarker[#]} & /@ regionNames, {0, 0, 0}]),
(Piecewise[{regionsMagnetization[#][[3]], ElementMarker == regionMarker[#]} & /@ regionNames, {0, 0, 0}])}, mesh]

(*setting up pde variables*)
u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]};
pde = -1/muDef[x, y, z] Laplacian[u, {x, y, z}] -
Cross[Grad[-1/(muDef[x, y, z]) , {x, y, z}], Curl[u, {x, y, z}]];
bcs = DirichletCondition[{ux[x, y, z] == 0, uy[x, y, z] == 0,
uz[x, y, z] == 0},
x == -20 || x == 20 || y == -10 || y == 10 || z == -10 || z == 10];

(*Solving PDE*)
{Ax, Ay, Az} = NDSolveValue[{bcs, pde == Curl[{magnetInterpolX[x, y, z], magnetInterpolY[x, y, z], magnetInterpolZ[x, y, z]}, {x, y, z}]}, {ux, uy, uz}, {x, y, z} \[Element] mesh]

(*Visualization of B Field*)
B[x_, y_, z_] =
Curl[{Ax[x, y, z], Ay[x, y, z], Az[x, y, z]}, {x, y, z}];

VectorPlot3D[B[x, y, z], {x, y, z} \[Element] mesh, VectorStyle -> Arrowheads[0.01], VectorPoints -> 10]
SliceVectorPlot3D[B[x, y, z], {x, y, z} \[Element] mesh, VectorStyle -> Arrowheads[0.01], VectorPoints -> 10]


Solution 2 (Magnetic Scalar Potential)

Needs["NDSolveFEM"]
Needs["OpenCascadeLink"]

(*Define Regions*)
diskRegion = Cylinder[{{0, 0, 3.097}, {0, 0, 3.097 + 6}}, 3];
airRegion = Cuboid[{-10, -10, -10}, {10, 10, 10}];
dummyRegion = Cuboid[{-5, -5, -2}, {5, 5, 2}];

(*Region and Mesh Properties*)
regionNames = {"Disk", "Air", "Dummy"};
regions = <|"Disk" -> diskRegion, "Air" -> airRegion, "Dummy" -> dummyRegion|>;

regionCoords = <|"Disk" -> RegionCentroid[regions["Disk"]], "Air" -> {9, -9, -9}, "Dummy" -> RegionCentroid[regions["Dummy"]]|>;

regionMarker = <|"Disk" -> 1, "Air" -> 2, "Dummy" -> 3|>;
regionsColour = <|"Disk" -> Red, "Air" -> Blue, "Dummy" -> Green|>;
regionsMagnetization = <|"Disk" -> {0, 0, -1.32}, "Air" -> {0, 0, 0},
"Dummy" -> {0, 0, 0}|>;
regionsMuR = <|"Disk" -> 1.1, "Air" -> 1, "Dummy" -> 1|>;
regionsResolution = <|"Disk" -> 0.1, "Air" -> 0.5, "Dummy" -> 0.05|>;
maxCellMeasure = 1;
maxBoundaryCellmeasure = 0.05;

mu0 = 1; (*Give Magnetization in Tesla Remanent Flux instead of \
Ampere/Meter*)

(*Definition of Mesh dependent µ*)
muDef[x_, y_, z_] :=
Piecewise[{regionsMuR[#]*mu0, ElementMarker == regionMarker[#]} & /@
regionNames, mu0]

(*Visualize System*)
Show[Graphics3D[{regionsColour[[#]], PointSize[0.02], Point[regionCoords[[#]]]} & /@ regionNames], bmesh["Wireframe"],
Boxed -> False]
mesh = ToElementMesh[bmesh, "RegionMarker" -> ({regionCoords[#], regionMarker[#], regionsResolution[#]} & /@ regionNames), MaxCellMeasure -> maxCellMeasure ,"MaxBoundaryCellMeasure " ->maxBoundaryCellmeasure ]

(*Function for Remanent Flux of Magnet dependent on Mesh*)
{magnetInterpolX, magnetInterpolY, magnetInterpolZ} = EvaluateOnElementMesh[{x, y, z}, {(Piecewise[{regionsMagnetization[#][[1]],         ElementMarker == regionMarker[#]} & /@ regionNames, {0, 0, 0}]),
(Piecewise[{regionsMagnetization[#][[2]], ElementMarker == regionMarker[#]} & /@ regionNames, {0, 0, 0}]),
(Piecewise[{regionsMagnetization[#][[3]], ElementMarker == regionMarker[#]} & /@ regionNames, {0, 0, 0}])}, mesh]

(*Setting up pde*)
pde = -Div[{{muDef[x, y, z], 0, 0}, {0, muDef[x, y, z], 0}, {0, 0, muDef[x, y, z]}} . (Grad[vm[x, y, z], {x, y, z}] + {magnetInterpolX[x, y, z], magnetInterpolY[x, y, z], magnetInterpolZ[x, y, z]}), {x, y, z}];
bcs = DirichletCondition[{vm[x, y, z] == 0}, x == -10 || x == 10 || y == -10 || y == 10 || z == -10 || z == 10];

(*Solving pde*)
vmOut = NDSolveValue[{bcs, pde == 0}, vm, {x, y, z} \[Element] mesh]

MuMesh[x_, y_, z_] = EvaluateOnElementMesh[{x, y, z}, muDef[x, y, z], mesh][x, y, z];
B[x_, y_, z_] = ((MuMesh[x, y, z])*Grad[vmOut[x, y, z], {x, y, z}] + MuMesh[x, y, z]*{magnetInterpolX[x, y, z], magnetInterpolY[x, y, z], MagnetInterpolZ[x, y, z]});

VectorPlot3D[B[x, y, z], {x, y, z} \[Element] mesh, VectorStyle -> Arrowheads[0.01], VectorPoints -> 10]
SliceVectorPlot3D[B[x, y, z], {x, y, z} \[Element] mesh, VectorStyle -> Arrowheads[0.01], VectorPoints -> 10]


The first difference was, that the solution with the scalar magnetic Potential was faster (51s vs 60 s). Anyway there are some problems with numeric results I was not able to counter until now. Please find below a comparison of the results with values given by CST Simulation Suite and an analytical solution provided by this Paper. You see the comparison of results along the z-Axis from -2cm to 2cm. The magnetic field is given in Tesla

1 Vector Potential:

2 Scalar Potential:

The results are in good agreement with the other sources. We see, that the numerical errors within the solution for scalar potential are higher. While we are very close to get comparable results there is still room for improvement. Does somebody know additional steps that are used to improve the results in terms of numeric stability and accuracy?

Update 3: By searching for the source of the numerical problems I encountered a Problem with the EvaluateOnElementMesh[] function. In the picture below we can see the z component of the magnetization as well as the calculated B-Field along the z Axis.

We can see clearly see that the calculated magnetization has not the wanted square shape. As the discontinuities correlate with the spikes on the B(z) component I guess they are most likely the reason for the numerical errors. Unfortunately the EvaluateOnElementMesh[] function is not documented yet and also Option[EvaluateOnElementMesh] gives back {}. If anyone has an idea how to fix this write me please. I try to circumvent the problem somehow.

• OK, I'll continue to look into this. For now: Here is a much more efficient way to generate the interpolating functions: {magnetInterpolX, magnetInterpolY, magnetInterpolZ} = EvaluateOnElementMesh[{x, y, z}, magnetization,mesh]. Please also see this Sep 7, 2022 at 11:31
• Yes. This function is faster. I adapted the code in my answer accordingly. Anyway I could not implement it in the way you did, because I got an error, that EvaluateOnElementMesh does not evaluate to an scalar. Sep 7, 2022 at 12:11
• Your answer brings the improvement, that we can get rid of the additional regions as well. Overall this feels much less hacky than before ... Sep 7, 2022 at 12:16
• I was able to work on this topic again. Please see my update to this answer. Nov 11, 2022 at 13:02
• I found a possible cause for the numeric errors. Please find update 3 in the answer above. Nov 15, 2022 at 14:04

@user21 and I have been working on new features to better model these problems. Such is a DiscontinuousInterpolatingFunction representing an approximate discontinuous function whose values are found by interpolation. It is very useful when modeling multi-material domains, such as in electromagnetics, when we have, for example, a magnetization or a relative permeability that is different for each material we have. According to theory, the B and H fields can have discontinuities at interfaces of different materials, either the normal or the tangential components. This new functionality addresses this behavior seen in the field at interfaces

Also, in this answer, I added a new way to generate the curl-curl equation of the magnetic vector potential. I must say that I don't recommend using the A-formulation in 3D because Mathematica does not have Edge/Curle/Vector elements, and these types of problems are better solved with these elements. Using the A-formulation in 3D with Mathematica elements (nodal elements) can lead to spurious results, and the specification of boundary conditions is not straightforward. So, if the problem is free of currents, it is better to use a scalar magnetic potential Vm; as this formulation is a Div-Grad equation, nodal elements can be used, and boundary conditions are naturally handled.

This answer does not use any approximation for specifying the magnetization vector, which I found great because the definition is simple, and it gives a good result.

References for the new discontinuous interpolating function can be found in the new documentation at FEMDocumentation/ref/DiscontinuousInterpolatingFunction or FEMDocumentation/ref/EvaluateOnElementMesh. This answer uses element markers because they are necessary for creating the discontinuous function previously mentioned.

Create the boundary mesh:

Needs["NDSolveFEM"]
cubeShape = OpenCascadeShape[Cuboid[{-20, -10, -10}, {20, 10, 10}]];


Define region markers and visualize:

regions = {"Magnet", "Air"};
regionCoords = <|"Magnet" -> {0, 0, 0}, "Air" -> {-19, -9, -9}|>;
regionMarker = <|"Magnet" -> 1, "Air" -> 2|>;
regionsColour = <|"Magnet" -> Orange, "Air" -> Blue|>;
Show[Graphics3D[{regionsColour[[#]], PointSize[0.02],
Point[regionCoords[[#]]]} & /@ regions], bmesh["Wireframe"],
Boxed -> False]


Create the mesh:

mesh = ToElementMesh[bmesh,
"RegionMarker" -> ({regionCoords[#], regionMarker[#]} & /@
regions)];


Define variables:

vars = {{Ax[x, y, z], Ay[x, y, z], Az[x, y,
z]}, {x, y, z}};


Define the vacuum permeability and the magnetization vector $$\vec{M}$$:

\[Mu]0 = 4.0*\[Pi]*10^-7;
mx = If[ElementMarker == 1, 1, 0];
bmx[x_, y_, z_] := Evaluate[{mx, 0, 0}]


Define the coefficients to generate the first term of the equation $$\nabla\times\left(\mu_0^{-1}\nabla \times \vec{A}\right)$$:

coefficients = {{{{0, 0, 0}, {0, 1/\[Mu]0, 0}, {0, 0, 1/\[Mu]0}}, {{0,
0, 0}, {-1/\[Mu]0, 0, 0}, {0, 0, 0}}, {{0, 0, 0}, {0, 0,
0}, {-1/\[Mu]0, 0, 0}}}, {{{0, -1/\[Mu]0, 0}, {0, 0, 0}, {0, 0,
0}}, {{1/\[Mu]0, 0, 0}, {0, 0, 0}, {0, 0, 1/\[Mu]0}}, {{0, 0,
0}, {0, 0, 0}, {0, -1/\[Mu]0, 0}}}, {{{0, 0, -1/\[Mu]0}, {0, 0,
0}, {0, 0, 0}}, {{0, 0, 0}, {0, 0, -1/\[Mu]0}, {0, 0,
0}}, {{1/\[Mu]0, 0, 0}, {0, 1/\[Mu]0, 0}, {0, 0, 0}}}};


The second term of the equation is$$-\nabla \times \left(\vec{M}\right)$$.

Generate the equation:

op = DiffusionPDETerm[vars,
coefficients] - {Inactive[
Div][{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}} . bmx[x, y, z], {x, y,
z}], Inactive[Div][{{0, 0, -1}, {0, 0, 0}, {1, 0, 0}} .
bmx[x, y, z], {x, y, z}],
Inactive[Div][{{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}} .
bmx[x, y, z], {x, y, z}]};


Define boundary conditions:

bcs = DirichletCondition[{Ax[x, y, z] == 0, Ay[x, y, z] == 0,
Az[x, y, z] == 0},
x == -20 || x == 20 || y == -10 || y == 10 || z == -10 || z == 10]


Solve the PDE:

{AxFun, AyFun, AzFun} =
NDSolveValue[{op == {0, 0, 0}, bcs}, {Ax, Ay,
Az}, {x, y, z} \[Element] mesh];


Generate the Discontinuous interpolating functions for each term with DisconitnuousInterpolatingFunction:

AxDIF = DiscontinuousInterpolatingFunction[AxFun];
AyDIF = DiscontinuousInterpolatingFunction[AyFun];
AzDIF = DiscontinuousInterpolatingFunction[AzFun];


In this case, the DIFs will prioritize the magnet region, so every time one evaluates a value at the interface, the value with material parameters of the magnet will be considered. This can be changed. Now that we have defined the DIFs, we can proceed to take the derivative to compute the B field.

Next, compute the discontinuous B field:

BField =
Curl[{AxDIF[x, y, z], AyDIF[x, y, z], AzDIF[x, y, z]}, {x, y, z}];


Visualize the field:

SliceVectorPlot3D[BField, {x, y, z} \[Element] mesh,
VectorStyle -> Arrowheads[0.01], VectorPoints -> 10,
BoxRatios -> {2, 1, 1}]


Other visualizations:

DensityPlot[
Evaluate[
Sqrt[BField[[1]]^2 + BField[[2]]^2 + BField[[3]]^2] /. {z ->
0}], {x, -20, 20}, {y, -10, 10}, PlotRange -> All,
AspectRatio -> Automatic, PlotLegends -> Automatic,
ColorFunction -> "Rainbow", PlotPoints -> 90]


VectorPlot[
Evaluate[{BField[[1]], BField[[2]]} /. {z -> 0}], {x, -20,
20}, {y, -10, 10}, VectorPoints -> Fine, AspectRatio -> Automatic]


VectorPlot[
Evaluate[{BField[[1]], BField[[3]]} /. {y -> 0}], {x, -20,
20}, {z, -10, 10}, VectorPoints -> Fine, AspectRatio -> Automatic]


These plots show the natural behavior of a permanent magnet. In the plots graphic, we can see that the B-Field lines are the lines of a magnet with magnetization of {1,0,0}.

• Thank you very much for this great post. Do I get it right, that this alternative formulation of the problem will not solve the problem with the discontinuities? Jan 17 at 15:27
• It solves the problem with the discontinuities to some extent. One can address discontinuities in the dependent variables and its derivatives. Also, if one would like to have a function for a material parameter that has different values for different regions one can use EvaluateOnElementMesh, and this also will reproduce a discontinuous interpolating function. You can use it for the magnetization M, and the jumps you saw in your plot of Mz will disappear. Jan 18 at 19:49
• There are other discontinuities in FEM that are not related to different materials. @Tschibi, you saw a few of them when plotting the magnetic field of your cylinder and comparing them with analytical solutions. These little jumps you see are discontinuities that are a consequence of taking derivatives of your primary solution. When you take derivatives you lose accuracy. If you had a 2-order interpolation function, the derivative of this function now will be 1st order. Jan 18 at 19:53
• @Tschibi, you can see a more detailed explanation of the Discontinuity of secondary unknown quantities at reference.wolfram.com/language/PDEModels/tutorial/… Jan 18 at 20:03
• Thank you very much for your explanation. In some of my use cases increasing the number of mesh cells is not suitable. Other simulation software do give better results with much less mesh cells / calculation time then. Is there a possibility that somewhere in the future higher orders of interpolation can be used to counter this type of problem in another way than increasing the mesh cells? Jan 19 at 20:53

The main problem is that bmx does not evaluate to what you think.

Create the boundary mesh:

Needs["NDSolveFEM"]
Needs["OpenCascadeLink"]
cubeShape = OpenCascadeShape[Cuboid[{-20, -10, -10}, {20, 10, 10}]];


Visualize

regions = {"Magnet", "Air"};
regionCoords = <|"Magnet" -> {0, 0, 0}, "Air" -> {-19, -9, -9}|>;
regionMarker = <|"Magnet" -> 1, "Air" -> 2|>;
regionsColour = <|"Magnet" -> Orange, "Air" -> Blue|>;
Show[
Graphics3D[{regionsColour[[#]], PointSize[0.02],
Point[regionCoords[[#]]]} & /@ regions],
bmesh["Wireframe"]
, Boxed -> False]


The mesh:

mesh = ToElementMesh[bmesh,
"RegionMarker" -> ({regionCoords[#], regionMarker[#]} & /@
regions)];


The PDE. Note the BC has not True as a predicate but the outer faces:

(*Setting up PDE and boundary conditions*)
u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]};
pde = Laplacian[u, {x, y, z}];
bcs = DirichletCondition[{ux[x, y, z] == 0, uy[x, y, z] == 0,
uz[x, y, z] == 0},
x == -20 || x == 20 || y == -10 || y == 10 || z == -10 ||
z == 10];


Also, note that I set bmx manually. Your bmx evaluates to {0,0,0};

mx = If[ElementMarker == 1, 1, 0];
bmx[x_, y_, z_] := Evaluate[{mx, 0, 0}]
bmx[x, y, z]
(* {If[ElementMarker == 1, 1, 0], 0, 0} *)


Solve:

{Ax, Ay, Az} =
NDSolveValue[{bcs, pde == -bmx[x, y, z]}, {ux, uy,
uz}, {x, y, z} \[Element] mesh];


Visualize:

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],
StreamPlot[{Ay[0, y, z], Az[0, y, z]}, {y, -5, 5}, {z, -5, 5}]}


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


• Great job, thank you very much. Looking at the simplicity of your solution I feel somehow stupid, that I could not figure out the problem earlier ... Anyway, do you have an Idea how to include the Curl[] from the magnetization in a working way? Without the Curl[], the solution will not match the natural behavior of a permanent magnet. In the VectorPlot3D[] graphic we can see, that the B-Field lines are not the lines of a magnet with magnetization of {1,0,0}. It looks more like the field a wire which has a current passed through in {1,0,0} direction. Sep 5, 2022 at 12:49
• I tried solve the curl symbolically and implement the solved form to the magnetization Term, but I had no success until now Sep 5, 2022 at 12:54
• Maybe I have to switch to a PDE using magnetic scalar potential Sep 5, 2022 at 19:47
• @Tschibi, the curl of something like this Curl[{f[x, y, z], 0, 0}, {x, y, z}] is {0, Derivative[0, 0, 1][f][x, y, z], -Derivative[0, 1, 0][f][x, y, z]}. Now, if have D[If[ElementMarker == 1, 1, 0], y] we get (correctly) 0. If you have something like D[If[ElementMarker == 1, y, 0], y] we get If[ElementMarker == 1, 1, 0]`. So I am having trouble understanding what you would like to do. Could you clarify this a bit (in your post, preferably) Sep 6, 2022 at 8:16
• Please find my answer in Update 4. If there is a way to chat directly we feel free to invite me. Sep 6, 2022 at 9:04