Merge pull request #15 from magpylib/dipole_force

Dipole force

Author: fslanovc

CHANGELOG.md

@@ -1,4 +1,8 @@
-# Unpublished
+# v0.4.0 (2025/02/03)
+
+- Implementation of Dipole as target
+- Bugfix for magnets with meshing parameter = 1
+- Additional tests for Dipole extension
 
 # v0.3.1 (2024/10/30)
 

src/magpylib_force/force.py

@@ -14,6 +14,7 @@
 from magpylib._src.obj_classes.class_magnet_Cylinder import Cylinder
 from magpylib._src.obj_classes.class_magnet_CylinderSegment import CylinderSegment
 from magpylib._src.obj_classes.class_magnet_Sphere import Sphere
+from magpylib._src.obj_classes.class_misc_Dipole import Dipole
 
 from magpylib_force.meshing import mesh_target
 from magpylib_force.utility import check_input_anchor, check_input_targets
@@ -35,8 +36,8 @@ def getFT(sources, targets, anchor=None, eps=1e-5, squeeze=True):
         or a 1D list of l sources and/or collection objects.
 
     targets: single object or 1D list of t objects that are Sphere, Cuboid, Polyline,
-        Cylinder, or CylinderSegment. Force and Torque acting on targets in the magnetic
-        field generated by the sources will be computed. A target must have a valid
+        Cylinder, CylinderSegment, or Dipole. Force and Torque acting on targets in the magnetic
+        field generated by the sources will be computed. A target (except of Dipole) must have a valid
         `meshing` parameter.
 
     anchor: array_like, default=None
@@ -65,14 +66,15 @@ def getFT(sources, targets, anchor=None, eps=1e-5, squeeze=True):
     n = len(targets)
 
     # split targets into lists of similar types
-    TARGET_TYPES = [Cuboid, Polyline, Sphere, Cylinder, CylinderSegment, Circle]
+    TARGET_TYPES = [Cuboid, Polyline, Sphere, Cylinder, CylinderSegment, Circle, Dipole]
     getFT_FUNCS = [
         getFTmagnet,
         getFTcurrent,
         getFTmagnet,
         getFTmagnet,
         getFTmagnet,
         getFTcurrent_circ,
+        getFTdipole,
     ]
     objects = [[] for _ in TARGET_TYPES]
     orders = [[] for _ in TARGET_TYPES]
@@ -202,6 +204,8 @@ def getFTmagnet(sources, targets, eps=1e-5, anchor=None):
             POSS[insti[i] : insti[i + 1], j] = mesh + ev + tgt.position
 
     BB = magpy.getB(sources, POSS, sumup=True)
+    if BB.ndim == 2:
+        BB = np.expand_dims(BB, axis=0)
     gradB = (BB[:, 1::2] - BB[:, 2::2]) / (2 * eps)
     gradB = np.swapaxes(gradB, 0, 1)
 
@@ -333,3 +337,71 @@ def getFTcurrent(sources, targets, anchor=None, eps=None):  # noqa:  ARG001
     )
 
     return np.array((F, -T))
+
+
+def getFTdipole(sources, targets, anchor=None, eps=None):
+    """
+    Compute force and torque acting on magnetic dipoles
+
+    Parameters
+    ----------
+    sources: source and collection objects or 1D list thereof
+        Sources that generate the magnetic field. Can be a single source (or collection)
+        or a 1D list of l sources and/or collection objects.
+
+    targets: Dipole object or 1D list of Dipole objects
+        Force and Torque acting on targets in the magnetic field generated by the sources
+        will be computed.
+
+    eps: float, default=1e-5
+        The magnetic field gradient is computed using finite differences (FD). eps is
+        the FD step size. A good value is 1e-5 * characteristic system size (magnet size,
+        distance between sources and targets, ...).
+
+    anchor: array_like, default=None
+        The Force adds to the Torque via the anchor point. For a freely floating magnet
+        this would be the barycenter. If `anchor=None`, this part of the Torque computation
+        is omitted.
+    """
+    # number of magnets
+    tgt_number = len(targets)
+
+    # field computation positions (1xfor B, 6x for gradB)
+    POSS = np.zeros((tgt_number, 7, 3))
+
+    # moment of each instance
+    MOM = np.zeros((tgt_number, 3))
+
+    # MISSING: eps should be defined relative to the sizes of the objects
+    eps_vec = np.array(
+        [
+            (0, 0, 0),
+            (eps, 0, 0),
+            (-eps, 0, 0),
+            (0, eps, 0),
+            (0, -eps, 0),
+            (0, 0, eps),
+            (0, 0, -eps),
+        ]
+    )
+
+    for i, tgt in enumerate(targets):
+        dipole_mom = tgt.orientation.apply(tgt.moment)
+        MOM[i] = dipole_mom
+
+        for j, ev in enumerate(eps_vec):
+            POSS[i, j] = ev + tgt.position
+
+    BB = magpy.getB(sources, POSS, sumup=True)
+    if BB.ndim == 2:
+        BB = np.expand_dims(BB, axis=0)
+    gradB = (BB[:, 1::2] - BB[:, 2::2]) / (2 * eps)
+    gradB = np.swapaxes(gradB, 0, 1)
+
+    F = np.sum((gradB * MOM), axis=2).T
+    # Ts = np.zeros((no_inst,3))
+    T = np.cross(BB[:, 0], MOM)
+    if anchor is not None:
+        T -= np.cross(POSS[:, 0] - anchor, F)
+
+    return np.array((F, -T))

src/magpylib_force/utility.py

@@ -13,6 +13,7 @@
 from magpylib._src.obj_classes.class_magnet_Cylinder import Cylinder
 from magpylib._src.obj_classes.class_magnet_CylinderSegment import CylinderSegment
 from magpylib._src.obj_classes.class_magnet_Sphere import Sphere
+from magpylib._src.obj_classes.class_misc_Dipole import Dipole
 
 
 def check_input_anchor(anchor):
@@ -43,7 +44,7 @@ def check_input_targets(targets):
         targets = [targets]
     for t in targets:
         if not isinstance(
-            t, Cuboid | Polyline | Sphere | Cylinder | CylinderSegment | Circle
+            t, Cuboid | Polyline | Sphere | Cylinder | CylinderSegment | Circle | Dipole
         ):
             msg = (
                 "Bad `targets` input for getFT."
@@ -52,16 +53,21 @@ def check_input_targets(targets):
                 f" Instead received type {type(t)} target."
             )
             raise ValueError(msg)
-        if not hasattr(t, "meshing"):
-            msg = (
-                "Missing input for getFT `targets`."
-                " `targets` must have the `meshing` parameter set."
-            )
-            raise ValueError(msg)
-        if not isinstance(t, Polyline) and np.isscalar(t.meshing) and t.meshing < 20:
-            warnings.warn(
-                "Input parameter `meshing` is set to a low value which will result in "
-                "inaccurate computation of force and torque.",
-                stacklevel=2,
-            )
+        if not isinstance(t, Dipole):
+            if not hasattr(t, "meshing"):
+                msg = (
+                    "Missing input for getFT `targets`."
+                    " `targets` must have the `meshing` parameter set."
+                )
+                raise ValueError(msg)
+            if (
+                not isinstance(t, Polyline)
+                and np.isscalar(t.meshing)
+                and t.meshing < 20
+            ):
+                warnings.warn(
+                    "Input parameter `meshing` is set to a low value which will result in "
+                    "inaccurate computation of force and torque.",
+                    stacklevel=2,
+                )
     return targets

tests/test_basics.py

@@ -2,6 +2,7 @@
 
 import magpylib as magpy
 import numpy as np
+from scipy.spatial.transform import Rotation as R
 
 from magpylib_force import getFT
 
@@ -57,3 +58,20 @@ def test_rotation2():
     FT = getFT(s1, [c1, c2], anchor=(0, 0, 0))
 
     np.testing.assert_allclose(FT[0], FT[1])
+
+
+def test_orientation():
+    """
+    test if dipole with orientation gives same result as rotated magnetic moment
+    """
+    mm, md = np.array((0.976, 4.304, 2.055)), np.array((0.878, -1.527, 2.918))
+    pm, pd = np.array((-1.248, 7.835, 9.273)), np.array((-2.331, 5.835, 0.578))
+
+    magnet = magpy.magnet.Cuboid(position=pm, dimension=(1, 2, 3), polarization=mm)
+    r = R.from_euler("xyz", (25, 65, 150), degrees=True)
+    dipole1 = magpy.misc.Dipole(position=pd, moment=md, orientation=r)
+    dipole2 = magpy.misc.Dipole(position=pd, moment=r.apply(md))
+
+    F = getFT(magnet, [dipole1, dipole2], anchor=(0, 0, 0))
+
+    np.testing.assert_allclose(F[0], F[1])

tests/test_physics.py

@@ -242,3 +242,36 @@ def test_physics_force_between_cocentric_loops():
     F_ana = pf * ((2 - k2) / (1 - k2) * ellipe(k**2) - 2 * ellipk(k**2))
 
     assert abs((F_num - F_ana) / (F_num + F_ana)) < 1e-5
+
+
+def test_physics_force_between_two_dipoles():
+    """
+    Compare force between two magnetic Dipoles with well-known formula
+    (e.g. https://en.wikipedia.org/wiki/Magnetic_dipole%E2%80%93dipole_interaction)
+    """
+    m1, m2 = np.array((0.976, 4.304, 2.055)), np.array((0.878, -1.527, 2.918))
+    p1, p2 = np.array((-1.248, 7.835, 9.273)), np.array((-2.331, 5.835, 0.578))
+
+    # numerical solution
+    dipole1 = magpy.misc.Dipole(position=p1, moment=m1)
+    dipole2 = magpy.misc.Dipole(position=p2, moment=m2)
+    F_num = getFT(dipole1, dipole2, anchor=(0, 0, 0))[0, :]
+    # F_num = getFT([dipole1], [dipole2, dipole2])
+
+    # analytical solution
+    r = p2 - p1
+    r_abs = np.linalg.norm(r)
+    r_unit = r / r_abs
+    F_ana = (
+        3
+        * magpy.mu_0
+        / (4 * np.pi * r_abs**4)
+        * (
+            np.cross(np.cross(r_unit, m1), m2)
+            + np.cross(np.cross(r_unit, m2), m1)
+            - 2 * r_unit * np.dot(m1, m2)
+            + 5 * r_unit * (np.dot(np.cross(r_unit, m1), np.cross(r_unit, m2)))
+        )
+    )
+
+    np.testing.assert_allclose(F_num, F_ana, rtol=1e-8, atol=1e-8)

tests/test_self_consistency.py

@@ -264,3 +264,27 @@ def test_consistency_polyline_circle():
     F2, T2 = getFT(src, loop2, anchor=(0, 0, 0))
     assert abs(np.linalg.norm(F1 - F2) / np.linalg.norm(F1 + F2)) < 1e-7
     assert abs(np.linalg.norm(T1 - T2) / np.linalg.norm(T1 + T2)) < 1e-7
+
+
+def test_consistency_sphere_dipole():
+    """
+    force on sphere and dipole should be the same in nearly homogeneous field
+    """
+
+    src = magpy.current.Circle(diameter=10, current=123)
+    pos = (0, 0, 0)
+    diameter = 0.5
+    magnetization_sphere = np.array((1e6, 2e6, 3e6))
+    moment_dipole = magnetization_sphere * diameter**3 / 6 * np.pi
+
+    sphere = magpy.magnet.Sphere(
+        diameter=diameter, magnetization=magnetization_sphere, position=pos
+    )
+    sphere.meshing = 100
+
+    dipole = magpy.misc.Dipole(position=pos, moment=moment_dipole)
+
+    FT_sphere = getFT(src, sphere, anchor=(0, 0, 0))
+    FT_dipole = getFT(src, dipole, anchor=(0, 0, 0))
+
+    np.testing.assert_allclose(FT_sphere, FT_dipole, rtol=1e-6, atol=1e-6)