How to simulate permanent magnets using the PDE Toolbox in R2025a

I am trying to model and visualize two permanent magnet spheres using the below code snippet:

function twoMagnets_femodel_demo()
clc; clear; close all;
% -------------------------------
% 1) Define geometry parameters
% -------------------------------
% Let's assume "1/2-inch magnets" => diameter = 0.5 inch => radius = 0.25 in
% 1 inch ~ 0.0254 m, so 0.25 in ~ 0.00635 m
R_sphere = 0.00635;   % [m]

% Bounding box of ~5 cm each side
boxDim = 0.05;        % [m]
% -------------------------------
% 2) Build geometry parts
% -------------------------------
% (a) Air region: bounding box
boxGm    = multicuboid(boxDim, boxDim, boxDim);
% (b) Sphere #1
sphere1Gm = multisphere(R_sphere);
% (c) Sphere #2
sphere2Gm = multisphere(R_sphere);
% Translate spheres so that they're centered in the box, side-by-side
% Let the box corners go from (0,0,0) to (0.05,0.05,0.05), so center ~ (0.025,0.025,0.025)
centerBox = [0.025, 0.025, 0.025];
offset1   = centerBox + [-R_sphere, 0, 0];  % shift left in x
offset2   = centerBox + [ R_sphere, 0, 0];  % shift right in x
sphere1Gm = translate(sphere1Gm, offset1);
sphere2Gm = translate(sphere2Gm, offset2);
% -------------------------------
% 3) Combine geometries
% -------------------------------
% addCell(...) is the new PDE Toolbox function that 
% keeps each shape as a separate "cell" in the final geometry.
gm = addCell(boxGm, sphere1Gm);
gm = addCell(gm,   sphere2Gm);
% Visualize & label
figure;
pdegplot(gm, FaceAlpha=0.2, CellLabels="on");
title("Bounding Box + Two Spheres");
axis equal
% At this point we should see "Cell #1 = bounding box, #2 = first sphere, #3 = second sphere"
% -------------------------------
% 4) Create femodel for magnetostatics
% -------------------------------
model = femodel(AnalysisType="magnetostatic", Geometry=gm);
% Specify vacuum permeability in SI (H/m)
model.VacuumPermeability = 4*pi*1e-7;    % ~ 1.2566370614e-6
% -------------------------------
% 5) Set material properties
% -------------------------------
% We have 3 cells:
%   1 => bounding box (air)
%   2 => sphere #1
%   3 => sphere #2
%
% For air and the spheres, assume relative permeability ~ 1
model.MaterialProperties(1) = materialProperties(RelativePermeability=1);
model.MaterialProperties(2) = materialProperties(RelativePermeability=1);
model.MaterialProperties(3) = materialProperties(RelativePermeability=1);
% -------------------------------
% 6) Assign magnetization
% -------------------------------
% In new PDE Toolbox, you can set magnetization on the "CellLoad" for a subdomain (cell).
% Suppose sphere #1 is magnetized in +z; sphere #2 in -z.
% For a typical neodymium magnet, you might set M ~ 1e5 or 1e6 A/m. 
M1 = [0; 0; 1e5];
M2 = [0; 0;-1e5];
model.CellLoad(2) = cellLoad(Magnetization=M1);   % sphere1
model.CellLoad(3) = cellLoad(Magnetization=M2);   % sphere2
% -------------------------------
% 7) Boundary conditions
% -------------------------------
% Typically, you'd set a "zero" magnetic potential (A = 0) on the outer faces 
% of the bounding box to approximate an open boundary. 
% Identify the face IDs of the bounding box:
% pdegplot(gm, FaceAlpha=0.1, CellLabels="on", FaceLabels="on")
% The bounding box typically has 6 faces with IDs [1..6]. 
% So we do:
model.FaceBC(1:6) = faceBC(MagneticPotential=[0;0;0]);
% (If you find different face numbering, adjust accordingly.)
% -------------------------------
% 8) Generate mesh
% -------------------------------
% We can specify "Hcell" to refine the mesh in certain cells, e.g. the magnets,
% or keep it uniform. For example, refine in cells 2 & 3 (the spheres):
internalCells = [2 3];
model = generateMesh(model, Hcell={internalCells, R_sphere/3});
% -------------------------------
% 9) Solve
% -------------------------------
R = solve(model);
% R is a MagnetostaticResults object with fields:
%   R.MagneticPotential, R.MagneticField, R.MagneticFluxDensity, R.Mesh, etc.
% -------------------------------
% 10) Postprocessing
% -------------------------------
% Extract Bx, By, Bz at the mesh nodes
Bx   = R.MagneticFluxDensity.Bx;
By   = R.MagneticFluxDensity.By;
Bz   = R.MagneticFluxDensity.Bz;
Bmag = sqrt(Bx.^2 + By.^2 + Bz.^2);
% Simple 3D plot of |B|
figure;
pdeplot3D(R.Mesh, ColorMapData=Bmag);
title("|B| flux density magnitude");
colorbar;
% Optionally, a vector plot:
figure;
pdeplot3D(R.Mesh, FlowData=[Bx By Bz]);
title("Magnetic Field Vectors");
end

While logically complete, the above code doesn't work. The current error I am struggling with is that addCell() only works inside the original geometry. And throws the following error here:

Error using pde.DiscreteGeometry/addCell

Added cells must be located strictly inside a single cell of the original geometry.

Error in

TwoMagnets (line 38)

gm = addCell(boxGm, sphere1Gm);

^^^^^^^^^^^^^^^^^^^^^^^^^

Is there an example script on how to simulate permanent magnets using the latest PDE Toolbox?

Answer 1

[Shivam Gothi]

Hello @Philippe,

I investigated the code attached by you. I found the cause behind the error. Let me explain you in step by step manner:

I will first show the geometries defined by you. Consider first part of your code:

clc; clear; close all;
% -------------------------------
% 1) Define geometry parameters
% -------------------------------
% Let's assume "1/2-inch magnets" => diameter = 0.5 inch => radius = 0.25 in
% 1 inch ~ 0.0254 m, so 0.25 in ~ 0.00635 m
R_sphere = 0.00635;   % [m]
% Bounding box of ~5 cm each side
boxDim = 0.05;        % [m]
% -------------------------------
% 2) Build geometry parts
% -------------------------------
% (a) Air region: bounding box
boxGm    = multicuboid(boxDim, boxDim, boxDim);
% (b) Sphere #1
sphere1Gm = multisphere(R_sphere);
% (c) Sphere #2
sphere2Gm = multisphere(R_sphere);
% Translate spheres so that they're centered in the box, side-by-side
% Let the box corners go from (0,0,0) to (0.05,0.05,0.05), so center ~ (0.025,0.025,0.025)
centerBox = [0.025, 0.025, 0.025];
offset1   = centerBox + [-R_sphere, 0, 0];  % shift left in x
offset2   = centerBox + [ R_sphere, 0, 0];  % shift right in x
sphere1Gm = translate(sphere1Gm, offset1);
sphere2Gm = translate(sphere2Gm, offset2);
% -------------------------------
% 3) Combine geometries
% -------------------------------
% addCell(...) is the new PDE Toolbox function that 
% keeps each shape as a separate "cell" in the final geometry.
% Visualize & label
figure;
pdegplot(sphere1Gm, FaceAlpha=0.2, CellLabels="on");
hold on;
pdegplot(sphere2Gm, FaceAlpha=0.2, CellLabels="on");
hold on;
pdegplot(boxGm, FaceAlpha=0.2, CellLabels="on");
title("Bounding Box + Two Spheres");
axis equal

From the above figure, it is clearly seen that the spheres are not lying inside the "boxGm". Therefore, you are recieving the error:

"Added cells must be located strictly inside a single cell of the original geometry."

I corrected the position of spheres so that they lie inside the box as shown below. I corrected "centerBox".

% -------------------------------
% 1) Define geometry parameters
% -------------------------------
% Let's assume "1/2-inch magnets" => diameter = 0.5 inch => radius = 0.25 in
% 1 inch ~ 0.0254 m, so 0.25 in ~ 0.00635 m
R_sphere = 0.00635;   % [m]
% Bounding box of ~5 cm each side
boxDim = 0.05;        % [m]
% -------------------------------
% 2) Build geometry parts
% -------------------------------
% (a) Air region: bounding box
boxGm    = multicuboid(boxDim, boxDim, boxDim);
% (b) Sphere #1
sphere1Gm = multisphere(R_sphere);
% (c) Sphere #2
sphere2Gm = multisphere(R_sphere);
% Translate spheres so that they're centered in the box, side-by-side
% Let the box corners go from (0,0,0) to (0.05,0.05,0.05), so center ~ (0.025,0.025,0.025)
centerBox = [0, 0, 0.025];
offset1   = centerBox + [-R_sphere, 0, 0];  % shift left in x
offset2   = centerBox + [ R_sphere, 0, 0];  % shift right in x
sphere1Gm = translate(sphere1Gm, offset1);
sphere2Gm = translate(sphere2Gm, offset2);
% -------------------------------
% 3) Combine geometries
% -------------------------------
% addCell(...) is the new PDE Toolbox function that 
% keeps each shape as a separate "cell" in the final geometry.
% Visualize & label
figure;
pdegplot(sphere1Gm, FaceAlpha=0.2, CellLabels="on");
hold on;
pdegplot(sphere2Gm, FaceAlpha=0.2, CellLabels="on");
hold on;
pdegplot(boxGm, FaceAlpha=0.2, CellLabels="on");
title("Bounding Box + Two Spheres");
axis equal

You can see that now the spheres are properly positioned inside the box.

Now, I executed the entire code with the corrected value of "centerBox". Also, I modified the command :

model = generateMesh(model, Hcell={internalCells, R_sphere/3});

to:

model = generateMesh(model);

because it was throwing error. In order to demonstrate the correctness of code and illustrate it in simple manner, I have simply generated the mesh without specifing any additional properties. You can generate mesh as per you requirements by referring to the documentation for generateMesh function: https://www.mathworks.com/help/pde/ug/pde.pdemodel.generatemesh.html

I am attaching the modified code below:

% -------------------------------
% 1) Define geometry parameters
% -------------------------------
% Let's assume "1/2-inch magnets" => diameter = 0.5 inch => radius = 0.25 in
% 1 inch ~ 0.0254 m, so 0.25 in ~ 0.00635 m
R_sphere = 0.00635;   % [m]
% Bounding box of ~5 cm each side
boxDim = 0.05;        % [m]
% -------------------------------
% 2) Build geometry parts
% -------------------------------
% (a) Air region: bounding box
boxGm    = multicuboid(boxDim, boxDim, boxDim);
% (b) Sphere #1
sphere1Gm = multisphere(R_sphere);
% (c) Sphere #2
sphere2Gm = multisphere(R_sphere);
% Translate spheres so that they're centered in the box, side-by-side
% Let the box corners go from (0,0,0) to (0.05,0.05,0.05), so center ~ (0.025,0.025,0.025)
centerBox = [0, 0, 0.025];
offset1   = centerBox + [-R_sphere, 0, 0];  % shift left in x
offset2   = centerBox + [ R_sphere, 0, 0];  % shift right in x
sphere1Gm = translate(sphere1Gm, offset1);
sphere2Gm = translate(sphere2Gm, offset2);
% -------------------------------
% 3) Combine geometries
% -------------------------------
% addCell(...) is the new PDE Toolbox function that 
% keeps each shape as a separate "cell" in the final geometry.
gm = addCell(boxGm, sphere1Gm);
gm = addCell(gm,   sphere2Gm);
% Visualize & label
figure;
pdegplot(gm, FaceAlpha=0.2, CellLabels="on");
title("Bounding Box + Two Spheres");
axis equal

% At this point we should see "Cell #1 = bounding box, #2 = first sphere, #3 = second sphere"
% -------------------------------
% 4) Create femodel for magnetostatics
% -------------------------------
model = femodel(AnalysisType="magnetostatic", Geometry=gm);
% Specify vacuum permeability in SI (H/m)
model.VacuumPermeability = 4*pi*1e-7;    % ~ 1.2566370614e-6
% -------------------------------
% 5) Set material properties
% -------------------------------
% We have 3 cells:
%   1 => bounding box (air)
%   2 => sphere #1
%   3 => sphere #2
%
% For air and the spheres, assume relative permeability ~ 1
model.MaterialProperties(1) = materialProperties(RelativePermeability=1);
model.MaterialProperties(2) = materialProperties(RelativePermeability=1);
model.MaterialProperties(3) = materialProperties(RelativePermeability=1);
% -------------------------------
% 6) Assign magnetization
% -------------------------------
% In new PDE Toolbox, you can set magnetization on the "CellLoad" for a subdomain (cell).
% Suppose sphere #1 is magnetized in +z; sphere #2 in -z.
% For a typical neodymium magnet, you might set M ~ 1e5 or 1e6 A/m. 
M1 = [0; 0; 1e5];
M2 = [0; 0;-1e5];
model.CellLoad(2) = cellLoad(Magnetization=M1);   % sphere1
model.CellLoad(3) = cellLoad(Magnetization=M2);   % sphere2
% -------------------------------
% 7) Boundary conditions
% -------------------------------
% Typically, you'd set a "zero" magnetic potential (A = 0) on the outer faces 
% of the bounding box to approximate an open boundary. 
% Identify the face IDs of the bounding box:
% pdegplot(gm, FaceAlpha=0.1, CellLabels="on", FaceLabels="on")
% The bounding box typically has 6 faces with IDs [1..6]. 
% So we do:
model.FaceBC(1:6) = faceBC(MagneticPotential=[0;0;0]);
% (If you find different face numbering, adjust accordingly.)
% -------------------------------
% 8) Generate mesh
% -------------------------------
% We can specify "Hcell" to refine the mesh in certain cells, e.g. the magnets,
% or keep it uniform. For example, refine in cells 2 & 3 (the spheres):
internalCells = [2 3];
model = generateMesh(model);
% -------------------------------
% 9) Solve
% -------------------------------
R = solve(model);
% R is a MagnetostaticResults object with fields:
%   R.MagneticPotential, R.MagneticField, R.MagneticFluxDensity, R.Mesh, etc.
% -------------------------------
% 10) Postprocessing
% -------------------------------
% Extract Bx, By, Bz at the mesh nodes
Bx   = R.MagneticFluxDensity.Bx;
By   = R.MagneticFluxDensity.By;
Bz   = R.MagneticFluxDensity.Bz;
Bmag = sqrt(Bx.^2 + By.^2 + Bz.^2);
% Simple 3D plot of |B|
figure;
pdeplot3D(R.Mesh, ColorMapData=Bmag);
title("|B| flux density magnitude");
colorbar;

% Optionally, a vector plot:
figure;
pdeplot3D(R.Mesh, FlowData=[Bx By Bz]);
title("Magnetic Field Vectors");

I hope this helps !

[Philippe]

Thank you so much @Shivam Gothi for taking the time to respond in such a thoughtful and structured manner. I noticed that the locations were off, but I didn't realize that they were the cause of this.

[Philippe]

PS: I assume that since GenerateMesh() doesn't offer an 'HCell' parameter, would the correct way of refining the mesh near the spheres be to instead identify the face and use 'HFace'?

[Arne]

Thank you for laying this out so nicely! Would you happen to know what units MagneticFluxDensity is in? I can't seem to find this information anywhere.

[Shivam Gothi]

Hello @Arne,

The SI unit of Magnetic Flux Density is tesla (Also called as weber per meter square). I also tried to find out the unit of "MagneticFluxDensity", but I was not able to find it.

But I was going through the below documentation of "MagnetostaticResults" object:

https://www.mathworks.com/help/pde/ug/pde.magnetostaticresults.html#mw_749e161b-f4f9-4209-87b6-279ea607f7d8

below given is the snapshot from the docmentation:

the highlighted equation is standard equation of magnetostatics, which gives the value of magnetic flux density "B" in Tesla. "MagneticFluxDensity" is a property of object"MagnetostaticResults"

Therefore, as per my point of view, the unit of " MagneticFluxDensity" is tesla. But I am not 100% sure.

[Shivam Gothi]

Hello @Philippe,

I am not much familiar with Mesh structure and formation. May be, you can try to understand it by the means of examples stated in the documentaion:

https://www.mathworks.com/help/pde/ug/pde.pdemodel.generatemesh.html

If you still have questions, I recommend you to ask a new question in MATLAB answer community with proper title and MATLAB answer tags. The community users, who are interested in magnetostatic domain, will be able to help you.

[Arne]

Thank you @Shivam Gothi! I think you are right! I have tried calibrating a very simple model of a cube magnet with some published results and it all makes sense. Magnetization is in A/m and MagneticFluxDensity is in T.