force computation for dipole

Author: fslanovc

magpylib_force/force.py

@@ -12,6 +12,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_current_Circle import Circle
+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
@@ -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
@@ -66,10 +67,10 @@ def getFT(sources, targets, anchor=None, eps=1e-5, squeeze=True):
 
     # split targets into lists of similar types
     TARGET_TYPES = [
-        Cuboid, Polyline, Sphere, Cylinder, CylinderSegment, Circle
+        Cuboid, Polyline, Sphere, Cylinder, CylinderSegment, Circle, Dipole
     ]
     getFT_FUNCS = [
-        getFTmagnet, getFTcurrent, getFTmagnet, getFTmagnet, getFTmagnet, getFTcurrent_circ
+        getFTmagnet, getFTcurrent, getFTmagnet, getFTmagnet, getFTmagnet, getFTcurrent_circ, getFTdipole
     ]
     objects = [[] for _ in TARGET_TYPES]
     orders  = [[] for _ in TARGET_TYPES]
@@ -184,6 +185,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)
 
@@ -300,3 +303,61 @@ def getFTcurrent(sources, targets, anchor=None, eps=None):
     F = np.array([np.sum(F[insti[i]:insti[i+1]],axis=0) for i in range(tgt_number)])
 
     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 ommitted.
+    """
+    # 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))

magpylib_force/utility.py

@@ -10,6 +10,7 @@
 from magpylib._src.obj_classes.class_magnet_Sphere import Sphere
 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_misc_Dipole import Dipole
 
 def check_input_anchor(anchor):
     """
@@ -35,23 +36,24 @@ def check_input_targets(targets):
     if not isinstance(targets, list):
         targets = [targets]
     for t in targets:
-        if not isinstance(t, (Cuboid, Polyline, Sphere, Cylinder, CylinderSegment, Circle)):
+        if not isinstance(t, (Cuboid, Polyline, Sphere, Cylinder, CylinderSegment, Circle, Dipole)):
             raise ValueError(
                 "Bad `targets` input for getFT."
                 " `targets` can only be Cuboids, Polylines, Spheres, Cylinders, "
                 " CylinderSegments, and Circles."
                 f" Instead receivd type {type(t)} target."
             )
-        if not hasattr(t, "meshing"):
-            raise ValueError(
-                "Missing input for getFT `targets`."
-                " `targets` must have the `meshing` parameter set."
-            )
-        if not isinstance(t, (Polyline, )):
-            if np.isscalar(t.meshing):
-                if t.meshing<20:
-                    warnings.warn(
-                        "Input parameter `meshing` is set to a low value which will result in "
-                        "inaccurate computation of force and torque."
-                    )
+        if not isinstance(t, (Dipole, )):
+            if not hasattr(t, "meshing"):
+                raise ValueError(
+                    "Missing input for getFT `targets`."
+                    " `targets` must have the `meshing` parameter set."
+                )
+            if not isinstance(t, (Polyline, )):
+                if np.isscalar(t.meshing):
+                    if t.meshing<20:
+                        warnings.warn(
+                            "Input parameter `meshing` is set to a low value which will result in "
+                            "inaccurate computation of force and torque."
+                        )
     return targets

tests/test_physics.py

@@ -239,3 +239,31 @@ 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)