Add current loop (#10)

* add Circle target

* changelog

* pylint

* version bump

Author: OrtnerMichael

CHANGELOG.md

@@ -1,5 +1,10 @@
 # Unpublished
 
+# v0.3.0 (2024/10/29)
+- Added the current.Cicle class to the force computation
+- Improved warning msg when meshing parameter not given
+- Added more physics tests
+
 # v0.2.0 (2024/10/11)
 - Added a warning when no anchor is given to getFT to avoid misuse.
 - Added Cylinder computation

magpylib_force/__init__.py

@@ -7,8 +7,8 @@
 from magpylib_force.force import getFT
 
 # module level dunders
-__version__ = "0.2.0"
-__author__ = "SAL"
+__version__ = "0.3.0"
+__author__ = "The Magpylib Team"
 __all__ = [
     "getFT",
 ]

magpylib_force/force.py

@@ -1,6 +1,7 @@
 """
 Force computation codes
 """
+import warnings
 
 import magpylib as magpy
 import numpy as np
@@ -10,6 +11,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_current_Circle import Circle
 
 from magpylib_force.meshing import mesh_target
 from magpylib_force.utility import check_input_anchor
@@ -63,8 +65,12 @@ 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]
-    getFT_FUNCS = [getFTmagnet, getFTcurrent, getFTmagnet, getFTmagnet, getFTmagnet]
+    TARGET_TYPES = [
+        Cuboid, Polyline, Sphere, Cylinder, CylinderSegment, Circle
+    ]
+    getFT_FUNCS = [
+        getFTmagnet, getFTcurrent, getFTmagnet, getFTmagnet, getFTmagnet, getFTcurrent_circ
+    ]
     objects = [[] for _ in TARGET_TYPES]
     orders  = [[] for _ in TARGET_TYPES]
 
@@ -192,6 +198,41 @@ def getFTmagnet(sources, targets, eps=1e-5, anchor=None):
 
     return np.array((F, -T))
 
+
+def getFTcurrent_circ(sources, targets, anchor=None, eps=None):
+    """
+    The current force computation is designed for Polylines.
+    To make use of this, Circle objects are transformed into Polylines.
+    Its a dirty solution that should be fixed at some point
+    """
+    new_targets = []
+    for tgt in targets:
+        n = tgt.meshing
+        if n<20:
+            warnings.warn(
+                "Circle meshing parameter with low value detected. "
+                "This will give bad results. Please increase meshing. "
+                "Circle meshing defines the number of points on the circle."
+                )
+        r = tgt.diameter/2
+        verts = np.zeros((3,n))
+        verts[2] = np.linspace(0, 2*np.pi, n)
+        verts[0] = r*np.cos(verts[2])
+        verts[1] = r*np.sin(verts[2])
+        verts[2] = 0
+
+        poly = magpy.current.Polyline(
+            vertices=verts.T,
+            current=tgt.current,
+            position=tgt.position,
+            orientation=tgt.orientation,
+        )
+        poly.meshing=1
+        new_targets.append(poly)
+
+    return getFTcurrent(sources, new_targets, anchor, eps)
+
+
 #pylint: disable=unused-argument
 def getFTcurrent(sources, targets, anchor=None, eps=None):
     """
@@ -203,7 +244,6 @@ def getFTcurrent(sources, targets, anchor=None, eps=None):
     segements = linear segments within Polyline objects
     instances = computation instances, each segment is split into `meshing` points
     """
-
     # number of Polylines
     tgt_number = len(targets)
 

magpylib_force/meshing.py

@@ -13,7 +13,6 @@
 
 #pylint: disable=too-many-locals
 
-
 def mesh_target(obj):
     """
     create mesh for target objects

magpylib_force/utility.py

@@ -6,6 +6,7 @@
 import numpy as np
 from magpylib._src.obj_classes.class_magnet_Cuboid import Cuboid
 from magpylib._src.obj_classes.class_current_Polyline import Polyline
+from magpylib._src.obj_classes.class_current_Circle import Circle
 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
@@ -34,16 +35,16 @@ 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)):
+        if not isinstance(t, (Cuboid, Polyline, Sphere, Cylinder, CylinderSegment, Circle)):
             raise ValueError(
                 "Bad `targets` input for getFT."
-                " `targets` can only be Cuboids, Polylines, Spheres, Cylinders, and"
-                " CylinderSegments."
+                " `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(
-                "Bad `targets` input for getFT."
+                "Missing input for getFT `targets`."
                 " `targets` must have the `meshing` parameter set."
             )
     return targets

tests/test_physics.py

@@ -2,6 +2,8 @@
 import numpy as np
 import magpylib as magpy
 from magpylib_force.force import getFT
+from scipy.special import ellipe
+from scipy.special import ellipk
 
 
 # def test_physics_loop_torque():
@@ -304,3 +306,29 @@ def test_torque_sign():
     _,T = getFT(mag, loop, anchor=(0,0,0))
 
     assert T[1] < 0
+
+
+def test_force_between_cocentric_loops():
+    """
+    compare the numerical solution against the analytical solution of the force between two
+    cocentric current loops.
+    See e.g. IEEE TRANSACTIONS ON MAGNETICS, VOL. 49, NO. 8, AUGUST 2013
+    """
+    z1, z2 = 0.123, 1.321
+    i1, i2 = 3.2, 5.1
+    r1, r2 = 1.2, 2.3
+
+    # numerical solution
+    loop1 = magpy.current.Circle(diameter=2*r1, current=i1, position=(0,0,z1))
+    loop2 = magpy.current.Circle(diameter=2*r2, current=i2, position=(0,0,z2))
+    loop2.meshing=1000
+    F_num = getFT(loop1, loop2, anchor=(0,0,0))[0,2]
+    
+    # analytical solution
+    k2 = 4*r1*r2 / ((r1+r2)**2+(z1-z2)**2)
+    k = np.sqrt(k2)
+    pf = magpy.mu_0*i1*i2*(z1-z2)*k / 4 / np.sqrt(r1*r2)
+    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
+        

tests/test_self_consistency.py

@@ -229,3 +229,41 @@ def test_consistency_cylinder_segment_cuboid():
         ft2 = getFT(cube, cyls, anchor=(0,0,0))
 
         assert np.amax(abs((ft1+ft2)/(ft1-ft2))) < 0.09
+
+
+def test_consistency_polyline_circle():
+    """
+    compare Polyline solution to circle solution
+    """
+
+    src = magpy.magnet.Sphere(diameter=1, polarization=(1,2,3), position=(0,0,-1))
+    
+    # circle
+    loop1 = magpy.current.Circle(diameter=3, current=123)
+    loop1.meshing=500
+    # polyline
+    rr = loop1.diameter/2
+    ii = loop1.current
+    phis = np.linspace(0,2*np.pi,500)
+    verts = [(rr*np.cos(p), rr*np.sin(p), 0) for p in phis]
+    loop2 = magpy.current.Polyline(current=ii, vertices=verts)
+    loop2.meshing=1
+    
+    F1,T1 = getFT(src, loop1, anchor=(0,0,0))
+    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
+
+    loop1.move((1.123,2.321,.123))
+    loop2.move((1.123,2.321,.123))
+    F1,T1 = getFT(src, loop1, anchor=(0,0,0))
+    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
+
+    loop1.rotate_from_angax(20, 'x')
+    loop2.rotate_from_angax(20, 'x')
+    F1,T1 = getFT(src, loop1, anchor=(0,0,0))
+    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