Linear Filament

Fields

cfsem.flux_density_linear_filament

flux_density_linear_filament(
    xyzp: Array3xN,
    xyzfil: Array3xN,
    dlxyzfil: Array3xN,
    ifil: NDArray[float64],
    wire_radius: float | NDArray[float64] = 0.0,
    par: bool = True,
    output: Literal["vector", "matrix"] = "vector",
) -> tuple[
    NDArray[float64], NDArray[float64], NDArray[float64]
]

Biot-Savart law calculation for B-field contributions from many filament segments to many observation points.

Parameters:

Name	Type	Description	Default
xyzp	Array3xN	[m] x,y,z coords of observation points	required
xyzfil	Array3xN	[m] x,y,z coords of filament segment start points	required
dlxyzfil	Array3xN	[m] x,y,z deltas from segment start to segment end	required
ifil	NDArray[float64]	[A] current in each filament segment	required
wire_radius	float | NDArray[float64]	[m] filament radius, scalar or array of length m	0.0
par	bool	Whether to use CPU parallelism	True
output	Literal['vector', 'matrix']	"vector" for contracted field values at each target point, or "matrix" for row-major (nobs, nfil) source-target interaction matrices	'vector'

Returns:

Type	Description
NDArray[float64]	[T] (Bx, By, Bz) magnetic flux density at observation points,
NDArray[float64]	or explicit (nobs, nfil) interaction matrices if output="matrix"

Source code in cfsem/bindings.py

def flux_density_linear_filament(
    xyzp: Array3xN,
    xyzfil: Array3xN,
    dlxyzfil: Array3xN,
    ifil: NDArray[float64],
    wire_radius: float | NDArray[float64] = 0.0,
    par: bool = True,
    output: Literal["vector", "matrix"] = "vector",
) -> tuple[NDArray[float64], NDArray[float64], NDArray[float64]]:
    """
    Biot-Savart law calculation for B-field contributions from many filament segments
    to many observation points.

    Args:
        xyzp: [m] x,y,z coords of observation points
        xyzfil: [m] x,y,z coords of filament segment start points
        dlxyzfil: [m] x,y,z deltas from segment start to segment end
        ifil: [A] current in each filament segment
        wire_radius: [m] filament radius, scalar or array of length `m`
        par: Whether to use CPU parallelism
        output: `"vector"` for contracted field values at each target point,
            or `"matrix"` for row-major `(nobs, nfil)` source-target interaction matrices

    Returns:
        [T] (Bx, By, Bz) magnetic flux density at observation points,
        or explicit `(nobs, nfil)` interaction matrices if `output="matrix"`
    """
    xyzp = _3tup_contig(xyzp)
    xyzfil = _3tup_contig(xyzfil)
    dlxyzfil = _3tup_contig(dlxyzfil)
    ifil = ascontiguousarray(ifil).ravel()
    if asarray(wire_radius).ndim == 0:
        wire_radius = full(ifil.size, float(wire_radius))
    wire_radius = ascontiguousarray(wire_radius).ravel()
    if output == "vector":
        return em_flux_density_linear_filament(xyzp, xyzfil, dlxyzfil, ifil, wire_radius, par)
    if output == "matrix":
        bx, by, bz = em_flux_density_linear_filament_matrix(xyzp, xyzfil, dlxyzfil, ifil, wire_radius, par)
        nobs = xyzp[0].size
        nfil = ifil.size
        return (
            bx.reshape((nobs, nfil)),
            by.reshape((nobs, nfil)),
            bz.reshape((nobs, nfil)),
        )
    raise ValueError("output must be 'vector' or 'matrix'")

cfsem.flux_density_linear_filament_hierarchical builtin

flux_density_linear_filament_hierarchical(
    xyzp: ArrayTriple,
    xyzfil: ArrayTriple,
    dlxyzfil: ArrayTriple,
    ifil: FloatArray,
    wire_radius: FloatArray,
    theta: float = 0.05,
    construction_method: str = "longest_axis",
    par: bool = True,
    out: ArrayTriple | None = None,
    extra_diagnostics: bool = False,
) -> SolveResult

Evaluate linear-filament flux density with the hierarchical solver from Python inputs.

cfsem.vector_potential_linear_filament

vector_potential_linear_filament(
    xyzp: Array3xN,
    xyzfil: Array3xN,
    dlxyzfil: Array3xN,
    ifil: NDArray[float64],
    wire_radius: float | NDArray[float64] = 0.0,
    par: bool = True,
    output: Literal["vector", "matrix"] = "vector",
) -> tuple[
    NDArray[float64], NDArray[float64], NDArray[float64]
]

Vector potential calculation for A-field contribution from many current filament segments to many observation points.

Parameters:

Name	Type	Description	Default
xyzp	Array3xN	[m] x,y,z coords of observation points	required
xyzfil	Array3xN	[m] x,y,z coords of filament segment start points	required
dlxyzfil	Array3xN	[m] x,y,z deltas from segment start to segment end	required
ifil	NDArray[float64]	[A] current in each filament segment	required
wire_radius	float | NDArray[float64]	[m] filament radius, scalar or array of length m	0.0
par	bool	Whether to use CPU parallelism	True
output	Literal['vector', 'matrix']	"vector" for contracted field values at each target point, or "matrix" for row-major (nobs, nfil) source-target interaction matrices	'vector'

Returns:

Type	Description
NDArray[float64]	[Wb/m] or [V-s/m] (Ax, Ay, Az) magnetic vector potential at observation points,
NDArray[float64]	or explicit (nobs, nfil) interaction matrices if output="matrix"

Source code in cfsem/bindings.py

def vector_potential_linear_filament(
    xyzp: Array3xN,
    xyzfil: Array3xN,
    dlxyzfil: Array3xN,
    ifil: NDArray[float64],
    wire_radius: float | NDArray[float64] = 0.0,
    par: bool = True,
    output: Literal["vector", "matrix"] = "vector",
) -> tuple[NDArray[float64], NDArray[float64], NDArray[float64]]:
    """
    Vector potential calculation for A-field contribution from many current filament
    segments to many observation points.

    Args:
        xyzp: [m] x,y,z coords of observation points
        xyzfil: [m] x,y,z coords of filament segment start points
        dlxyzfil: [m] x,y,z deltas from segment start to segment end
        ifil: [A] current in each filament segment
        wire_radius: [m] filament radius, scalar or array of length `m`
        par: Whether to use CPU parallelism
        output: `"vector"` for contracted field values at each target point,
            or `"matrix"` for row-major `(nobs, nfil)` source-target interaction matrices

    Returns:
        [Wb/m] or [V-s/m] (Ax, Ay, Az) magnetic vector potential at observation points,
        or explicit `(nobs, nfil)` interaction matrices if `output="matrix"`
    """
    xyzp = _3tup_contig(xyzp)
    xyzfil = _3tup_contig(xyzfil)
    dlxyzfil = _3tup_contig(dlxyzfil)
    ifil = ascontiguousarray(ifil).ravel()
    if asarray(wire_radius).ndim == 0:
        wire_radius = full(ifil.size, float(wire_radius))
    wire_radius = ascontiguousarray(wire_radius).ravel()
    if output == "vector":
        return em_vector_potential_linear_filament(xyzp, xyzfil, dlxyzfil, ifil, wire_radius, par)
    if output == "matrix":
        ax, ay, az = em_vector_potential_linear_filament_matrix(
            xyzp, xyzfil, dlxyzfil, ifil, wire_radius, par
        )
        nobs = xyzp[0].size
        nfil = ifil.size
        return (
            ax.reshape((nobs, nfil)),
            ay.reshape((nobs, nfil)),
            az.reshape((nobs, nfil)),
        )
    raise ValueError("output must be 'vector' or 'matrix'")

cfsem.vector_potential_linear_filament_hierarchical builtin

vector_potential_linear_filament_hierarchical(
    xyzp: ArrayTriple,
    xyzfil: ArrayTriple,
    dlxyzfil: ArrayTriple,
    ifil: FloatArray,
    wire_radius: FloatArray,
    theta: float = 0.05,
    construction_method: str = "longest_axis",
    par: bool = True,
    out: ArrayTriple | None = None,
    extra_diagnostics: bool = False,
) -> SolveResult

Evaluate linear-filament vector potential with the hierarchical solver from Python inputs.

Force

cfsem.body_force_density_linear_filament

body_force_density_linear_filament(
    xyzfil: Array3xN,
    dlxyzfil: Array3xN,
    ifil: NDArray[float64],
    obs: Array3xN,
    j: Array3xN,
    wire_radius: float | NDArray[float64] = 0.0,
    par: bool = True,
) -> Array3xN

JxB (Lorentz) body force density (per volume) due to a linear current filament segment at an observation point with some current density (per area).

Parameters:

Name	Type	Description	Default
xyzfil	Array3xN	[m] x,y,z coords of current filament origins (start of segment)	required
dlxyzfil	Array3xN	[m] x,y,z length delta of current filaments	required
ifil	NDArray[float64]	[A] filament current	required
obs	Array3xN	[m] x,y,z coords of observation locations	required
j	Array3xN	[A/m^2] current density vector at observation locations	required
wire_radius	float | NDArray[float64]	[m] filament radius, scalar or array of length m	0.0
par	bool	Whether to use CPU parallelism	True

Returns:

Type	Description
Array3xN	[N/m^3] body force density

Source code in cfsem/bindings.py

def body_force_density_linear_filament(
    xyzfil: Array3xN,
    dlxyzfil: Array3xN,
    ifil: NDArray[float64],
    obs: Array3xN,
    j: Array3xN,
    wire_radius: float | NDArray[float64] = 0.0,
    par: bool = True,
) -> Array3xN:
    """
    JxB (Lorentz) body force density (per volume) due to a linear current
    filament segment at an observation point with some current density (per area).

    Args:
        xyzfil: [m] x,y,z coords of current filament origins (start of segment)
        dlxyzfil: [m] x,y,z length delta of current filaments
        ifil: [A] filament current
        obs: [m] x,y,z coords of observation locations
        j: [A/m^2] current density vector at observation locations
        wire_radius: [m] filament radius, scalar or array of length `m`
        par: Whether to use CPU parallelism

    Returns:
        [N/m^3] body force density
    """
    xyzfil = _3tup_contig(xyzfil)
    dlxyzfil = _3tup_contig(dlxyzfil)
    ifil = ascontiguousarray(ifil).ravel()
    obs = _3tup_contig(obs)
    j = _3tup_contig(j)
    if asarray(wire_radius).ndim == 0:
        wire_radius = full(ifil.size, float(wire_radius))
    wire_radius = ascontiguousarray(wire_radius).ravel()
    jxbx, jxby, jxbz = em_body_force_density_linear_filament(
        xyzfil, dlxyzfil, ifil, obs, j, wire_radius, par
    )  # [N/m^3]

    return jxbx, jxby, jxbz

Paths

cfsem.filament_coil

filament_coil(
    r: float,
    z: float,
    w: float,
    h: float,
    nt: float,
    nr: int,
    nz: int,
) -> NDArray

Create an array of filaments from coil cross-section, evenly spaced inside the winding pack. No filaments are coincident with the coil surface.

Parameters:

Name	Type	Description	Default
r	float	[m] radius, coil center	required
z	float	[m] axial position, coil center	required
w	float	[m] width of coil pack	required
h	float	[m] height of coil pack	required
nt	float	turns	required
nr	int	radial discretizations	required
nz	int	axial discretizations	required

Returns:

Type	Description
NDArray	(nr*nz) x 3, (r,z,n) of each filament

Source code in cfsem/__init__.py

def filament_coil(r: float, z: float, w: float, h: float, nt: float, nr: int, nz: int) -> NDArray:
    """
    Create an array of filaments from coil cross-section, evenly spaced
    _inside_ the winding pack. No filaments are coincident with the coil surface.

    Args:
        r: [m] radius, coil center
        z: [m] axial position, coil center
        w: [m] width of coil pack
        h: [m] height of coil pack
        nt: turns
        nr: radial discretizations
        nz: axial discretizations

    Returns:
        (nr*nz) x 3, (r,z,n) of each filament
    """

    # Build a 2D mesh of points evenly spaced in the interior of the bounding rectangle
    rs = np.linspace(r - w * (nr - 1) / nr / 2, r + w * (nr - 1) / nr / 2, nr)  # [m]
    zs = np.linspace(z - h * (nz - 1) / nz / 2, z + h * (nz - 1) / nz / 2, nz)  # [m]
    rmesh, zmesh = np.meshgrid(rs, zs, indexing="ij")  # [m]

    # Number of turns attributed to each point is not necessarily an integer
    n = np.full_like(rmesh.ravel(), float(nt) / (nr * nz))

    # Pack filament locations and number of turns into an array together
    filaments = np.dstack([rmesh.ravel(), zmesh.ravel(), n]).reshape(nr * nz, 3)

    return filaments  # [m], [m], [dimensionless]

cfsem.filament_helix_path

filament_helix_path(
    path: Array3xN,
    helix_start_offset: tuple[float, float, float],
    twist_pitch: float,
    angle_offset: float,
) -> Array3xN

Filamentize a helix about an arbitrary piecewise-linear path.

Assumes angle between sequential path segments is small and will fail if that angle approaches or exceeds 90 degrees.

The helix initial position vector, helix_start_offset, must be in a plane normal to the first path segment in order to produce good results. If it is not in-plane, it will be projected on to that plane and then scaled to the magnitude of its original length s.t. the distance from the helix to the path center is preserved but its orientation is not.

Description of the method:

1. Translate [filament segment n-1] to the base of [path segment n] and call it [filament segment n]
2. Take cross product of [path segment n] with [path segment n-1]
3. Rotate [filament segment n] segment about the axis of that cross product to bring it into the plane defined by [path segment n] as a normal vector
4. Rotate [filament seg. n] about [path seg. n] to continue the helix orbit

Parameters:

Name	Type	Description	Default
path	Array3xN	[m] 3xN Centerline points	required
helix_start_offset	tuple[float, float, float]	[m] (3x1) Initial position of helix rel. to centerline path	required
twist_pitch	float	[m] (scalar) Centerline length per helix orbit	required
angle_offset	float	[rad] (scalar) Initial rotation offset about centerline	required

Returns:

Type	Description
Array3xN	[m] 3xN array of points on the helix that twists around the path

Source code in cfsem/bindings.py

def filament_helix_path(
    path: Array3xN,
    helix_start_offset: tuple[float, float, float],
    twist_pitch: float,
    angle_offset: float,
) -> Array3xN:
    """
    Filamentize a helix about an arbitrary piecewise-linear path.

    Assumes angle between sequential path segments is small and will fail
    if that angle approaches or exceeds 90 degrees.

    The helix initial position vector, helix_start_offset, must be in a plane normal to
    the first path segment in order to produce good results. If it is not in-plane,
    it will be projected on to that plane and then scaled to the magnitude of its
    original length s.t. the distance from the helix to the path center is preserved
    but its orientation is not.

    Description of the method:

    1. Translate [filament segment n-1] to the base of [path segment n]
        and call it [filament segment n]
    2. Take cross product of [path segment n] with [path segment n-1]
    3. Rotate [filament segment n] segment about the axis of that cross product
        to bring it into the plane defined by [path segment n] as a normal vector
    4. Rotate [filament seg. n] about [path seg. n] to continue the helix orbit

    Args:
        path: [m] 3xN Centerline points
        helix_start_offset: [m] (3x1) Initial position of helix rel. to centerline path
        twist_pitch: [m] (scalar) Centerline length per helix orbit
        angle_offset: [rad] (scalar) Initial rotation offset about centerline

    Returns:
        [m] 3xN array of points on the helix that twists around the path
    """

    # Make sure input is contiguous, reallocating only if necessary
    path = ascontiguousarray(path)

    # Allocate output
    helix = zeros_like(path)  # [m]

    # Calculate, mutating output
    em_filament_helix_path(
        (*path,),
        helix_start_offset,
        twist_pitch,
        angle_offset,
        (*helix,),
    )

    return helix  # [m]

cfsem.rotate_filaments_about_path

rotate_filaments_about_path(
    path: Array3xN, angle_offset: float, fils: Array3xN
) -> Array3xN

Rotate a path of point about another path.

Intended for rotating a helix generated by filament_helix_path about the centerline that was used to generate it.

Parameters:

Name	Type	Description	Default
path	Array3xN	[m] x,y,z Centerline points	required
angle_offset	float	[rad] (scalar) Initial rotation offset about centerline	required
fils	Array3xN	[m] x,y,z Filaments to rotate around centerline	required

Returns:

Type	Description
Array3xN	[m] 3xN array of points on the helix that twists around the path

Source code in cfsem/bindings.py

def rotate_filaments_about_path(path: Array3xN, angle_offset: float, fils: Array3xN) -> Array3xN:
    """
    Rotate a path of point about another path.

    Intended for rotating a helix generated by [`filament_helix_path`][cfsem.filament_helix_path]
    about the centerline that was used to generate it.

    Args:
        path: [m] x,y,z Centerline points
        angle_offset: [rad] (scalar) Initial rotation offset about centerline
        fils: [m] x,y,z Filaments to rotate around centerline

    Returns:
        [m] 3xN array of points on the helix that twists around the path
    """

    # Make sure input is contiguous, reallocating only if necessary
    path = ascontiguousarray(path)

    new_fils = ascontiguousarray(fils).copy()

    em_rotate_filaments_about_path(
        (*path,),
        angle_offset,
        (*new_fils,),
    )

    return new_fils  # [m]