Python API

oersted

Python bindings for oersted

BHCurve dataclass

B-H curve for a nonlinear magnetic material

Source code in src/oersted/materials.py

@dataclass
class BHCurve:
    """B-H curve for a nonlinear magnetic material"""

    h_values: NDArray[float64]
    b_values: NDArray[float64]

    def lookup(self, h: float) -> float:
        """Linearly interpolate from the B-H curve"""

        return float(interp(h, self.h_values, self.b_values))

lookup

lookup(h: float) -> float

Linearly interpolate from the B-H curve

Source code in src/oersted/materials.py

def lookup(self, h: float) -> float:
    """Linearly interpolate from the B-H curve"""

    return float(interp(h, self.h_values, self.b_values))

CentroidMesh

A finite element mesh represented solely by the centroidal values of the elements

This is used in the point source calculations. It is an approximation, but extremely fast and accurate for far field or force calculations.

Source code in src/oersted/mesh.py

class CentroidMesh:
    """A finite element mesh represented solely by the centroidal values of the elements

    This is used in the `point source` calculations. It is an approximation, but extremely
    fast and accurate for far field or force calculations.
    """

    # Topology data
    _centroids: NDArray[float64]
    _volumes: NDArray[float64]

    def __init__(self, centroids: NDArray[float64], volumes: NDArray[float64]):
        self._centroids = centroids
        self._volumes = volumes

    @property
    def num_elems(self) -> int:
        return self._centroids.shape[0]

    @property
    def centroids(self):
        return self._centroids

    @property
    def volumes(self):
        return self._volumes

LinearMaterial dataclass

Bases: Material

A linear magnetic material that has a constant mu_r for all values of applied H-field.

Source code in src/oersted/materials.py

@dataclass
class LinearMaterial(Material):
    """A linear magnetic material that has a constant mu_r for all
    values of applied H-field.
    """

    _mu_r: float

    def __init__(self, _mu_r):
        self._mu_r = _mu_r

    def mu_r(self, h: float):
        return self._mu_r

    def to_mh_curve(self) -> tuple[NDArray[float64], NDArray[float64]]:
        h_values = array([0.0, 1e10])
        m_values = array([(self._mu_r - 1.0) * h for h in h_values])
        return (h_values, m_values)

Mesh

A continuous finite element mesh made of tet4 elements

Source code in src/oersted/mesh.py

class Mesh:
    """A continuous finite element mesh made of tet4 elements"""

    # Basic mesh topology data
    _nodes: NDArray[float64]
    _connectivity: NDArray[uint32]

    # Information that is compute on-demand for the user
    _edges: NDArray[uint32] | None
    _faces: NDArray[uint32] | None
    _centroids: NDArray[float64] | None
    _volumes: NDArray[float64] | None

    # Surface mesh data
    _surface: SurfaceMesh | None

    def __init__(self, nodes: NDArray[float64], connectivity: NDArray[uint32]):
        assert len(nodes.shape) == 2
        assert nodes.shape[1] == 3
        assert len(connectivity.shape) == 2
        assert connectivity.shape[1] == 4
        self._nodes = nodes
        self._connectivity = connectivity

        self._edges = None
        self._faces = None
        self._centroids = None
        self._volumes = None
        self._surface = None

    @property
    def nodes(self) -> NDArray[float64]:
        """Returns an (N,3) array of nodal coordinates in the mesh"""
        return self._nodes

    @property
    def connectivity(self) -> NDArray[uint32]:
        """Returns an (N,4) array of the node numbers associated with each element

        Node numbers are indices into the self._nodes array
        """
        return self._connectivity

    @property
    def num_nodes(self) -> int:
        """Returns the number of nodes in the model"""
        return self._nodes.shape[0]

    @property
    def num_elems(self) -> int:
        """Returns the number of elements in the model"""
        return self._connectivity.shape[0]

    def to_centroid_mesh(self) -> CentroidMesh:
        """Create a centroid mesh from a tet4 mesh"""
        return CentroidMesh(self.centroids, self.volumes)

    @property
    def edges(self):  # -> NDArray[uint32]:
        """Returns an (N,2) array of edges in the model

        Each value in the array is a node number associated with that edge.
        The first node is the start node, the second is the end node. This
        provides directionality for the edge.
        """
        if self._edges is None:
            # self._edges = _mesh_edges(self.nodes, self.connectivity)
            pass
        return self._edges

    @property
    def faces(self):
        """Returns an (N,3) array of nodes associated with each element face
        in the model

        Nodes are ordered such that the right hand rule forms the face normal.
        """
        if self._faces is None:
            # self._faces = mesh_faces(self.nodes, self.connectivity)
            pass
        return self._faces

    @property
    def centroids(self) -> NDArray[float64]:
        """Returns an (N,3) array of all element centroids in the mesh"""
        if self._centroids is None:
            self._centroids = mesh_centroids(
                np.ascontiguousarray(self.nodes),
                np.ascontiguousarray(self.connectivity),
            )

        return self._centroids

    @property
    def volumes(self) -> NDArray[float64]:
        """Return an (N,) array of the volume of each element in the mesh"""
        if self._volumes is None:
            self._volumes = mesh_volumes(
                np.ascontiguousarray(self.nodes),
                np.ascontiguousarray(self.connectivity),
            )
        return self._volumes

    @property
    def surface(self):
        if self._surface is None:
            faces = mesh_surface_faces(self.connectivity)
            self._surface = SurfaceMesh(self.nodes.copy(), faces)

        return self._surface

centroids property

centroids: NDArray[float64]

Returns an (N,3) array of all element centroids in the mesh

connectivity property

connectivity: NDArray[uint32]

Returns an (N,4) array of the node numbers associated with each element

Node numbers are indices into the self._nodes array

edges property

edges

Returns an (N,2) array of edges in the model

Each value in the array is a node number associated with that edge. The first node is the start node, the second is the end node. This provides directionality for the edge.

faces property

faces

Returns an (N,3) array of nodes associated with each element face in the model

Nodes are ordered such that the right hand rule forms the face normal.

nodes property

nodes: NDArray[float64]

Returns an (N,3) array of nodal coordinates in the mesh

num_elems property

num_elems: int

Returns the number of elements in the model

num_nodes property

num_nodes: int

Returns the number of nodes in the model

volumes property

volumes: NDArray[float64]

Return an (N,) array of the volume of each element in the mesh

to_centroid_mesh

to_centroid_mesh() -> CentroidMesh

Create a centroid mesh from a tet4 mesh

Source code in src/oersted/mesh.py

def to_centroid_mesh(self) -> CentroidMesh:
    """Create a centroid mesh from a tet4 mesh"""
    return CentroidMesh(self.centroids, self.volumes)

NonlinearMaterial dataclass

Bases: Material

A nonlinear magnetic material that has a B-H curve

Source code in src/oersted/materials.py

@dataclass
class NonlinearMaterial(Material):
    """A nonlinear magnetic material that has a B-H curve"""

    curve: BHCurve

    def __init__(self, curve: BHCurve):
        self.curve = curve

    def mu_r(self, h: float) -> float:
        b: float = self.curve.lookup(h)

        if h != 0.0:
            return (b / oersted.MU0) / h - 1.0
        else:
            return 1.0

    def to_mh_curve(self) -> tuple[NDArray[float64], NDArray[float64]]:
        return (self.curve.h_values, self.curve.b_values / MU0 - self.curve.h_values)

b_field

b_field(
    source: Mesh | CentroidMesh,
    j_density: NDArray[float64],
    targets: NDArray[float64],
    solver: DirectSolver | OctreeSolver | None = None,
) -> NDArray[float64]

Compute the magnetic flux density at a collection of target points using the specific source mesh and solver options, assuming the target points are in free space

Source code in src/oersted/biotsavart.py

def b_field(
    source: Mesh | CentroidMesh, j_density: NDArray[float64], targets: NDArray[float64], solver: DirectSolver | OctreeSolver | None = None
) -> NDArray[float64]:
    """Compute the magnetic flux density at a collection of target points using the specific source mesh and
    solver options, assuming the target points are in free space
    """
    return MU0 * h_field(source, j_density, targets, solver)

h_field

h_field(
    source: Mesh | CentroidMesh,
    j_density: NDArray[float64],
    targets: NDArray[float64],
    solver: DirectSolver | OctreeSolver | None = None,
) -> NDArray[float64]

Compute the magnetic field strength at a collection of target points using a current-carrying source mesh.

Source code in src/oersted/biotsavart.py

def h_field(
    source: Mesh | CentroidMesh, j_density: NDArray[float64], targets: NDArray[float64], solver: DirectSolver | OctreeSolver | None = None
) -> NDArray[float64]:
    """Compute the magnetic field strength at a collection of target points using a current-carrying source mesh."""

    if solver is None:
        solver = DirectSolver()

    j_density: NDArray[float64] = ascontiguousarray(j_density, dtype=float64)
    tgt_pts: NDArray[float64] = ascontiguousarray(targets, dtype=float64)

    if isinstance(source, CentroidMesh):
        src_pts = ascontiguousarray(source.centroids, dtype=float64)
        src_vol = ascontiguousarray(source.volumes, dtype=float64)

        if isinstance(solver, DirectSolver):
            return (1.0 / MU0) * b_current_point_direct(src_pts, src_vol, j_density, tgt_pts, solver.n_threads)

        elif isinstance(solver, OctreeSolver):
            return h_current_point_octree(src_pts, src_vol, j_density, tgt_pts, solver.theta, solver.leaf_threshold, solver.n_threads)

        else:
            raise TypeError(f"Unsupported source/solver combination: {type(source)}, {type(solver)}")

    elif isinstance(source, Mesh):
        src_nodes = ascontiguousarray(source.nodes, dtype=float64)
        src_connectivity = ascontiguousarray(source.connectivity, dtype=uint32)
        if isinstance(solver, DirectSolver):
            return h_current_tet4_direct(src_nodes, src_connectivity, j_density, tgt_pts, solver.n_threads)

        elif isinstance(solver, OctreeSolver):
            return h_current_tet4_octree(src_nodes, src_connectivity, j_density, tgt_pts, solver.theta, solver.leaf_threshold, solver.n_threads)

        else:
            raise TypeError(f"Unsupported source/solver combination: {type(source)}, {type(solver)}")

    else:
        raise TypeError(f"Unsupported source/solver combination: {type(source)}, {type(solver)}")

h_mag

h_mag(
    source: Mesh | CentroidMesh,
    m_field: NDArray[float64],
    targets: NDArray[float64],
    solver: DirectSolver | OctreeSolver | None = None,
) -> NDArray[float64]

Compute the magnetic field strength using a magnetized mesh as the source

Source code in src/oersted/biotsavart.py

def h_mag(
    source: Mesh | CentroidMesh, m_field: NDArray[float64], targets: NDArray[float64], solver: DirectSolver | OctreeSolver | None = None
) -> NDArray[float64]:
    """Compute the magnetic field strength using a magnetized mesh as the source"""

    if solver is None:
        solver = DirectSolver()

    m_field: NDArray[float64] = ascontiguousarray(m_field, dtype=float64)
    targets: NDArray[float64] = ascontiguousarray(targets, dtype=float64)

    if isinstance(source, CentroidMesh):
        raise NotImplementedError("Centroid mesh not yet implemented")

    elif isinstance(source, Mesh):
        src_nodes = ascontiguousarray(source.nodes, dtype=float64)
        src_connectivity = ascontiguousarray(source.connectivity, dtype=uint32)
        if isinstance(solver, DirectSolver):
            theta = 0.0
            leaf_threshold: uint32 = uint32(0)
            return h_mag_tet4(src_nodes, src_connectivity, m_field, targets, theta, leaf_threshold, solver.n_threads)

        elif isinstance(solver, OctreeSolver):
            return h_mag_tet4(src_nodes, src_connectivity, m_field, targets, solver.theta, solver.leaf_threshold, solver.n_threads)

        else:
            raise TypeError(f"Unsupported source/solver combination: {type(source)}, {type(solver)}")

    else:
        raise TypeError(f"Unsupported source/solver combination: {type(source)}, {type(solver)}")

mesh_step

mesh_step(
    infile: str,
    min_size: float,
    max_size: float,
    scale=0.001,
) -> Mesh

Mesh a step file using gmsh

Source code in src/oersted/mesh.py

def mesh_step(infile: str, min_size: float, max_size: float, scale=1e-3) -> Mesh:
    """Mesh a step file using gmsh"""

    mshfile = infile.split(".")[0] + ".msh"
    nodes: NDArray[float64]
    connectivity: NDArray[uint32]

    try:
        import gmsh

        gmsh.initialize()
        gmsh.option.setNumber("General.Terminal", 0)  # suppress output
        gmsh.model.occ.importShapes(infile)
        gmsh.model.occ.synchronize()
        gmsh.option.setNumber("Mesh.CharacteristicLengthMin", min_size)
        gmsh.option.setNumber("Mesh.CharacteristicLengthMax", max_size)
        gmsh.model.mesh.generate(3)  # mesh 3d elements
        gmsh.write(mshfile)

        print(f"Wrote gmsh mesh to `{mshfile}")

        # Get all nodes: returns (tags, coords, parametricCoords)
        node_tags, coords, _ = gmsh.model.mesh.getNodes()

        # coords is flat [x0,y0,z0,x1,y1,z1,...], reshape to (Nn, 3)
        nodes = np.array(coords).reshape(-1, 3) * scale

        # Build compact renumbering: gmsh tags can be sparse/non-sequential
        tag_to_compact = {tag: i for i, tag in enumerate(node_tags)}

        # Get tet elements (type 4 = 4-node tetrahedra)
        elem_tags, elem_node_tags = gmsh.model.mesh.getElementsByType(4)

        # elem_node_tags is flat [n0,n1,n2,n3, n0,n1,n2,n3, ...], reshape to (Ne, 4)
        raw_connectivity = np.array(elem_node_tags).reshape(-1, 4)

        # Renumber to compact 0-based indices
        connectivity = np.array([[tag_to_compact[tag] for tag in elem] for elem in raw_connectivity], dtype=np.uint32)
        gmsh.finalize()

        mesh_out: Mesh = Mesh(nodes, connectivity)
        print(f"Nodes: {mesh_out.num_nodes}, Elems: {mesh_out.num_elems}\n")

        return mesh_out

    except ImportError:
        raise RuntimeError(f"Error - gmsh is not installed. Could not mesh file `{infile}`") from None

oersted.testing

Utility functions for tests

bz_finite_length_solenoid

bz_finite_length_solenoid(
    jmag: float,
    length: float,
    r: float,
    dr: float,
    z: float,
) -> float

Compute the magnetic field on the axis of a finite-length solenoid

This function assumes that the solenoid dr dimension is small relative to the radius and thickness, and does not correct for finite radial thickness.

Parameters:

Name	Type	Description	Default
jmag	float	(A/m2) magnitude of current density in the solenoid	required
length	float	(m) length of the solenoid	required
r	float	(m) representative radius of the solenoid	required
dr	float	(m) thickness of the solenoid cross section	required
z	float	(m) position along the axis of the solenoid at which the field should be calculated	required

Returns:

Type	Description
float	(T) axial magnetic field

Reference

https://en.wikipedia.org/wiki/Solenoid#Finite_continuous_solenoid (with modifications for current density and finite thickness)

Source code in src/oersted/testing.py

def bz_finite_length_solenoid(jmag: float, length: float, r: float, dr: float, z: float) -> float:
    """Compute the magnetic field on the axis of a finite-length solenoid

    This function assumes that the solenoid `dr` dimension is small relative to the
    radius and thickness, and does not correct for finite radial thickness.

    Args:
        jmag: (A/m2) magnitude of current density in the solenoid
        length: (m) length of the solenoid
        r: (m) representative radius of the solenoid
        dr: (m) thickness of the solenoid cross section
        z: (m) position along the axis of the solenoid at which the field should be calculated

    Returns:
        (T) axial magnetic field

    Reference:
        <https://en.wikipedia.org/wiki/Solenoid#Finite_continuous_solenoid>
        (with modifications for current density and finite thickness)
    """

    a: float = 0.5 * MU0 * jmag * dr
    b: float = (z + 0.5 * length) / np.sqrt(r**2 + (z + 0.5 * length) ** 2)
    c: float = (z - 0.5 * length) / np.sqrt(r**2 + (z - 0.5 * length) ** 2)

    return a * (b - c)

make_helmholtz

make_helmholtz(
    filename: str,
    size: float,
    jmag: None | float = None,
    scale=0.001,
) -> tuple[Mesh, NDArray[float64]]

Make the helmholtz coil test problem

Source code in src/oersted/testing.py

def make_helmholtz(filename: str, size: float, jmag: None | float = None, scale=1e-3) -> tuple[Mesh, NDArray[float64]]:
    """Make the helmholtz coil test problem"""

    # datafile: str = "ring"
    # package_root: Path = Path(__file__).parent.parent.parent.absolute()  # tests is 2 levels up
    ring_mesh: Mesh = mesh_step(filename, size, size, scale)

    # The current mesh is centered on the xy plane and is only one circular ring
    # We need to split the single ring into two rings and assign current densities to the elements
    if jmag is None:
        jmag: float = 100.0e3 / (0.02 * 0.02)
    nodes_upper: NDArray[float64] = ring_mesh.nodes.copy()
    nodes_upper[:, 2] += 0.1  # shift upper coil up
    nodes_lower: NDArray[float64] = nodes_upper.copy()
    nodes_lower[:, 2] -= 0.2  # flip to lower side
    nodes = np.vstack((nodes_upper, nodes_lower))
    connectivity_upper: NDArray[uint32] = ring_mesh.connectivity.copy()
    connectivity_lower: NDArray[uint32] = ring_mesh.connectivity.copy() + uint32(ring_mesh.num_nodes)

    helmholtz_mesh = Mesh(nodes, np.vstack((connectivity_upper, connectivity_lower)))

    jdensity = np.zeros((helmholtz_mesh.num_elems, 3))
    phi = np.atan2(helmholtz_mesh.centroids[:, 1], helmholtz_mesh.centroids[:, 0])
    jdensity[:, 0] = -jmag * np.sin(phi)
    jdensity[:, 1] = jmag * np.cos(phi)

    return (helmholtz_mesh, jdensity)

mean_absolute_error

mean_absolute_error(
    baseline: NDArray[float64],
    measurement: NDArray[float64],
) -> float

Compute the mean absolute error of measurement against baseline

Source code in src/oersted/testing.py

def mean_absolute_error(baseline: NDArray[float64], measurement: NDArray[float64]) -> float:
    """Compute the mean absolute error of `measurement` against `baseline`"""

    assert baseline.shape == measurement.shape
    assert len(baseline.shape) == 1
    n: int = baseline.shape[0]
    assert n > 0

    return (1 / n) * np.sum(baseline - measurement)

mean_relative_error

mean_relative_error(
    baseline: NDArray[float64],
    measurement: NDArray[float64],
) -> float

Compute the mean relative error of measurement against baseline

Source code in src/oersted/testing.py

def mean_relative_error(baseline: NDArray[float64], measurement: NDArray[float64]) -> float:
    """Compute the mean *relative* error of `measurement` against `baseline`"""

    assert baseline.shape == measurement.shape
    assert len(baseline.shape) == 1
    n: int = baseline.shape[0]
    assert n > 0

    relative_diff = (measurement - baseline) / baseline
    return (1 / n) * np.sum(np.abs(relative_diff))

mean_squared_error

mean_squared_error(
    baseline: NDArray[float64],
    measurement: NDArray[float64],
) -> float

Compute the mean squared error of measurement measured against baseline

Error is computed according to this reference: https://en.wikipedia.org/wiki/Mean_squared_error#Predictor

This is an absolute, not a relative error measurement.

---

baseline: N-length array of 'ground-truth' values
measurement: N-length array of values to which determine error (deviation) from baseline

Source code in src/oersted/testing.py

def mean_squared_error(baseline: NDArray[float64], measurement: NDArray[float64]) -> float:
    """Compute the mean squared error of `measurement` measured against `baseline`

    Error is computed according to this reference:
    https://en.wikipedia.org/wiki/Mean_squared_error#Predictor

    This is an absolute, not a relative error measurement.

    Args:
    ---
        baseline: N-length array of 'ground-truth' values
        measurement: N-length array of values to which determine error (deviation) from baseline

    """

    assert baseline.shape == measurement.shape
    assert len(baseline.shape) == 1
    n: int = baseline.shape[0]
    assert n > 0

    return (1 / n) * np.sum((baseline - measurement) ** 2)

smape

smape(
    baseline: NDArray[float64],
    measurement: NDArray[float64],
) -> float

Compute the symmetric mean absolute percentage error of measurement against baseline

SMAPE is defined here: https://en.wikipedia.org/wiki/Symmetric_mean_absolute_percentage_error

Source code in src/oersted/testing.py

def smape(baseline: NDArray[float64], measurement: NDArray[float64]) -> float:
    """Compute the symmetric mean absolute percentage error of `measurement` against `baseline`

    SMAPE is defined here:
    https://en.wikipedia.org/wiki/Symmetric_mean_absolute_percentage_error
    """

    assert baseline.shape == measurement.shape
    assert len(baseline.shape) == 1
    n: int = baseline.shape[0]
    assert n > 0

    numerator = np.abs(measurement - baseline)
    denominator = np.abs(measurement) + np.abs(baseline)
    return (2 / n) * np.sum(numerator / denominator)