update to magpylib v5 (#4)

* update to magpylib v5

* fix tests

* fix isort black conflict

---------

Co-authored-by: Boisselet Alexandre (IFAT DC ATV SC D TE2) <[email protected]>

.pre-commit-config.yaml

@@ -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

docs/examples/cuboids_demagnetization.md

@@ -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
@@ -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)
 ```
@@ -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(
@@ -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,
@@ -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.

docs/examples/soft_magnets.md

@@ -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
@@ -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}"
 
@@ -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)
@@ -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
@@ -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,
@@ -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.

magpylib_material_response/demag.py

@@ -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
@@ -255,7 +255,7 @@ def apply_demag(
 ):
     """
     Computes the interaction between all collection magnets and fixes their
-    magnetization.
+    polarization.
 
     Parameters
     ----------
@@ -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
@@ -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,
@@ -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