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["NDSolve`FEM`"]
(*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 = Simplify`PWToUnitStep@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 = Simplify`PWToUnitStep@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
Thanks in advance.
Update: I tried another approach by generating the mesh with OpencascadeLink. Unfortunately the problem persists.
Needs["NDSolve`FEM`"]
Needs["OpenCascadeLink`"]
core = OpenCascadeShapeImport[
FileNameJoin[{NotebookDirectory[], "Magnet.stl"}]];
air = OpenCascadeShape[Cuboid[{-20, -10, -10}, {20, 10, 10}]];
intersection = OpenCascadeShapeIntersection[core, air];
shape = OpenCascadeShapeSewing[{air, intersection}];
(*Create Boundary Mesh*)
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[shape];
(*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["NDSolve`FEM`"]
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 = Simplify`PWToUnitStep@
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 = Simplify`PWToUnitStep@
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.
Truewill do that, but perhaps you only want the bc at the outer faces of the cube? See this $\endgroup$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}]` $\endgroup$