Skip to content

Python API

oersted

Python bindings for oersted

BHCurve dataclass

B-H curve for a nonlinear magnetic material

Source code in src/oersted/materials.py 
11
12
13
14
15
16
17
18
19
20
21
@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 
18
19
20
21
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 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
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 the number of elements in the mesh"""
        return self._centroids.shape[0]

    @property
    def centroids(self) -> NDArray[float64]:
        """Return an (N,3) array of the centroid position of each element in the mesh"""
        return self._centroids

    @property
    def volumes(self) -> NDArray[float64]:
        """Return an (N,) array of the volume of each element"""
        return self._volumes

centroids property 

centroids: NDArray[float64]

Return an (N,3) array of the centroid position of each element in the mesh

num_elems property 

num_elems: int

Return the number of elements in the mesh

volumes property 

volumes: NDArray[float64]

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

DirectSolver

Bases: Solver

Controls solver options for using the direct (full integration) solution routines

Source code in src/oersted/solver.py 
46
47
48
49
50
51
52
53
54
55
56
57
class DirectSolver(Solver):
    """Controls solver options for using the direct (full integration)
    solution routines"""

    def __init__(
        self, n_threads: int = 0, max_iterations=100, tol=1.0, alpha=0.5, edge=False
    ):
        self._n_threads = n_threads
        self._max_iterations = max_iterations
        self._tol = tol
        self._alpha = alpha
        self._edge = edge

FreeSpace dataclass

Bases: Material

Source code in src/oersted/materials.py 
38
39
40
41
42
@dataclass
class FreeSpace(Material):
    def mu_r(self, h: float) -> float:
        """Free space (vacuum) has a relative permeability of exactly 1.0"""
        return 1.0

mu_r 

mu_r(h: float) -> float

Free space (vacuum) has a relative permeability of exactly 1.0

Source code in src/oersted/materials.py 
40
41
42
def mu_r(self, h: float) -> float:
    """Free space (vacuum) has a relative permeability of exactly 1.0"""
    return 1.0

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 
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
@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 the relative magnetic permeability"""
        return self._mu_r

    def to_mh_curve(self) -> tuple[NDArray[float64], NDArray[float64]]:
        """Convert the B-H curve to an M-H curve, which is more convenient for
        magnetization calculations"""
        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)

mu_r 

mu_r(h: float)

Return the relative magnetic permeability

Source code in src/oersted/materials.py 
56
57
58
def mu_r(self, h: float):
    """Return the relative magnetic permeability"""
    return self._mu_r

to_mh_curve 

to_mh_curve() -> tuple[NDArray[float64], NDArray[float64]]

Convert the B-H curve to an M-H curve, which is more convenient for magnetization calculations

Source code in src/oersted/materials.py 
60
61
62
63
64
65
def to_mh_curve(self) -> tuple[NDArray[float64], NDArray[float64]]:
    """Convert the B-H curve to an M-H curve, which is more convenient for
    magnetization calculations"""
    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 
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
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
        raise NotImplementedError
        # 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) -> SurfaceMesh:
        """Return the surface mesh associated with the volumetric mesh"""
        if self._surface is None:
            faces = mesh_surface_faces(ascontiguousarray(self.connectivity))
            self._surface = SurfaceMesh(self.nodes.copy(), faces)

        return self._surface

    def plot(self, filename: str | None = None):
        """Convenience function for plotting just the mesh itself"""
        plot_mesh(self, filename)

    @classmethod
    def from_step(
        cls,
        filename: str,
        mesh_size: float,
        mesh_size_scale: float = 1e3,
        part_size_scale: float = 1e-3,
    ) -> Mesh:
        """Create a Mesh from a step file

        !!! note
            `oersted` needs mesh dimensions in meters, but gmsh typically works in mm.
            This function scales the input mesh size *up* to be in mm, and scales the
            gmsh resultant *down* to be in m. Adjust these parameters if the mesh isn't
            working properly.

        Args:
            filename: STEP file to mesh
            mesh_size: (m) nominal element size to use for the mesh
            mesh_size_scale: (mm/m) adjust if the model units are in mm and not m
            part_size_scale: (m/mm) adjust if the mesh units are in mm and not m

        Returns:
            volumetric tet4 mesh of the STEP file
        """
        return mesh_step(
            filename,
            mesh_size * mesh_size_scale,
            mesh_size * mesh_size_scale,
            part_size_scale,
        )

    def append(self, mesh: Mesh) -> Mesh:
        """Convenience function for appending two meshes together."""

        nodes = np.vstack((self.nodes, mesh.nodes))
        connectivity = np.vstack(
            (self.connectivity, mesh.connectivity + uint32(self.num_nodes))
        )
        return Mesh(nodes, connectivity)

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

surface property 

surface: SurfaceMesh

Return the surface mesh associated with the volumetric mesh

volumes property 

volumes: NDArray[float64]

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

append 

append(mesh: Mesh) -> Mesh

Convenience function for appending two meshes together.

Source code in src/oersted/mesh.py 
274
275
276
277
278
279
280
281
def append(self, mesh: Mesh) -> Mesh:
    """Convenience function for appending two meshes together."""

    nodes = np.vstack((self.nodes, mesh.nodes))
    connectivity = np.vstack(
        (self.connectivity, mesh.connectivity + uint32(self.num_nodes))
    )
    return Mesh(nodes, connectivity)

from_step classmethod 

from_step(
    filename: str,
    mesh_size: float,
    mesh_size_scale: float = 1000.0,
    part_size_scale: float = 0.001,
) -> Mesh

Create a Mesh from a step file

Note

oersted needs mesh dimensions in meters, but gmsh typically works in mm. This function scales the input mesh size up to be in mm, and scales the gmsh resultant down to be in m. Adjust these parameters if the mesh isn't working properly.

Parameters:

Name Type Description Default
filename str

STEP file to mesh

 required
mesh_size float

(m) nominal element size to use for the mesh

 required
mesh_size_scale float

(mm/m) adjust if the model units are in mm and not m

 1000.0
part_size_scale float

(m/mm) adjust if the mesh units are in mm and not m

 0.001

Returns:

Type Description
Mesh

volumetric tet4 mesh of the STEP file
Source code in src/oersted/mesh.py 
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
@classmethod
def from_step(
    cls,
    filename: str,
    mesh_size: float,
    mesh_size_scale: float = 1e3,
    part_size_scale: float = 1e-3,
) -> Mesh:
    """Create a Mesh from a step file

    !!! note
        `oersted` needs mesh dimensions in meters, but gmsh typically works in mm.
        This function scales the input mesh size *up* to be in mm, and scales the
        gmsh resultant *down* to be in m. Adjust these parameters if the mesh isn't
        working properly.

    Args:
        filename: STEP file to mesh
        mesh_size: (m) nominal element size to use for the mesh
        mesh_size_scale: (mm/m) adjust if the model units are in mm and not m
        part_size_scale: (m/mm) adjust if the mesh units are in mm and not m

    Returns:
        volumetric tet4 mesh of the STEP file
    """
    return mesh_step(
        filename,
        mesh_size * mesh_size_scale,
        mesh_size * mesh_size_scale,
        part_size_scale,
    )

plot 

plot(filename: str | None = None)

Convenience function for plotting just the mesh itself

Source code in src/oersted/mesh.py 
238
239
240
def plot(self, filename: str | None = None):
    """Convenience function for plotting just the mesh itself"""
    plot_mesh(self, filename)

to_centroid_mesh 

to_centroid_mesh() -> CentroidMesh

Create a centroid mesh from a tet4 mesh

Source code in src/oersted/mesh.py 
178
179
180
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 
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
@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:
        """Return the relative magnetic permeability at a given H-field strength"""
        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]]:
        """Convert the B-H curve to an M-H curve, which is more convenient for
        magnetization calculations"""
        return (self.curve.h_values, self.curve.b_values / MU0 - self.curve.h_values)

mu_r 

mu_r(h: float) -> float

Return the relative magnetic permeability at a given H-field strength

Source code in src/oersted/materials.py 
77
78
79
80
81
82
83
84
def mu_r(self, h: float) -> float:
    """Return the relative magnetic permeability at a given H-field strength"""
    b: float = self.curve.lookup(h)

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

to_mh_curve 

to_mh_curve() -> tuple[NDArray[float64], NDArray[float64]]

Convert the B-H curve to an M-H curve, which is more convenient for magnetization calculations

Source code in src/oersted/materials.py 
86
87
88
89
def to_mh_curve(self) -> tuple[NDArray[float64], NDArray[float64]]:
    """Convert the B-H curve to an M-H curve, which is more convenient for
    magnetization calculations"""
    return (self.curve.h_values, self.curve.b_values / MU0 - self.curve.h_values)

OctreeSolver

Bases: Solver

Controls solver options for using the octree (Barnes-Hut) solution routines

This uses the "new" 3-zone octree solver.

Warning

There is currently a bug with leaf_threshold > 1, so the solver will always use a leaf threshold of 1.

Source code in src/oersted/solver.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
class OctreeSolver(Solver):
    """Controls solver options for using the octree (Barnes-Hut) solution routines

    This uses the "new" 3-zone octree solver.

    !!! warning
        There is currently a bug with `leaf_threshold` > 1, so the solver will
        always use a leaf threshold of 1.
    """

    _theta: float
    _leaf_threshold: int
    _n_threads: int
    _max_iterations: int
    _tol: float
    _alpha: float
    _edge: bool

    def __init__(
        self,
        theta: float = 0.25,
        alpha: float = 3.0,
        leaf_threshold: int = 1,
        n_threads: int = 0,
        max_iterations=100,
        tol=1.0,
    ):
        self._theta = theta
        self._alpha = alpha
        self._n_threads = n_threads
        self._max_iterations = max_iterations
        self._tol = tol
        self._leaf_threshold = 1
        self._edge = False

    @property
    def theta(self):
        """Returns the Barnes-Hut angle-opening criteria (accuracy control)"""
        return self._theta

    @property
    def alpha(self):
        """Returns the near/mid-field element radius ratio (accuracy control)"""
        return self._alpha

    @property
    def leaf_threshold(self):
        """Returns the number of sources that will be evaluated individually at the
        octree leaf level"""
        return self._leaf_threshold

alpha property 

alpha

Returns the near/mid-field element radius ratio (accuracy control)

leaf_threshold property 

leaf_threshold

Returns the number of sources that will be evaluated individually at the octree leaf level

theta property 

theta

Returns the Barnes-Hut angle-opening criteria (accuracy control)

OctreeSolver2Zone

Bases: Solver

Controls solver options for using the octree (Barnes-Hut) solution routines

This uses the "old" 2-zone octree solver.

Source code in src/oersted/solver.py 
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
class OctreeSolver2Zone(Solver):
    """Controls solver options for using the octree (Barnes-Hut) solution routines

    This uses the "old" 2-zone octree solver.
    """

    _theta: float
    _leaf_threshold: int
    _n_threads: int
    _max_iterations: int
    _tol: float
    _alpha: float
    _edge: bool

    def __init__(
        self,
        theta: float = 0.25,
        leaf_threshold: int = 16,
        n_threads: int = 0,
        max_iterations=100,
        tol=1.0,
        alpha=0.5,
        edge=False,
    ):
        self._theta = theta
        self._leaf_threshold = leaf_threshold
        self._n_threads = n_threads
        self._max_iterations = max_iterations
        self._tol = tol
        self._alpha = alpha
        self._edge = False

    @property
    def theta(self):
        """Returns the Barnes-Hut angle-opening criteria (accuracy control)"""
        return self._theta

    @property
    def leaf_threshold(self):
        """Returns the number of sources that will be evaluated individually at the
        octree leaf level"""
        return self._leaf_threshold

leaf_threshold property 

leaf_threshold

Returns the number of sources that will be evaluated individually at the octree leaf level

theta property 

theta

Returns the Barnes-Hut angle-opening criteria (accuracy control)

Solver

Defines a generic solver with options for iterative and multi-threaded solutions

Source code in src/oersted/solver.py 
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
class Solver:
    """Defines a generic solver with options for iterative and multi-threaded
    solutions"""

    _n_threads: int
    _max_iterations: int
    _tol: float
    _alpha: float
    _edge: bool

    @property
    def n_threads(self):
        """Number of threads for solving."""
        return self._n_threads

    @property
    def max_iterations(self):
        """Maximum number of iterations for iterative solves."""
        return self._max_iterations

    @property
    def tol(self):
        """Convergence tolerance for iterative solves.

        !!! note
            This is typically a value of the H-field for magnetization calculations;
            the default value is roughly 1.3e-6 T, which may be more than necessary
            for most applications.
        """
        return self._tol

    @property
    def alpha(self):
        """Under-relaxation factor for iterative solves."""
        return self._alpha

    @property
    def edge(self):
        """Whether or not the edge-based element calculation should be used"""
        return self._edge

alpha property 

alpha

Under-relaxation factor for iterative solves.

edge property 

edge

Whether or not the edge-based element calculation should be used

max_iterations property 

max_iterations

Maximum number of iterations for iterative solves.

n_threads property 

n_threads

Number of threads for solving.

tol property 

tol

Convergence tolerance for iterative solves.

Note

This is typically a value of the H-field for magnetization calculations; the default value is roughly 1.3e-6 T, which may be more than necessary for most applications.

b_field 

b_field(
    source: Mesh | CentroidMesh,
    j_density: NDArray[float64],
    targets: NDArray[float64],
    solver: DirectSolver
    | OctreeSolver
    | OctreeSolver2Zone
    | 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

Parameters:

Name Type Description Default
source Mesh | CentroidMesh

mesh to use as the field source

 required
j_density NDArray[float64]

(A/m^2) (N,3) array of current density vectors at each of the element centroids

 required
targets NDArray[float64]

(m) (N,3) array of target point positions in 3D space

 required
solver DirectSolver | OctreeSolver | OctreeSolver2Zone | None

selects the solution settings

 None

Returns:

Type Description
NDArray[float64]

(T) (N,3) array of magnetic flux density (B) vectors at each target position
Source code in src/oersted/biotsavart.py 
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
def b_field(
    source: Mesh | CentroidMesh,
    j_density: NDArray[float64],
    targets: NDArray[float64],
    solver: DirectSolver | OctreeSolver | OctreeSolver2Zone | 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

    Args:
        source: mesh to use as the field source
        j_density: (A/m^2) (N,3) array of current density vectors at each of the
            element centroids
        targets: (m) (N,3) array of target point positions in 3D space
        solver: selects the solution settings

    Returns:
        (T) (N,3) array of magnetic flux density (B) vectors at each target position
    """
    return MU0 * h_field(source, j_density, targets, solver)

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 
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
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)

bz_loop_axis 

bz_loop_axis(
    current: float, radius: float, z: float
) -> float

Compute the vertical field Bz at the center of a current-carrying loop

Parameters:

Name Type Description Default
current float

(A) total electrical current in the loop

 required
radius float

(m) centerline radius of the loop

 required
z float

(m) height of the target point along the loop axis

 required

Returns:

Type Description
float

(T) magnetic flux density, oriented along the loop axis
Source code in src/oersted/testing.py 
170
171
172
173
174
175
176
177
178
179
180
181
def bz_loop_axis(current: float, radius: float, z: float) -> float:
    """Compute the vertical field Bz at the center of a current-carrying loop

    Args:
        current: (A) total electrical current in the loop
        radius: (m) centerline radius of the loop
        z: (m) height of the target point along the loop axis

    Returns:
        (T) magnetic flux density, oriented along the loop axis
    """
    return MU0 * current * (radius**2) / (2.0 * (z**2 + radius**2) ** 1.5)

dbzdz_loop_axis 

dbzdz_loop_axis(
    current: float, radius: float, z: float
) -> float

Compute the vertical field gradient dBz/dz at the center of a current-carrying loop

Parameters:

Name Type Description Default
current float

(A) total electrical current in the loop

 required
radius float

(m) centerline radius of the loop

 required
z float

(m) height of the target point along the loop axis

 required

Returns:

Type Description
float

(T/m) magnetic flux density gradient, oriented along the loop axis
Source code in src/oersted/testing.py 
184
185
186
187
188
189
190
191
192
193
194
195
196
def dbzdz_loop_axis(current: float, radius: float, z: float) -> float:
    """Compute the vertical field gradient dBz/dz at the center of a
    current-carrying loop

    Args:
        current: (A) total electrical current in the loop
        radius: (m) centerline radius of the loop
        z: (m) height of the target point along the loop axis

    Returns:
        (T/m) magnetic flux density gradient, oriented along the loop axis
    """
    return -1.5 * MU0 * current * (radius**2) * z / (z**2 + radius**2) ** 2.5

demag_solve 

demag_solve(
    mesh: Mesh,
    material: Material,
    h_external: NDArray[float64],
    solver: DirectSolver | OctreeSolver | OctreeSolver2Zone,
) -> tuple[NDArray[float64], NDArray[float64]]

Compute magnetization field M and the total H field at element centroids, given a background field

Uses simple fixed-point iteration and therefore only converges for low-permeable materials.

Parameters:

Name Type Description Default
mesh Mesh

finite element mesh on which to evaluate the demagnetizing field

 required
material Material

linear or nonlinear magnetic maaterial properties

 required
h_external NDArray[float64]

(A/m) an (Ne,3) array of external field at each element centroid

 required
solver DirectSolver | OctreeSolver | OctreeSolver2Zone

solution parameters for the problem, including iteration method

 required

Returns:

Type Description
(M, Htotal)

each (Ne, 3), magnetization field M(Htotal) and total H field
NDArray[float64]

at element centroids. These can be summed to give B = mu0 * (Htotal + M).
Source code in src/oersted/magnetization.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
def demag_solve(
    mesh: Mesh,
    material: Material,
    h_external: NDArray[float64],
    solver: DirectSolver | OctreeSolver | OctreeSolver2Zone,
) -> tuple[NDArray[float64], NDArray[float64]]:
    """Compute magnetization field M and the total H field at element centroids,
        given a background field

    Uses simple fixed-point iteration and therefore only converges for
        low-permeable materials.

    Args:
        mesh: finite element mesh on which to evaluate the demagnetizing field
        material: linear or nonlinear magnetic maaterial properties
        h_external: (A/m) an (Ne,3) array of external field at each element centroid
        solver: solution parameters for the problem, including iteration method

    Returns:
        (M, Htotal): each (Ne, 3), magnetization field M(Htotal) and total H field
        at element centroids. These can be summed to give B = mu0 * (Htotal + M).
    """

    theta: float
    leaf_threshold: uint32

    if isinstance(solver, DirectSolver):
        theta = 0.5
        leaf_threshold = uint32(0)
    elif isinstance(solver, OctreeSolver):
        raise NotImplementedError(
            f"{type(solver)} not implemented yet for demag solve."
        )
    else:
        theta = solver.theta
        leaf_threshold = solver.leaf_threshold

    return magnetization_tet4(
        ascontiguousarray(mesh.nodes),
        ascontiguousarray(mesh.connectivity),
        material.chi(1.0),
        ascontiguousarray(h_external),
        solver.tol,
        solver.max_iterations,
        theta,
        leaf_threshold,
        solver.alpha,
        solver.n_threads,
        solver.edge,
    )

h_field 

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

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

Parameters:

Name Type Description Default
source Mesh | CentroidMesh

mesh to use as the field source

 required
j_density NDArray[float64]

(A/m^2) (N,3) array of current density vectors at each of the element centroids

 required
targets NDArray[float64]

(m) (N,3) array of target point positions in 3D space

 required
solver DirectSolver | OctreeSolver | OctreeSolver2Zone | None

selects the solution settings

 None

Returns:

Type Description
NDArray[float64]

(T) (N,3) array of magnetic field strength (H) vectors at each target position
Source code in src/oersted/biotsavart.py 
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
def h_field(
    source: Mesh | CentroidMesh,
    j_density: NDArray[float64],
    targets: NDArray[float64],
    solver: DirectSolver | OctreeSolver | OctreeSolver2Zone | None = None,
    edge: bool = False,
) -> NDArray[float64]:
    """Compute the magnetic field strength at a collection of target points using
    a current-carrying source mesh.

    Args:
        source: mesh to use as the field source
        j_density: (A/m^2) (N,3) array of current density vectors at each of the
            element centroids
        targets: (m) (N,3) array of target point positions in 3D space
        solver: selects the solution settings

    Returns:
        (T) (N,3) array of magnetic field strength (H) vectors at each target position
    """

    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 h_current_point_direct(
                src_pts, src_vol, j_density, tgt_pts, solver.n_threads
            )

        elif isinstance(solver, OctreeSolver2Zone):
            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, edge
            )

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

        elif isinstance(solver, OctreeSolver):
            return h_current_octree(
                src_nodes,
                src_connectivity,
                tgt_pts,
                j_density,
                uint32(solver.leaf_threshold),
                solver.alpha,
                solver.theta,
                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
    | OctreeSolver2Zone
    | None = None,
) -> NDArray[float64]

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

Parameters:

Name Type Description Default
source Mesh | CentroidMesh

mesh to use as the field source

 required
m_field NDArray[float64]

(A/m) (N,3) array of magnetization field vectors at each of the element centroids

 required
targets NDArray[float64]

(m) (N,3) array of target point positions in 3D space

 required
solver DirectSolver | OctreeSolver | OctreeSolver2Zone | None

selects the solution settings

 None

Returns:

Type Description
NDArray[float64]

(T) (N,3) array of magnetic field strength (H) vectors at each target position
Source code in src/oersted/biotsavart.py 
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
def h_mag(
    source: Mesh | CentroidMesh,
    m_field: NDArray[float64],
    targets: NDArray[float64],
    solver: DirectSolver | OctreeSolver | OctreeSolver2Zone | None = None,
) -> NDArray[float64]:
    """Compute the magnetic field strength using a magnetized mesh as the source

    Args:
        source: mesh to use as the field source
        m_field: (A/m) (N,3) array of magnetization field vectors at each of the
            element centroids
        targets: (m) (N,3) array of target point positions in 3D space
        solver: selects the solution settings

    Returns:
        (T) (N,3) array of magnetic field strength (H) vectors at each target position
    """

    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):
        src_centroids = ascontiguousarray(source.centroids, dtype=float64)
        src_volumes = ascontiguousarray(source.volumes, dtype=float64)

        assert source.centroids.shape[0] == m_field.shape[0]

        if isinstance(solver, DirectSolver):
            theta = 0.0
            leaf_threshold: uint32 = uint32(0)
            use_octree = False

            return h_mag_point(
                src_centroids,
                src_volumes,
                m_field,
                targets,
                theta,
                leaf_threshold,
                solver.n_threads,
                use_octree,
            )

        elif isinstance(solver, OctreeSolver2Zone):
            use_octree = True
            return h_mag_point(
                src_centroids,
                src_volumes,
                m_field,
                targets,
                solver.theta,
                solver.leaf_threshold,
                solver.n_threads,
                use_octree,
            )
        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)

        assert src_connectivity.shape[0] == m_field.shape[0]
        if isinstance(solver, DirectSolver):
            theta = 0.0
            leaf_threshold: uint32 = uint32(0)
            use_octree = False
            return h_mag_tet4(
                src_nodes,
                src_connectivity,
                m_field,
                targets,
                theta,
                leaf_threshold,
                solver.n_threads,
                use_octree,
                solver.edge,
            )

        elif isinstance(solver, OctreeSolver2Zone):
            use_octree = True
            return h_mag_tet4(
                src_nodes,
                src_connectivity,
                m_field,
                targets,
                solver.theta,
                solver.leaf_threshold,
                solver.n_threads,
                use_octree,
                solver.edge,
            )

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

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

kelvin_force_density 

kelvin_force_density(
    mesh: Mesh,
    m_field_centroids: NDArray[float64],
    b_field_nodes: NDArray[float64],
) -> NDArray[float64]

Compute the Kelvin force density acting on a magnetized mesh.

Note

The h_field must be calculated at the element nodes, while the magnetization field must be known at the centroids.

Parameters:

Name Type Description Default
mesh Mesh

the mesh on which to evaluate the forces

 required
m_field_centroids NDArray[float64]

(A/m) the magnetization field (M field), evaluated at the centroids of the elements in the mesh

 required
b_field_nodes NDArray[float64]

(T) the magnetic flux density, evaluated at the nodes within the mesh

 required

Returns:

Type Description
NDArray[float64]

(N/m^3) an (N,3) array of the force density vector acting on each element centroid
Source code in src/oersted/results.py 
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def kelvin_force_density(
    mesh: Mesh, m_field_centroids: NDArray[float64], b_field_nodes: NDArray[float64]
) -> NDArray[float64]:
    """Compute the Kelvin force density acting on a magnetized mesh.

    !!! note
        The h_field must be calculated at the element nodes, while the magnetization
        field must be known at the centroids.

    Args:
        mesh: the mesh on which to evaluate the forces
        m_field_centroids: (A/m) the magnetization field (M field), evaluated at the
            centroids of the elements in the mesh
        b_field_nodes: (T) the magnetic flux density, evaluated at the nodes within
            the mesh

    Returns:
        (N/m^3) an (N,3) array of the force density vector acting on each element
            centroid
    """

    # Check that the inputs are defined properly
    assert mesh.centroids.shape == m_field_centroids.shape
    assert mesh.nodes.shape == b_field_nodes.shape

    return mesh_kelvin_force_density(
        ascontiguousarray(mesh.nodes),
        ascontiguousarray(mesh.connectivity),
        ascontiguousarray(m_field_centroids),
        ascontiguousarray(b_field_nodes),
    )

kelvin_forces 

kelvin_forces(
    mesh: Mesh,
    m_field_centroids: NDArray[float64],
    b_field_nodes: NDArray[float64],
) -> NDArray[float64]

Compute the Kelvin forces acting on a magnetized mesh.

Note

The h_field must be calculated at the element nodes, while the magnetization field must be known at the centroids.

Parameters:

Name Type Description Default
mesh Mesh

the mesh on which to evaluate the forces

 required
m_field_centroids NDArray[float64]

(A/m) the magnetization field (M field), evaluated at the centroids of the elements in the mesh

 required
b_field_nodes NDArray[float64]

(T) the magnetic flux density, evaluated at the nodes within the mesh

 required

Returns:

Type Description
NDArray[float64]

(N) an (N,3) array of the force vector acting on each element centroid
Source code in src/oersted/results.py 
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def kelvin_forces(
    mesh: Mesh, m_field_centroids: NDArray[float64], b_field_nodes: NDArray[float64]
) -> NDArray[float64]:
    """Compute the Kelvin forces acting on a magnetized mesh.

    !!! note
        The h_field must be calculated at the element nodes, while the magnetization
        field must be known at the centroids.

    Args:
        mesh: the mesh on which to evaluate the forces
        m_field_centroids: (A/m) the magnetization field (M field), evaluated at the
            centroids of the elements in the mesh
        b_field_nodes: (T) the magnetic flux density, evaluated at the nodes within
            the mesh

    Returns:
        (N) an (N,3) array of the force vector acting on each element centroid
    """

    # Check that the inputs are defined properly
    assert mesh.centroids.shape == m_field_centroids.shape
    assert mesh.nodes.shape == b_field_nodes.shape

    return (
        kelvin_force_density(mesh, m_field_centroids, b_field_nodes)
        * mesh.volumes[:, newaxis]
    )

lorentz_force_density 

lorentz_force_density(
    j_density: NDArray[float64], b_field: NDArray[float64]
) -> NDArray[float64]

Compute the Lorentz forces acting on a mesh

Parameters:

Name Type Description Default
j_density NDArray[float64]

(A/m^2) an (N,3) array containing the current density vector at every element

 required
b_field NDArray[float64]

(T) an (N,3) array containing the magnetic flux density vector at every element

 required

Returns:

Type Description
NDArray[float64]

(N/m^3) an (N,3) array of the force density vectors acting on the mesh at every element
Source code in src/oersted/results.py 
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
def lorentz_force_density(
    j_density: NDArray[float64], b_field: NDArray[float64]
) -> NDArray[float64]:
    """Compute the Lorentz forces acting on a mesh

    Args:
        j_density: (A/m^2) an (N,3) array containing the current density vector
            at every element
        b_field: (T) an (N,3) array containing the magnetic flux density vector
            at every element

    Returns:
        (N/m^3) an (N,3) array of the force density vectors acting on the mesh
            at every element
    """

    assert j_density.shape == b_field.shape
    assert j_density.shape[1] == 3

    # This might introduce a copy, though shouldn't be a big hit to the user's
    # observed performance
    jxb: NDArray[float64] = cross(j_density, b_field).astype(float64)

    return jxb

lorentz_forces 

lorentz_forces(
    mesh: Mesh | CentroidMesh,
    j_density: NDArray[float64],
    b_field: NDArray[float64],
    total: bool = False,
) -> NDArray[float64]

Compute the Lorentz forces acting on a mesh

Parameters:

Name Type Description Default
mesh Mesh | CentroidMesh

the finite element mesh on which the fields and forces are calculated

 required
j_density NDArray[float64]

(A/m^2) an (N,3) array containing the current density vector at every element

 required
b_field NDArray[float64]

(T) an (N,3) array containing the magnetic flux density vector at every element

 required
total bool

if true, return only the total force (default: false)

 False

Returns:

Type Description
NDArray[float64]

(N) an (N,3) array of the force vectors acting on the mesh at every element
Source code in src/oersted/results.py 
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
def lorentz_forces(
    mesh: Mesh | CentroidMesh,
    j_density: NDArray[float64],
    b_field: NDArray[float64],
    total: bool = False,
) -> NDArray[float64]:
    """Compute the Lorentz forces acting on a mesh

    Args:
        mesh: the finite element mesh on which the fields and forces are calculated
        j_density: (A/m^2) an (N,3) array containing the current density vector
            at every element
        b_field: (T) an (N,3) array containing the magnetic flux density vector
            at every element
        total: if true, return only the total force (default: false)

    Returns:
        (N) an (N,3) array of the force vectors acting on the mesh at every element
    """

    assert j_density.shape[0] == mesh.num_elems

    forces: NDArray[float64] = (
        lorentz_force_density(j_density, b_field) * mesh.volumes[:, newaxis]
    )

    if total:
        return sum(forces, axis=0)

    else:
        return forces

make_helmholtz 

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

Make the helmholtz coil test problem

The geometry of the problem is defined as:

  • Two circular coils of radius R=0.2m
  • Distance between the coils is d=0.2m
  • Coils are aligned with the z-axis and symmetric about the xy plane
  • Currents in the currents are flowing in the same direction

Parameters:

Name Type Description Default
filename str

STEP file containing the geometry of a single loop in the Helmholtz pair

 required
size float

(m) mesh size

 required
jmag None | float

(A/m^2) magnitude of current density in each loop

 None
scale float

(mm/m) scale factor for meshing

 0.001

Returns:

Type Description
tuple[Mesh, NDArray[float64]]

(mesh, current density) the volumetric mesh and current density vectors of the Helmholtz coil pair
Source code in src/oersted/testing.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
def make_helmholtz(
    filename: str, size: float, jmag: None | float = None, scale: float = 1e-3
) -> tuple[Mesh, NDArray[float64]]:
    """Make the helmholtz coil test problem

    The geometry of the problem is defined as:

    * Two circular coils of radius R=0.2m
    * Distance between the coils is d=0.2m
    * Coils are aligned with the z-axis and symmetric about the xy plane
    * Currents in the currents are flowing in the same direction

    Args:
        filename: STEP file containing the geometry of a single loop in the Helmholtz
            pair
        size: (m) mesh size
        jmag: (A/m^2) magnitude of current density in each loop
        scale: (mm/m) scale factor for meshing

    Returns:
        (mesh, current density) the volumetric mesh and current density vectors of the
            Helmholtz coil pair

    """

    ring_mesh: Mesh = Mesh.from_step(filename, size, 1e3, scale)

    # The current mesh is centered on the xy plane and is only one circular ring
    # 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)

maxwell_forces 

maxwell_forces(
    mesh: SurfaceMesh, b_field: NDArray[float64]
) -> NDArray[float64]

Compute the maxwell stress tensor and determine the force vector acting on each surface face centroid. Returns an (N,3) array of the force vector

Warning

This form of the Maxwell stress tensor is numerically unstable, especially in the presence of large background fields. Perform a mesh convergence study to identify an appropriate mesh.

Parameters:

Name Type Description Default
mesh SurfaceMesh

a surface mesh on which to integrate the Maxwell stress tensor

 required
b_field NDArray[float64]

(T) magnetic flux density evaluated at each of the centroids of the surface mesh elements

 required

Returns:

Type Description
NDArray[float64]

(N) an (N,3) array of the force acting on each of the surface mesh centroids
Source code in src/oersted/results.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
def maxwell_forces(mesh: SurfaceMesh, b_field: NDArray[float64]) -> NDArray[float64]:
    """Compute the maxwell stress tensor and determine the force vector acting on each
    surface face centroid. Returns an (N,3) array of the force vector

    !!! warning
        This form of the Maxwell stress tensor is numerically unstable, especially in
        the presence of large background fields. Perform a mesh convergence study to
        identify an appropriate mesh.

    Args:
        mesh: a surface mesh on which to integrate the Maxwell stress tensor
        b_field: (T) magnetic flux density evaluated at each of the centroids of the
            surface mesh elements

    Returns:
        (N) an (N,3) array of the force acting on each of the surface mesh centroids
    """

    assert mesh.centroids.shape == b_field.shape

    return mesh_surface_forces(mesh.areas, mesh.normals, b_field)

mesh_step 

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

Mesh a step file using gmsh

Note

This requires gmsh to be installed: pip install gmsh

Parameters:

Name Type Description Default
infile str

path to the STEP file to mesh

 required
max_size float

(m) maximum allowable element size

 required
min_size float

(m) minimum allowable element size

 0.0
scale float

(mm/m) adjust if the part or mesh is scaled incorrectly

 0.001

Returns:

Type Description
Mesh

a tet4 (volumetric) mesh of the component
Source code in src/oersted/mesh.py 
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
def mesh_step(
    infile: str, max_size: float, min_size: float = 0.0, scale: float = 1e-3
) -> Mesh:
    """Mesh a step file using gmsh

    !!! note
        This requires `gmsh` to be installed: `pip install gmsh`

    Args:
        infile: path to the STEP file to mesh
        max_size: (m) maximum allowable element size
        min_size: (m) minimum allowable element size
        scale: (mm/m) adjust if the part or mesh is scaled incorrectly

    Returns:
        a tet4 (volumetric) mesh of the component
    """

    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

plot_mesh 

plot_mesh(
    mesh: Mesh,
    filename: str | None = None,
    scalars: NDArray[float64] | None = None,
    centroids: NDArray[float64] | None = None,
    vectors: NDArray[float64] | None = None,
    vector_scale: float | None = None,
    transparency: bool = False,
)

Make a 3D plot of the mesh

Note

This function requires pyvista: pip install pyvista.

Parameters:

Name Type Description Default
mesh Mesh

the tet4 mesh to plot

 required
filename str | None

if this argument is passed, save to file only

 None
scalars NDArray[float64] | None

(Ne,) array of scalar values to color the mesh, defined at element centroids

 None
centroids NDArray[float64] | None

(N,3) array of vector positions for plotting vector values on the mesh

 None
vectors NDArray[float64] | None

(N,3) array of vector magnitudes for plotting vector values on the mesh; must be same length as centroids

 None
vector_scale float | None

adjust for setting vector length

 None
transparency bool

set to True to plot a wireframe of the mesh

 False
Source code in src/oersted/mesh.py 
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
def plot_mesh(
    mesh: Mesh,
    filename: str | None = None,
    scalars: NDArray[float64] | None = None,
    centroids: NDArray[float64] | None = None,
    vectors: NDArray[float64] | None = None,
    vector_scale: float | None = None,
    transparency: bool = False,
):
    """Make a 3D plot of the mesh

    !!! note
        This function requires `pyvista`: `pip install pyvista`.

    Args:
        mesh: the tet4 mesh to plot
        filename: if this argument is passed, save to file only
        scalars: (Ne,) array of scalar values to color the mesh, defined at element
            centroids
        centroids: (N,3) array of vector positions for plotting vector values on the
            mesh
        vectors: (N,3) array of vector magnitudes for plotting vector values on the
            mesh; must be same length as `centroids`
        vector_scale: adjust for setting vector length
        transparency: set to `True` to plot a wireframe of the mesh
    """

    try:
        import pyvista as pv

        cells = np.hstack(
            [np.full((mesh.num_elems, 1), 4, dtype=np.int64), mesh.connectivity]
        )
        celltypes = np.full(mesh.num_elems, pv.CellType.TETRA)
        pv_mesh = pv.UnstructuredGrid(cells.ravel(), celltypes, mesh.nodes)

        pl = pv.Plotter(off_screen=filename is not None)
        if transparency:
            pl.add_mesh(pv_mesh, style="wireframe", color="black", line_width=0.5)
        else:
            pl.add_mesh(
                pv_mesh,
                scalars=scalars,
                show_edges=True,
                line_width=0.5,
                scalar_bar_args={"title": "magnitude", "vertical": True},
            )

        factor = 1.0
        if vector_scale is not None:
            factor = vector_scale

        if centroids is not None and vectors is not None:
            assert centroids.shape == vectors.shape
            arrow_mesh = pv.PolyData(centroids)
            arrow_mesh["vectors"] = vectors
            arrow_mesh["magnitude"] = np.linalg.norm(vectors, axis=1)
            arrows = arrow_mesh.glyph(orient="vectors", scale=False, factor=factor)
            pl.add_mesh(
                arrows, scalars="magnitude", cmap="viridis", show_scalar_bar=False
            )

        if filename is not None:
            pl.save_graphic(filename)
        else:
            pl.show()

    except ImportError:
        print("`pyvista` is not installed.")
        print("Please install before continuing: `pip install pyvista")

smape 

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

Compute the symmetric mean absolute percentage error of measurement against baseline

Parameters:

Name Type Description Default
baseline NDArray[float64]

an (N,) array of data values to compute error against

 required
measurement NDArray[float64]

an (N,) array of data values of which to compute error

 required

Returns:

Type Description
float

error associated with this comparison

Reference

https://en.wikipedia.org/wiki/Symmetric_mean_absolute_percentage_error

Source code in src/oersted/testing.py 
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
def smape(baseline: NDArray[float64], measurement: NDArray[float64]) -> float:
    """Compute the symmetric mean absolute percentage error of `measurement`
        against `baseline`

    Args:
        baseline: an (N,) array of data values to compute error against
        measurement: an (N,) array of data values of which to compute error

    Returns:
        error associated with this comparison

    !!! reference
        <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)