# 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 == #[],
"MeshElement" -> "MeshElements",
"ElementMeshDirective" ->
Directive[EdgeForm[], FaceForm[#[]]]]] &, {{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[], B[]} /. 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

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

Create the boundary mesh:

Needs["NDSolveFEM"]
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[], B[]} /. 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

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

(*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][], B[x, y, 0][]}, {x, -20, 20}, {y, -10, 10}, VectorPoints -> Fine, AspectRatio -> Automatic]
VectorPlot[{B[x, 0, z][], B[x, 0, z][]}, {x, -20, 20}, {z, -10, 10}, VectorPoints -> Fine, AspectRatio -> Automatic]
VectorPlot[{B[0, y, z][], B[0, y, z][]}, {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"]
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[#][],         ElementMarker == regionMarker[#]} & /@ regionNames, {0, 0, 0}]),
(Piecewise[{regionsMagnetization[#][], ElementMarker == regionMarker[#]} & /@ regionNames, {0, 0, 0}]),
(Piecewise[{regionsMagnetization[#][], 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"]

(*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[#][],         ElementMarker == regionMarker[#]} & /@ regionNames, {0, 0, 0}]),
(Piecewise[{regionsMagnetization[#][], ElementMarker == regionMarker[#]} & /@ regionNames, {0, 0, 0}]),
(Piecewise[{regionsMagnetization[#][], 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