vector Line code, replaces also Line classic

Author: OrtnerMichael

magpylib/_lib/classes/currents.py

@@ -24,7 +24,7 @@
 
 from numpy import array, float64, ndarray
 from magpylib._lib.mathLib import angleAxisRotation_priv
-from magpylib._lib.fields.Current_Line import Bfield_CurrentLine
+from magpylib._lib.fields.Current_Line_vector import Bfield_CurrentLineV
 from magpylib._lib.fields.Current_CircularLoop import Bfield_CircularCurrentLoop
 from magpylib._lib.classes.base import LineCurrent
 
@@ -238,7 +238,7 @@ def getB(self, pos):  # Particular Line current B field calculation. Check RCS f
         # rotate this vector into the CS of the magnet (inverse rotation)
         rotatedPos = angleAxisRotation_priv(self.angle, -self.axis, posRel) # pylint: disable=invalid-unary-operand-type
         # rotate field vector back
-        BCm = angleAxisRotation_priv(self.angle, self.axis, Bfield_CurrentLine(rotatedPos, self.vertices, self.current))
+        BCm = angleAxisRotation_priv(self.angle, self.axis, Bfield_CurrentLineV(self.vertices, self.current,rotatedPos))
         # BCm is the obtained magnetic field in Cm
         # the field is well known in the magnet coordinates.
         return BCm

magpylib/_lib/fields/Current_Line.py renamed to magpylib/_lib/fields/Current_Line_vector.py

@@ -25,6 +25,8 @@
 from numpy import array, NaN
 from magpylib._lib.mathLib import fastSum3D, fastNorm3D, fastCross3D
 from warnings import warn
+import numpy as np
+from numpy.linalg import norm
 
 # %% CURRENT LINE
 # Describes the magnetic field of a line current. The line is given by a set of
@@ -38,6 +40,96 @@
 # source: http://www.phys.uri.edu/gerhard/PHY204/tsl216.pdf
 # FieldOfLineCurrent
 
+# VECTORIZED VERSION
+# here we only use vectorized code as this function will primarily be
+#   used to call on multiple line segments. The vectorized code was
+#   developed based on the depreciated version below.
+
+def Bfield_LineSegmentV(p0, p1, p2, I0):
+    ''' private
+    base function determines the fields of given line segments
+    p0 = observer position
+    p1->p2 = current flows from vertex p1 to vertex p2
+    I0 = current in [A]
+    '''
+
+    N = len(p0)
+    fields = np.zeros((N,3)) # default values for mask0 and mask1
+
+    # Check for zero-length segment
+    mask0 = np.all(p1==p2,axis=1)
+
+    # projection of p0 onto line p1-p2
+    nm0 = np.invert(mask0)
+    p1p2 = (p1[nm0]-p2[nm0])
+    p4 = p1[nm0]+(p1p2.T*np.sum((p0[nm0]-p1[nm0])*p1p2,axis=1)/np.sum(p1p2**2,axis=1)).T
+
+    # determine anchorrect normal vector to surface spanned by triangle
+    cross0 = np.cross(-p1p2, p0[nm0]-p4)
+    norm_cross0 = norm(cross0,axis=1)
+
+    # on-line cases (include when position is on current path)
+    mask1 = (norm_cross0 == 0)
+
+    # normal cases
+    nm1 = np.invert(mask1)
+    eB = (cross0[nm1].T/norm_cross0[nm1]) #field direction
+
+    # not mask0 and not mask1
+    NM = np.copy(nm0)
+    NM[NM==True] = nm1
+
+    norm_04 = norm(p0[NM] -p4[nm1],axis=1)
+    norm_01 = norm(p0[NM] -p1[NM],axis=1)
+    norm_02 = norm(p0[NM] -p2[NM],axis=1)
+    norm_12 = norm(p1[NM] -p2[NM],axis=1)
+    norm_41 = norm(p4[nm1]-p1[NM],axis=1)
+    norm_42 = norm(p4[nm1]-p2[NM],axis=1)
+
+    sinTh1 = norm_41/norm_01
+    sinTh2 = norm_42/norm_02
+
+    deltaSin = np.empty((N))[NM]
+
+    # determine how p1,p2,p4 are sorted on the line (to get sinTH signs)
+    # both points below
+    mask2 = ((norm_41>norm_12) * (norm_41>norm_42))
+    deltaSin[mask2] = abs(sinTh1[mask2]-sinTh2[mask2])
+    # both points above
+    mask3 = ((norm_42>norm_12) * (norm_42>norm_41))
+    deltaSin[mask3] = abs(sinTh2[mask3]-sinTh1[mask3])
+    # one above one below or one equals p4
+    mask4 = np.invert(mask2)*np.invert(mask3)
+    deltaSin[mask4] = abs(sinTh1[mask4]+sinTh2[mask4])
+
+    # missing 10**-6 from m->mm conversion #T->mT conversion
+    fields[NM] = (I0[NM]*deltaSin/norm_04*eB).T/10
+
+    return fields
+
+
+
+def Bfield_CurrentLineV(VERT,i0,poso):
+    ''' private
+    determine total field from a multi-segment line current
+    '''
+
+    N = len(VERT)-1
+    P0 = np.outer(np.ones((N)),poso)
+    P1 = VERT[:-1]
+    P2 = VERT[1:]
+    I0 = np.ones((N))*i0
+
+    Bv = Bfield_LineSegmentV(P0,P1,P2,I0)
+
+    return np.sum(Bv,axis=0)
+
+
+
+
+
+''' DEPRECHIATED VERSION (non-vectorized)
+
 # observer at p0
 # current I0 flows in straight line from p1 to p2
 def Bfield_LineSegment(p0, p1, p2, I0):
@@ -90,15 +182,4 @@ def Bfield_LineSegment(p0, p1, p2, I0):
     B = I0*deltaSin/norm_04*eB/10
 
     return B
-
-# determine total field from multiple segments
-
-
-def Bfield_CurrentLine(p0, possis, I0):
-
-    B = array([0., 0., 0.])
-    for i in range(len(possis)-1):
-        p1 = possis[i]
-        p2 = possis[i+1]
-        B += Bfield_LineSegment(p0, p1, p2, I0)
-    return B
+'''

magpylib/_lib/getBvector.py

@@ -26,6 +26,7 @@
 from magpylib._lib.fields.PM_Sphere_vector import Bfield_SphereV
 from magpylib._lib.fields.Moment_Dipole_vector import Bfield_DipoleV
 from magpylib._lib.fields.Current_CircularLoop_vector import Bfield_CircularCurrentLoopV
+from magpylib._lib.fields.Current_Line_vector import Bfield_LineSegmentV
 from magpylib._lib.mathLib_vector import QconjV, QrotationV, QmultV, getRotQuatV, angleAxisRotationV_priv
 import numpy as np
 
@@ -121,7 +122,7 @@ def getBv_current(type,CURR,DIM,POSm,POSo,ANG=[],AX=[],ANCH=[]):
 
     DIM : NxY numpy array float [mm]
         vector of N dimensions for each evaluation. The form of this vector depends
-        on the source type. Y=3/2/1 for box/cylinder/sphere
+        on the source type. Y=1/3x3 for circular/line.
 
     POSo : Nx3 numpy array float [mm]
         vector of N positions of the observer.
@@ -162,6 +163,8 @@ def getBv_current(type,CURR,DIM,POSm,POSo,ANG=[],AX=[],ANCH=[]):
     # calculate field
     if type == 'circular':
         Brot = Bfield_CircularCurrentLoopV(CURR, DIM, POSrot)
+    elif type == 'line':
+        Brot = Bfield_LineSegmentV(POSrot,DIM[:,0],DIM[:,1],CURR)
     else:
         print('Bad type')
         return 0

tests/test_magpyVector.py

@@ -3,7 +3,7 @@
 import numpy as np
 from magpylib.source.magnet import Box, Cylinder, Sphere
 from magpylib.source.moment import Dipole
-from magpylib.source.current import Circular
+from magpylib.source.current import Circular, Line
 from magpylib.vector import getBv_magnet, getBv_current, getBv_moment
 from magpylib.math import axisFromAngles
 from magpylib.math import angleAxisRotationV
@@ -156,3 +156,37 @@ def test_vectorCurrentCircular():
     Bv = getBv_current('circular',I,D,Pm,Po)
 
     assert np.amax(abs(Bc-Bv))<1e-10
+
+
+def test_vectorLine():
+
+    #general cases
+    NN=100
+    V = np.random.rand(NN,3)-.5
+    V1 = V[:-1]
+    V2 = V[1:]
+    DIM = np.array([[v1,v2] for v1,v2 in zip(V1,V2)])
+    I0 = np.ones([NN-1])
+    Po = np.ones((NN-1,3))
+    Pm = np.zeros((NN-1,3))
+    
+    Bc = []
+    for dim,i0,po in zip(DIM,I0,Po):
+        s = Line(curr=i0,vertices=dim)
+        Bc += [s.getB(po)]
+    Bc = np.array(Bc)
+
+    Bv = getBv_current('line',I0,DIM,Pm,Po)
+    assert np.amax(abs(Bc-Bv))<1e-12
+
+    #special cases
+    V = np.array([[0,0,0],[1,2,3],[1,2,3],[3,3,3],[1,1,1],[1,1,2],[1,1,0]])
+    V1 = V[:-1]
+    V2 = V[1:]
+    DIM = np.array([[v1,v2] for v1,v2 in zip(V1,V2)])
+    I0 = np.ones([6])
+    Po = np.ones((6,3))
+    Pm = np.zeros((6,3))
+
+    Bv = getBv_current('line',I0,DIM,Pm,Po)
+    assert [all(b == 0) for b in Bv] == [False, True, False, True, True, True]