Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .pre-commit-config.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ repos:
hooks:
- id: isort
name: isort (python)
args: ["--profile", "black", "--filter-files"]

- repo: https://github.com/asottile/pyupgrade
rev: v3.15.2
Expand Down
39 changes: 23 additions & 16 deletions docs/examples/cuboids_demagnetization.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ jupytext:
extension: .md
format_name: myst
format_version: 0.13
jupytext_version: 1.14.5
jupytext_version: 1.16.1
kernelspec:
display_name: Python 3 (ipykernel)
language: python
Expand All @@ -28,32 +28,41 @@ import numpy as np
import pandas as pd
import plotly.express as px
from loguru import logger
from magpylib_material_response.data import get_dataset
from magpylib_material_response.demag import apply_demag
from magpylib_material_response.meshing import mesh_all

if magpy.__version__.split(".")[0] != "5":
raise RuntimeError(
f"Magpylib version must be >=5, (installed: {magpy.__version__})"
)

magpy.defaults.display.backend = "plotly"

# some low quality magnets with different susceptibilities
cube1 = magpy.magnet.Cuboid(magnetization=(0, 0, 1000), dimension=(1, 1, 1))
cube1.move((-1.5, 0, 0))
cube1 = magpy.magnet.Cuboid(polarization=(0, 0, 1), dimension=(0.001, 0.001, 0.001))
cube1.move((-0.0015, 0, 0))
cube1.xi = 0.3 # µr=1.3
cube1.style.label = f"Cuboid, xi={cube1.xi}"

cube2 = magpy.magnet.Cuboid(magnetization=(900, 0, 0), dimension=(1, 1, 1))
cube2.rotate_from_angax(-45, "y").move((0, 0, 0.2))
cube2 = magpy.magnet.Cuboid(polarization=(0.9, 0, 0), dimension=(0.001, 0.001, 0.001))
cube2.rotate_from_angax(-45, "y").move((0, 0, 0.0002))
cube2.xi = 1.0 # µr=2.0
cube2.style.label = f"Cuboid, xi={cube2.xi}"

mx, my = 600 * np.sin(30 / 180 * np.pi), 600 * np.cos(30 / 180 * np.pi)
cube3 = magpy.magnet.Cuboid(magnetization=(mx, my, 0), dimension=(1, 1, 2))
cube3.move((1.6, 0, 0.5)).rotate_from_angax(30, "z")
mx = 0.6 * np.sin(np.deg2rad(30))
my = 0.6 * np.cos(np.deg2rad(30))
cube3 = magpy.magnet.Cuboid(polarization=(mx, my, 0), dimension=(0.001, 0.001, 0.002))
cube3.move((0.0016, 0, 0.0005)).rotate_from_angax(30, "z")
cube3.xi = 0.5 # µr=1.5
cube3.style.label = f"Cuboid, xi={cube3.xi}"

# collection of all cells
coll = magpy.Collection(cube1, cube2, cube3, style_label="No demag")

sensor = magpy.Sensor(position=np.linspace((-4, 0, -1), (4, 0, -1), 301))
sensor = magpy.Sensor(
position=np.linspace((-0.004, 0, -0.001), (0.004, 0, -0.001), 301)
)

magpy.show(*coll, sensor)
```
Expand Down Expand Up @@ -107,8 +116,6 @@ def get_magpylib_dataframe(collection, sensors):
return df


from magpylib_material_response.data import get_dataset

sim_ANSYS = get_dataset("FEMdata_test_cuboids")

df = pd.concat(
Expand All @@ -118,20 +125,20 @@ df = pd.concat(
]
).sort_values(["computation", "path"])

df["Distance [mm]"] = sensor.position[df["path"]][:, 0]
df["Distance [mm]"] -= df["Distance [mm]"].min()
df["Distance [m]"] = sensor.position[df["path"]][:, 0]
df["Distance [m]"] -= df["Distance [m]"].min()
```

```{code-cell} ipython3
px_kwargs = dict(
x="Distance [mm]",
x="Distance [m]",
y=B_cols,
facet_col="variable",
color="computation",
line_dash="computation",
height=400,
facet_col_spacing=0.05,
labels={Bk: f"{Bk} [mT]" for Bk in B_cols},
labels={**{Bk: f"{Bk} [T]" for Bk in B_cols}, "value": "value [T]"},
)
fig1 = px.line(
df,
Expand All @@ -157,4 +164,4 @@ fig2.update_yaxes(matches=None, showticklabels=True)
display(fig1, fig2)
```

As shown above, already with a low number of mesh elements, the result is approaching the reference FEM values.
As shown above, already with a low number of mesh elements, the result is approaching the reference FEM values and improves while refining the mesh.
24 changes: 12 additions & 12 deletions docs/examples/soft_magnets.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ jupytext:
extension: .md
format_name: myst
format_version: 0.13
jupytext_version: 1.14.5
jupytext_version: 1.16.1
kernelspec:
display_name: Python 3 (ipykernel)
language: python
Expand Down Expand Up @@ -39,14 +39,14 @@ from magpylib_material_response.meshing import mesh_all
magpy.defaults.display.backend = "plotly"

# hard magnet
cube1 = magpy.magnet.Cuboid(magnetization=(0, 0, 1000), dimension=(1, 1, 2))
cube1.move((0, 0, 0.5))
cube1 = magpy.magnet.Cuboid(polarization=(0, 0, 1), dimension=(0.001, 0.001, 0.002))
cube1.move((0, 0, 0.0005))
cube1.xi = 0.5 # µr=1.5
cube1.style.label = f"Hard cuboid magnet, xi={cube1.xi}"

# soft magnet
cube2 = magpy.magnet.Cuboid(magnetization=(0, 0, 0), dimension=(1, 1, 1))
cube2.rotate_from_angax(angle=45, axis="y").move((1.5, 0, 0))
cube2 = magpy.magnet.Cuboid(polarization=(0, 0, 0), dimension=(0.001, 0.001, 0.001))
cube2.rotate_from_angax(angle=45, axis="y").move((0.0015, 0, 0))
cube2.xi = 3999 # µr=4000
cube2.style.label = f"Soft cuboid magnet, xi={cube2.xi}"

Expand All @@ -56,10 +56,10 @@ coll = magpy.Collection(cube1, cube2, style_label="No demag")
# add sensors
sensors = [
magpy.Sensor(
position=np.linspace((-4, 0, z), (6, 0, z), 1001),
style_label=f"Sensor, z={z}mm",
position=np.linspace((-0.004, 0, z), (0.006, 0, z), 1001),
style_label=f"Sensor, z={z}m",
)
for z in (-1, -3, -5)
for z in (-0.001, -0.003, -0.005)
]

magpy.show(*coll, *sensors)
Expand Down Expand Up @@ -130,8 +130,8 @@ df = pd.concat(
).sort_values(["computation", "path"])


df["Distance [mm]"] = sensors[0].position[df["path"]][:, 0]
df["Distance [mm]"] -= df["Distance [mm]"].min()
df["Distance [m]"] = sensors[0].position[df["path"]][:, 0]
df["Distance [m]"] -= df["Distance [m]"].min()
```

```{code-cell} ipython3
Expand All @@ -144,7 +144,7 @@ px_kwargs = dict(
line_dash="computation",
height=600,
facet_col_spacing=0.05,
labels={Bk: f"{Bk} [mT]" for Bk in B_cols},
labels={**{Bk: f"{Bk} [T]" for Bk in B_cols}, "value": "value [T]"},
)
fig1 = px.line(
df,
Expand Down Expand Up @@ -172,4 +172,4 @@ display(fig1, fig2)

+++ {"user_expressions": []}

As shown above, the demagnetized collection outputs are approaching the reference FEM values.
As shown above, the demagnetized collection outputs are approaching the reference FEM values while refining the mesh.
64 changes: 32 additions & 32 deletions magpylib_material_response/demag.py
Original file line number Diff line number Diff line change
Expand Up @@ -105,44 +105,44 @@ def demag_tensor(
rot0 = R.from_quat(rotQ0)

H_point = []
for unit_mag in [(1, 0, 0), (0, 1, 0), (0, 0, 1)]:
mag_all = rot0.inv().apply(unit_mag)
for unit_pol in [(1, 0, 0), (0, 1, 0), (0, 0, 1)]:
pol_all = rot0.inv().apply(unit_pol)
# point matching field and demag tensor
with timelog(f"getH with unit_mag={unit_mag}", min_log_time=min_log_time):
with timelog(f"getH with unit_pol={unit_pol}", min_log_time=min_log_time):
if pairs_matching or max_dist != 0:
magnetization = np.repeat(mag_all, len(src_list), axis=0)[mask_inds]
polarization = np.repeat(pol_all, len(src_list), axis=0)[mask_inds]
H_unique = magpy.getH(
"Cuboid", magnetization=magnetization, **getH_params
"Cuboid", polarization=polarization, **getH_params
)
if max_dist != 0:
H_temp = np.zeros((len(src_list) ** 2, 3))
H_temp[mask_inds] = H_unique
H_unit_mag = H_temp
H_unit_pol = H_temp
else:
H_unit_mag = H_unique[unique_inv_inds]
H_unit_pol = H_unique[unique_inv_inds]
else:
for src, mag in zip(src_list, mag_all):
src.magnetization = mag
for src, pol in zip(src_list, pol_all):
src.polarization = pol
if split > 1:
src_list_split = np.array_split(src_list, split)
with logger.contextualize(
task="Splitting field calculation", split=split
):
H_unit_mag = []
H_unit_pol = []
for split_ind, src_list_subset in enumerate(src_list_split):
logger.info(
f"Sources subset {split_ind+1}/{len(src_list_split)}"
)
if src_list_subset.size > 0:
H_unit_mag.append(
H_unit_pol.append(
magpy.getH(src_list_subset.tolist(), pos0)
)
H_unit_mag = np.concatenate(H_unit_mag, axis=0)
H_unit_pol = np.concatenate(H_unit_pol, axis=0)
else:
H_unit_mag = magpy.getH(src_list, pos0)
H_point.append(H_unit_mag) # shape (n_cells, n_pos, 3_xyz)
H_unit_pol = magpy.getH(src_list, pos0)
H_point.append(H_unit_pol) # shape (n_cells, n_pos, 3_xyz)

# shape (3_unit_mag, n_cells, n_pos, 3_xyz)
# shape (3_unit_pol, n_cells, n_pos, 3_xyz)
T = np.array(H_point).reshape((3, nof_src, nof_src, 3))

return T
Expand Down Expand Up @@ -255,7 +255,7 @@ def apply_demag(
):
"""
Computes the interaction between all collection magnets and fixes their
magnetization.
polarization.

Parameters
----------
Expand Down Expand Up @@ -332,11 +332,11 @@ def apply_demag(
)
with timelog(demag_msg, min_log_time=min_log_time):
# set up mr
mag_magnets = [
src.orientation.apply(src.magnetization) for src in magnets_list
pol_magnets = [
src.orientation.apply(src.polarization) for src in magnets_list
] # ROTATION CHECK
mag_magnets = np.reshape(
mag_magnets, (3 * n, 1), order="F"
pol_magnets = np.reshape(
pol_magnets, (3 * n, 1), order="F"
) # shape ii = x1, ... xn, y1, ... yn, z1, ... zn

# set up S
Expand All @@ -349,7 +349,7 @@ def apply_demag(
)
S = np.diag(np.tile(xi, 3)) # shape ii, jj

# set up T (3 mag unit, n cells, n positions, 3 Bxyz)
# set up T (3 pol unit, n cells, n positions, 3 Bxyz)
with timelog("Demagnetization tensor calculation", min_log_time=min_log_time):
T = demag_tensor(
magnets_list,
Expand All @@ -358,32 +358,32 @@ def apply_demag(
max_dist=max_dist,
)

T *= 4 * np.pi / 10
T *= magpy.mu_0
T = T.swapaxes(2, 3).reshape((3 * n, 3 * n)).T # shape ii, jj

mag_tolal = mag_magnets
pol_tolal = pol_magnets

if currents_list:
with timelog(
"Add current sources contributions", min_log_time=min_log_time
):
pos = np.array([src.position for src in magnets_list])
mag_currents = magpy.getB(currents_list, pos, sumup=True)
mag_currents = np.reshape(mag_currents, (3 * n, 1), order="F")
mag_tolal += np.matmul(S, mag_currents)
pol_currents = magpy.getB(currents_list, pos, sumup=True)
pol_currents = np.reshape(pol_currents, (3 * n, 1), order="F")
pol_tolal += np.matmul(S, pol_currents)

# set up Q
Q = np.eye(3 * n) - np.matmul(S, T)

# determine new magnetization vectors
# determine new polarization vectors
with timelog("Solving of linear system", min_log_time=1):
mag_new = np.linalg.solve(Q, mag_tolal)
pol_new = np.linalg.solve(Q, pol_tolal)

mag_new = np.reshape(mag_new, (n, 3), order="F")
# mag_new *= .4*np.pi
pol_new = np.reshape(pol_new, (n, 3), order="F")
# pol_new *= .4*np.pi

for s, mag in zip(collection.sources_all, mag_new):
s.magnetization = s.orientation.inv().apply(mag) # ROTATION CHECK
for s, pol in zip(collection.sources_all, pol_new):
s.polarization = s.orientation.inv().apply(pol) # ROTATION CHECK

if not inplace:
return collection
Loading