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

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 
 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
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
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

    @nodes.setter
    def nodes(self, new_nodes: NDArray[float64]):
        assert self._nodes.shape == new_nodes.shape
        self._nodes = new_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]

    @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 | Path,
        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(
            str(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 writable 

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 
245
246
247
248
249
250
251
252
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 | Path,
    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 | Path

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 
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
@classmethod
def from_step(
    cls,
    filename: str | Path,
    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(
        str(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 
209
210
211
def plot(self, filename: str | None = None):
    """Convenience function for plotting just the mesh itself"""
    plot_mesh(self, filename)

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)

SolverSettings dataclass

Define settings to use during a solve

Attributes:

Name Type Description
method Method

use either "direct" O(N^2) or "octree" O(N log(N)) integration for Biot-Savart law integration
integration Integration

use either "element" for full integration (high accuracy, slow) or "point" for a point-source approximation (low accuracy near the source, but ~40x faster)
n_threads int

number of cpu threads to use for solution; default is 0 for all available
theta float

the Barnes-Hut angle-opening criteria (theta >= 0.0)
near_field_ratio float

(deprecated) ratio of target distance to source element size; defines the 'mid-field' during a barnes-hut solve (alpha >= 0.0)
max_leaf_size int

defines the maximum leaf size before splitting in the octree (>=1)
batch_size int

defines the number of target points to process at a time (for Barnes-Hut solves))
multipole_order MultipoleOrder

defines the multipole expansion order for Barnes-Hut solves; currently supported are "monopole" and "dipole"
max_iterations int

maximum number of iterations before solution returns early
atol float

absolute tolerance for convergence criteria
under_relaxation_factor float

provides solution stability for fixed-point iteration
Source code in src/oersted/solver.py 
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
44
45
46
47
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
@dataclass(frozen=True, slots=True, kw_only=True)
class SolverSettings:
    """Define settings to use during a solve

    Attributes:
        method: use either "direct" O(N^2) or "octree" O(N log(N)) integration for
            Biot-Savart law integration
        integration: use either "element" for full integration (high
            accuracy, slow) or "point" for a point-source approximation (low accuracy
            near the source, but ~40x faster)
        n_threads: number of cpu threads to use for solution; default is 0 for all
            available
        theta: the Barnes-Hut angle-opening criteria (theta >= 0.0)
        near_field_ratio: (deprecated) ratio of target distance to source element size;
            defines the 'mid-field' during a barnes-hut solve (alpha >= 0.0)
        max_leaf_size: defines the maximum leaf size before splitting in the octree
            (>=1)
        batch_size: defines the number of target points to process at a time (for
            Barnes-Hut solves))
        multipole_order: defines the multipole expansion order for Barnes-Hut solves;
            currently supported are "monopole" and "dipole"
        max_iterations: maximum number of iterations before solution returns early
        atol: absolute tolerance for convergence criteria
        under_relaxation_factor: provides solution stability for fixed-point iteration
    """

    # General settings
    method: Method = "direct"
    integration: Integration = "element"
    n_threads: int = 0

    # Octree settings
    theta: float = 0.5
    max_leaf_size: int = 16
    multipole_order: MultipoleOrder = "dipole"
    batch_size: int = 10_000

    # Iterative solve settings
    max_iterations: int = 100
    atol: float = 1e-6
    under_relaxation_factor: float = 0.5
    verbose: bool = False

    def __post_init__(self):
        if self.method not in get_args(Method):
            raise ValueError(
                f"Solution method must be one of {get_args(Method)}, \
                    got {self.method!r}"
            )
        if self.integration not in get_args(Integration):
            raise ValueError(
                f"Near-field integration method must be one of \
                    {get_args(Integration)}, got `{self.integration!r}"
            )
        if self.n_threads < 0:
            raise ValueError(
                f"Number of cpu threads must be greater than 1, got `{self.n_threads}`"
            )
        if self.theta < 0.0:
            raise ValueError(f"theta must be > 0, got {self.theta}")
        if self.max_leaf_size < 1:
            raise ValueError(f"Max leaf size must be >= 1, got {self.max_leaf_size}")
        if self.max_iterations < 1:
            raise ValueError(
                f"Max solver iterations must be >= 1, got `{self.max_iterations}"
            )
        if self.atol < 0:
            raise ValueError(f"Relative error tolerance must be > 0, got `{self.atol}")
        if self.under_relaxation_factor > 1.0 or self.under_relaxation_factor < 0.0:
            raise ValueError(
                f"Under relaxation factor must be in range [0,1], got \
                    `{self.under_relaxation_factor}"
            )

    def expansion_order(self) -> int:
        if self.multipole_order == "dipole":
            return 1
        elif self.multipole_order == "octupole":
            return 2
        else:
            return 0

a_field 

a_field(
    src_mesh: Mesh,
    targets: NDArray[float64],
    *,
    jdensity: NDArray[float64] | None = None,
    magnetization: NDArray[float64] | None = None,
    settings: SolverSettings = DEFAULT_SETTINGS,
) -> NDArray[float64]

Compute the magnetic vector potential (A field) at a collection of target points

Parameters:

Name Type Description Default
src_mesh Mesh

mesh to use as the field source

 required
targets NDArray[float64]

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

 required
jdensity NDArray[float64] | None

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

 None
magnetization NDArray[float64] | None

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

 None
settings SolverSettings

selects the solution settings

 DEFAULT_SETTINGS

Returns:

Type Description
NDArray[float64]

(T-m) (N,3) array of magnetic vector potential (A) vectors at each target position
Source code in src/oersted/biotsavart.py 
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
def a_field(
    src_mesh: Mesh,
    targets: NDArray[float64],
    *,  # Force remaining variables to be passed by name
    jdensity: NDArray[float64] | None = None,
    magnetization: NDArray[float64] | None = None,
    settings: SolverSettings = DEFAULT_SETTINGS,
) -> NDArray[float64]:
    """Compute the magnetic vector potential (A field) at a collection of
        target points

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

    Returns:
        (T-m) (N,3) array of magnetic vector potential (A) vectors at each
            target position
    """

    return _evaluate_fields(
        src_mesh,
        targets,
        _RequestedField.AFIELD,
        jdensity=jdensity,
        magnetization=magnetization,
        settings=settings,
    )

b_field 

b_field(
    src_mesh: Mesh,
    targets: NDArray[float64],
    *,
    jdensity: NDArray[float64] | None = None,
    magnetization: NDArray[float64] | None = None,
    settings: SolverSettings = DEFAULT_SETTINGS,
) -> NDArray[float64]

Compute the magnetic flux density (B field) at a collection of target points, assuming the target points are in free space

Parameters:

Name Type Description Default
src_mesh Mesh

mesh to use as the field source

 required
targets NDArray[float64]

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

 required
jdensity NDArray[float64] | None

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

 None
magnetization NDArray[float64] | None

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

 None
settings SolverSettings

selects the solution settings

 DEFAULT_SETTINGS

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 
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
def b_field(
    src_mesh: Mesh,
    targets: NDArray[float64],
    *,  # Force remaining variables to be passed by name
    jdensity: NDArray[float64] | None = None,
    magnetization: NDArray[float64] | None = None,
    settings: SolverSettings = DEFAULT_SETTINGS,
) -> NDArray[float64]:
    """Compute the magnetic flux density (B field) at a collection of target points,
    assuming the target points are in free space

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

    Returns:
        (T) (N,3) array of magnetic flux density (B) vectors at each target position
    """

    if magnetization is not None:
        warnings.warn(
            "Computing magnetic flux density using a magnetized mesh as the source.\n"
            "This calculation is only valid in free space. Ensure that target points "
            "are in free space and not within a magnetized mesh.",
            UserWarning,
            stacklevel=2,
        )

    return MU0 * h_field(
        src_mesh,
        targets,
        jdensity=jdensity,
        magnetization=magnetization,
        settings=settings,
    )

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 
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
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 
225
226
227
228
229
230
231
232
233
234
235
236
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)

curl 

curl(
    fx: NDArray[float64],
    fy: NDArray[float64],
    fz: NDArray[float64],
    spacing: tuple[float, float, float],
) -> tuple[
    NDArray[float64], NDArray[float64], NDArray[float64]
]

Compute the curl of a vector-value function in 3D space on a uniform grid

Parameters:

Name Type Description Default
fx NDArray[float64]

shape (nx, ny, nz) value of function in x-direction at all points in 3d grid

 required
fy NDArray[float64]

shape (nx, ny, nz) value of function in y-direction at all points in 3d grid

 required
fz NDArray[float64]

shape (nx, ny, nz) value of function in z-direction at all points in 3d grid

 required
spacing tuple[float, float, float]

(dx, dy, dz) step size in each direction

 required

Returns:

Type Description
tuple[NDArray[float64], NDArray[float64], NDArray[float64]]

curl_x, curl_y, curl_z: each shape (nx, ny, nz)
Source code in src/oersted/testing.py 
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
282
283
284
285
286
287
def curl(
    fx: NDArray[float64],
    fy: NDArray[float64],
    fz: NDArray[float64],
    spacing: tuple[float, float, float],
) -> tuple[NDArray[float64], NDArray[float64], NDArray[float64]]:
    """Compute the curl of a vector-value function in 3D space on a uniform grid

    Args:
        fx: shape (nx, ny, nz) value of function in x-direction at all points in 3d grid
        fy: shape (nx, ny, nz) value of function in y-direction at all points in 3d grid
        fz: shape (nx, ny, nz) value of function in z-direction at all points in 3d grid
        spacing: (dx, dy, dz) step size in each direction

    Returns:
        curl_x, curl_y, curl_z: each shape (nx, ny, nz)
    """

    (dx, dy, dz) = spacing

    dfz_dy = np.gradient(fz, dy, axis=1)
    dfy_dz = np.gradient(fy, dz, axis=2)

    dfx_dz = np.gradient(fx, dz, axis=2)
    dfz_dx = np.gradient(fz, dx, axis=0)

    dfy_dx = np.gradient(fy, dx, axis=0)
    dfx_dy = np.gradient(fx, dy, axis=1)

    curl_x = dfz_dy - dfy_dz
    curl_y = dfx_dz - dfz_dx
    curl_z = dfy_dx - dfx_dy

    return curl_x, curl_y, curl_z

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 
239
240
241
242
243
244
245
246
247
248
249
250
251
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],
    settings: SolverSettings,
) -> 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

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
def demag_solve(
    mesh: Mesh,
    material: Material,
    h_external: NDArray[float64],
    settings: SolverSettings,
) -> 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).
    """

    element_integration = settings.integration == "element"
    use_octree = settings.method == "octree"

    return magnetization_solve(
        ascontiguousarray(mesh.nodes),
        ascontiguousarray(mesh.connectivity),
        ascontiguousarray(mesh.centroids),
        material.chi(1.0),
        ascontiguousarray(h_external),
        element_integration,
        settings.n_threads,
        settings.atol,
        settings.max_iterations,
        settings.under_relaxation_factor,
        use_octree,
        settings.theta,
        settings.expansion_order(),
        settings.max_leaf_size,
        settings.batch_size,
        settings.verbose,
    )

h_field 

h_field(
    src_mesh: Mesh,
    targets: NDArray[float64],
    *,
    jdensity: NDArray[float64] | None = None,
    magnetization: NDArray[float64] | None = None,
    settings: SolverSettings = DEFAULT_SETTINGS,
) -> NDArray[float64]

Compute the magnetic field strength (H field) at a collection of target points

Parameters:

Name Type Description Default
src_mesh Mesh

mesh to use as the field source

 required
targets NDArray[float64]

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

 required
jdensity NDArray[float64] | None

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

 None
magnetization NDArray[float64] | None

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

 None
settings SolverSettings

selects the solution settings

 DEFAULT_SETTINGS

Returns:

Type Description
NDArray[float64]

(A/m) (N,3) array of magnetic field strength (H) vectors at each target position
Source code in src/oersted/biotsavart.py 
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
def h_field(
    src_mesh: Mesh,
    targets: NDArray[float64],
    *,  # Force remaining variables to be passed by name
    jdensity: NDArray[float64] | None = None,
    magnetization: NDArray[float64] | None = None,
    settings: SolverSettings = DEFAULT_SETTINGS,
) -> NDArray[float64]:
    """Compute the magnetic field strength (H field) at a collection of target points

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

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

    return _evaluate_fields(
        src_mesh,
        targets,
        _RequestedField.HFIELD,
        jdensity=jdensity,
        magnetization=magnetization,
        settings=settings,
    )

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,
    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

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,
    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 | Path,
    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 | Path

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 
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
def make_helmholtz(
    filename: str | Path, 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)

make_ring 

make_ring(
    mesh_size: float = 0.015, jmag: float = 100000000.0
) -> tuple[Mesh, NDArray[float64]]

Make a mesh of a current-carrying ring for testing

Source code in src/oersted/testing.py 
180
181
182
183
184
185
186
187
188
189
190
191
def make_ring(
    mesh_size: float = 15e-3, jmag: float = 1e8
) -> tuple[Mesh, NDArray[float64]]:
    """Make a mesh of a current-carrying ring for testing"""

    mesh = Mesh.from_step("tests/data/ring.stp", mesh_size)
    phi = np.atan2(mesh.centroids[:, 1], mesh.centroids[:, 0])
    jdensity = np.zeros_like(mesh.centroids)
    jdensity[:, 0] = -jmag * np.sin(phi)
    jdensity[:, 1] = jmag * np.cos(phi)

    return mesh, jdensity

max_verr 

max_verr(
    measurement: NDArray[float64],
    baseline: NDArray[float64],
    eps: float = 1e-06,
) -> float

Compute max relative vector error, with the denominator bounded by eps to prevent divide by zero issues

Source code in src/oersted/testing.py 
38
39
40
41
42
43
44
45
def max_verr(
    measurement: NDArray[float64], baseline: NDArray[float64], eps: float = 1e-6
) -> float:
    """Compute max relative vector error, with the denominator bounded by `eps` to
    prevent divide by zero issues
    """

    return float(np.max(verr(measurement, baseline, eps)))

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)

mean_verr 

mean_verr(
    measurement: NDArray[float64],
    baseline: NDArray[float64],
    eps: float = 1e-06,
) -> float

Compute mean relative vector error, with the denominator bounded by eps to prevent divide by zero issues

Source code in src/oersted/testing.py 
28
29
30
31
32
33
34
35
def mean_verr(
    measurement: NDArray[float64], baseline: NDArray[float64], eps: float = 1e-6
) -> float:
    """Compute mean relative vector error, with the denominator bounded by `eps` to
    prevent divide by zero issues
    """

    return float(np.mean(verr(measurement, baseline, eps)))

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 
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
354
355
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
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 
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
282
283
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
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 
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
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)
    if denominator.any() < RTOL:
        raise ValueError("Both measurement and baseline are near zero.")

    return (2 / n) * np.sum(numerator / denominator)

uniform_3d_grid 

uniform_3d_grid(
    xrange: tuple[float, float],
    yrange: tuple[float, float],
    zrange: tuple[float, float],
    n: tuple[int, int, int],
) -> tuple[NDArray[float64], tuple[float, float, float]]

Create a uniform 3D grid

Args xrange: (xmin, xmax) yrange: (ymin, ymax) zrange: (zmin, zmax) n: (nx, ny, nz)

Returns:

Name Type Description
pts NDArray[float64]

shape (nxnynz,3), locations of points in 3d space
(dx, dy, dz)

length increment along each axis
Source code in src/oersted/testing.py 
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
def uniform_3d_grid(
    xrange: tuple[float, float],
    yrange: tuple[float, float],
    zrange: tuple[float, float],
    n: tuple[int, int, int],
) -> tuple[NDArray[float64], tuple[float, float, float]]:
    """Create a uniform 3D grid

    Args
        xrange: (xmin, xmax)
        yrange: (ymin, ymax)
        zrange: (zmin, zmax)
        n: (nx, ny, nz)

    Returns:
        pts: shape (nx*ny*nz,3), locations of points in 3d space
        (dx, dy, dz): length increment along each axis
    """

    (xmin, xmax) = xrange
    (ymin, ymax) = yrange
    (zmin, zmax) = zrange
    (nx, ny, nz) = n

    x = np.linspace(xmin, xmax, nx)
    y = np.linspace(ymin, ymax, ny)
    z = np.linspace(zmin, zmax, nz)
    X, Y, Z = np.meshgrid(x, y, z, indexing="ij")
    pts = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()])

    dx, dy, dz = x[1] - x[0], y[1] - y[0], z[1] - z[0]

    return pts, (dx, dy, dz)

verr 

verr(
    measurement: NDArray[float64],
    baseline: NDArray[float64],
    eps: float = 1e-06,
) -> NDArray[float64]

Compute simple relative vector error, with the denominator bounded by eps to prevent divide by zero issues

Returns:

Type Description
NDArray[float64]

element-wise error at each output point
Source code in src/oersted/testing.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
def verr(
    measurement: NDArray[float64], baseline: NDArray[float64], eps: float = 1e-6
) -> NDArray[float64]:
    """Compute simple relative vector error, with the denominator bounded by `eps` to
    prevent divide by zero issues

    Returns:
        element-wise error at each output point
    """

    return np.linalg.norm(measurement - baseline, axis=1) / (
        np.linalg.norm(baseline, axis=1) + eps
    )