Skip to content

Solenoid Stress & Strain

This package includes three complementary layers:

  • a 1D finite-difference radial stress solver for winding-pack models with zero rz shear,
  • a 2D quadrilateral FEM solver with axisymmetric and plane-strain formulations, matrix-free load assembly, explicit sparse operator exports, and cached Rust-side sparse-LU solves,
  • analytic reference formulas used for validation and convergence studies.

1D Finite-Difference Solver

cfsem.solenoid_stress.SolenoidStress1D

Bases: NumpyModel

Source code in cfsem/solenoid_stress/solenoid_1d.py 
 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
class SolenoidStress1D(NumpyModel):
    model_config = ConfigDict(validate_assignment=True, frozen=True, extra="forbid")

    rgrid: NpNDArray
    """[m] 1D grid of r-coordinates"""
    elasticity_modulus: float
    """[Pa] diagonal terms in material property matrix"""
    poisson_ratio: float
    """[dimensionless] factor determining off-diagonal terms in material property matrix"""
    order: Literal[2, 4] = 4
    """Finite-difference stencil polynomial order.
       Higher order operators produce excessive numerical error under typical use."""
    direct_inverse: bool = False
    """Whether to generate fully-dense direct inverse of the system, which
    can be useful as a linear operator. Alternatively, the system can be solved
    using an LU solver with reduced memory usage and better numerical conditioning."""

    @cached_property
    def operators(self) -> SolenoidStress1DOperators:
        """
        Linear operators for solving stress and strain in a pancake coil
        following Iwasa 2e section 3.6.
        """
        return solenoid_1d_structural_operators(
            np.array(self.rgrid), self.elasticity_modulus, self.poisson_ratio, self.order, self.direct_inverse
        )

    @cached_property
    def displacement_solver(self) -> Callable[[NDArray], NDArray]:
        """LU solver for load-displacement relation (A_ub)
        as an alternative to taking a direct inverse of A_bu"""
        return factorized(self.operators.a_bu)

direct_inverse class-attribute instance-attribute 

direct_inverse: bool = False

Whether to generate fully-dense direct inverse of the system, which can be useful as a linear operator. Alternatively, the system can be solved using an LU solver with reduced memory usage and better numerical conditioning.

displacement_solver cached property 

displacement_solver: Callable[[NDArray], NDArray]

LU solver for load-displacement relation (A_ub) as an alternative to taking a direct inverse of A_bu

elasticity_modulus instance-attribute 

elasticity_modulus: float

[Pa] diagonal terms in material property matrix

operators cached property 

operators: SolenoidStress1DOperators

Linear operators for solving stress and strain in a pancake coil following Iwasa 2e section 3.6.

order class-attribute instance-attribute 

order: Literal[2, 4] = 4

Finite-difference stencil polynomial order. Higher order operators produce excessive numerical error under typical use.

poisson_ratio instance-attribute 

poisson_ratio: float

[dimensionless] factor determining off-diagonal terms in material property matrix

rgrid instance-attribute 

rgrid: NpNDArray

[m] 1D grid of r-coordinates

cfsem.solenoid_stress.SolenoidStress1DOperators dataclass

Linear operators for solving stress and strain in a pancake coil following Iwasa 2e section 3.6.

A_bu, (n x n) sparse operator mapping displacement to the RHS like A @ u_r = -c * j * bz A_ub, (n x n) fully-dense direct inverse of A_bu mapping RHS to displacement A_eu (2n x n), A_eu_radial (n x n), A_eu_hoop (n x n), sparse operators mapping displacement to strain * First entry is combined operator producing both strain components * Second and third entries are split operators, which are equivalent because they are fully decoupled A_se (2n x 2n), sparse operator mapping strain to stress

Source code in cfsem/solenoid_stress/solenoid_1d.py 
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
@dataclass(frozen=True)
class SolenoidStress1DOperators:
    """
    Linear operators for solving stress and strain in a pancake coil
    following Iwasa 2e section 3.6.

    A_bu, (n x n) sparse operator mapping displacement to the RHS like A @ u_r = -c * j * bz
    A_ub, (n x n) fully-dense direct inverse of A_bu mapping RHS to displacement
    A_eu (2n x n), A_eu_radial (n x n), A_eu_hoop (n x n), sparse operators mapping displacement to strain
        * First entry is combined operator producing both strain components
        * Second and third entries are split operators, which are equivalent because they are fully decoupled
    A_se (2n x 2n), sparse operator mapping strain to stress
    """

    a_bu: CSC
    """(n x n) sparse operator mapping displacement to the RHS like A @ u_r = -c * j * bz"""
    a_ub: NDArray | None
    """(n x n) fully-dense direct inverse of A_bu mapping RHS to displacement.
        Only generated if `direct_inverse` flag is set."""
    a_eu: CSR
    """(2n x n), sparse operator mapping displacement to strain; contains both radial and hoop components"""
    a_eu_radial: CSR
    """(n x n), sparse operators mapping displacement to strain; radial component only"""
    a_eu_hoop: CSR
    """(n x n), sparse operators mapping displacement to strain; hoop component only"""
    a_se: CSR
    """(2n x 2n), sparse operator mapping strain to stress"""

    def write_mat(self, dst: str | Path) -> str:
        """Write the collection of operators in .mat format.

        Args:
            dst: Target directory to place the file named "stress_operators.mat"

        Raises:
            IOError: If the directory does not exist
        """
        # Check directory
        dst = Path(dst).absolute()
        fpath = dst / "stress_operators.mat"
        getLogger("cfsem").info(f"Saving stress operator data to {fpath}")

        to_save = {
            "A_bu": self.a_bu,
            "A_ub": self.a_ub,
            "A_eu": self.a_eu,
            "A_eu_radial": self.a_eu_radial,
            "A_eu_hoop": self.a_eu_hoop,
            "A_se": self.a_se,
        }

        if self.a_ub is None:  # savemat fails on None value
            to_save.pop("A_ub")

        #    Note this will implicitly convert all CSR matrices to CSC, which is .mat's preferred I/O
        io.savemat(fpath, to_save)

        return f"{fpath}"

a_bu instance-attribute 

a_bu: csc_matrix

(n x n) sparse operator mapping displacement to the RHS like A @ u_r = -c * j * bz

a_eu instance-attribute 

a_eu: csr_matrix

(2n x n), sparse operator mapping displacement to strain; contains both radial and hoop components

a_eu_hoop instance-attribute 

a_eu_hoop: csr_matrix

(n x n), sparse operators mapping displacement to strain; hoop component only

a_eu_radial instance-attribute 

a_eu_radial: csr_matrix

(n x n), sparse operators mapping displacement to strain; radial component only

a_se instance-attribute 

a_se: csr_matrix

(2n x 2n), sparse operator mapping strain to stress

a_ub instance-attribute 

a_ub: NDArray | None

(n x n) fully-dense direct inverse of A_bu mapping RHS to displacement. Only generated if direct_inverse flag is set.

write_mat 

write_mat(dst: str | Path) -> str

Write the collection of operators in .mat format.

Parameters:

Name Type Description Default
dst str | Path

Target directory to place the file named "stress_operators.mat"

 required

Raises:

Type Description
IOError

If the directory does not exist
Source code in cfsem/solenoid_stress/solenoid_1d.py 
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
def write_mat(self, dst: str | Path) -> str:
    """Write the collection of operators in .mat format.

    Args:
        dst: Target directory to place the file named "stress_operators.mat"

    Raises:
        IOError: If the directory does not exist
    """
    # Check directory
    dst = Path(dst).absolute()
    fpath = dst / "stress_operators.mat"
    getLogger("cfsem").info(f"Saving stress operator data to {fpath}")

    to_save = {
        "A_bu": self.a_bu,
        "A_ub": self.a_ub,
        "A_eu": self.a_eu,
        "A_eu_radial": self.a_eu_radial,
        "A_eu_hoop": self.a_eu_hoop,
        "A_se": self.a_se,
    }

    if self.a_ub is None:  # savemat fails on None value
        to_save.pop("A_ub")

    #    Note this will implicitly convert all CSR matrices to CSC, which is .mat's preferred I/O
    io.savemat(fpath, to_save)

    return f"{fpath}"

cfsem.solenoid_stress.solenoid_1d_structural_factor 

solenoid_1d_structural_factor(
    elasticity_modulus: float, poisson_ratio: float
) -> float

Structural factor applied to RHS of solenoid stress solve

Source code in cfsem/solenoid_stress/solenoid_1d.py 
42
43
44
45
def solenoid_1d_structural_factor(elasticity_modulus: float, poisson_ratio: float) -> float:
    """Structural factor applied to RHS of solenoid stress solve"""
    c = (1.0 - poisson_ratio**2) / elasticity_modulus  # [m/N]
    return c

cfsem.solenoid_stress.solenoid_1d_structural_rhs 

solenoid_1d_structural_rhs(
    c: float,
    j: NDArray | list[float],
    bz: NDArray | list[float],
    pi: float = 0.0,
    po: float = 0.0,
) -> NDArray

Right-hand-side for solenoid stress solve, including zero values at the BCs.

From Iwasa 2e eqn. 3.64a

Recommend padding the grid with a dummy value at either end to make room for the BCs without losing accounting of nonzero current density at the inner/outer radius.

Padding for BCs can be done like: rgrid = np.array([r0 - 1e-6] + rgrid.tolist() + [r1 + 1e-6])

Padding region is ultimately treated as structural material, so the padded region should be small to avoid introducing error, but not so small that it causes numerical error in the finite difference scheme.

Parameters:

Name Type Description Default
c float

[m/N] scalar structural factor; see solenoid_1d_structural_factor()

 required
j NDArray | list[float]

[A/m^2] with shape (n x 1), current density at each point in the r-grid

 required
bz NDArray | list[float]

[T] with shape (n x 1), Z-axis B-field at each point in the r-grid

 required
pi float

[Pa] scalar pressure on inner wall, defined in +r direction

 0.0
po float

[Pa] scalar pressure on outer wall, defined in -r direction

 0.0

Returns:

Type Description
NDArray

-c * j * bz, [1/m^2] with shape (n x 1), the right-hand side of the solenoid stress PDE
Source code in cfsem/solenoid_stress/solenoid_1d.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
def solenoid_1d_structural_rhs(
    c: float,
    j: NDArray | list[float],
    bz: NDArray | list[float],
    pi: float = 0.0,
    po: float = 0.0,
) -> NDArray:
    """
    Right-hand-side for solenoid stress solve,
    including zero values at the BCs.

    From Iwasa 2e eqn. 3.64a

    Recommend padding the grid with a dummy value at either end
    to make room for the BCs without losing accounting of nonzero
    current density at the inner/outer radius.

    Padding for BCs can be done like:
    `rgrid = np.array([r0 - 1e-6] + rgrid.tolist() + [r1 + 1e-6])`

    Padding region is ultimately treated as structural material,
    so the padded region should be small to avoid introducing error,
    but not so small that it causes numerical error in the finite difference scheme.

    Args:
        c: [m/N] scalar structural factor; see `solenoid_1d_structural_factor()`
        j: [A/m^2] with shape (n x 1), current density at each point in the r-grid
        bz: [T] with shape (n x 1), Z-axis B-field at each point in the r-grid
        pi: [Pa] scalar pressure on inner wall, defined in +r direction
        po: [Pa] scalar pressure on outer wall, defined in -r direction

    Returns:
        -c * j * bz, [1/m^2] with shape (n x 1), the right-hand side of the solenoid stress PDE
    """
    # Guarantee arrays
    j = np.array(j)
    bz = np.array(bz)

    # RHS without BCs
    rhs = -c * j * bz

    # BC for r-stress at inner and outer radius
    # is the fluid or mechanical pressure, which is usually going to be set to zero
    # to represent an unsupported system, but could be set to a nonzero value
    # to represent a surface load.
    rhs[0] = -pi  # Sign convention: compression is negative stress
    rhs[-1] = -po

    return rhs

2D FEM

The FEM path supports:

  • axisymmetric and plane-strain structural formulations,
  • quad4, inferred quad9, and explicit quad9 elements,
  • gl3 and gl4 quadrature,
  • optional per-element in-plane material orientation angles,
  • optional threaded stiffness assembly with par=True,
  • explicit reduced-space operator exports for body force, pressure, traction, and nodal-temperature thermal strain,
  • explicit location-based sparse operator exports for interpolation, strain, and stress,
  • matrix-free location-based strain, stress, thermal-strain, and thermal-stress recovery,
  • direct sparse-LU reduced-system solves,
  • float64 numeric storage; floating input arrays must already have dtype float64,
  • model-owned Dirichlet constraints applied during assembly.

The intended workflow is:

  1. call assemble_structural_2d(...) once with mesh, materials, load topology, and prescribed Dirichlet values,
  2. build each reduced load vector with matrix-free model.build_rhs(...) or user-owned sparse operator exports,
  3. solve with model.solve(rhs), using the cached sparse-LU factorization,
  4. recover fields with model.strain(locations, displacement), model.stress(locations, displacement), model.thermal_strain(locations, temperature), or model.thermal_stress(locations, temperature), using locations from model.quadrature().locations, model.locate_points(...), or model.locate_points_in_elements(...).

By default, model.solve(rhs) uses the direct sparse-LU path and returns the full displacement array. The sparse-LU factorization is built lazily on the first solve and then cached on the model for repeated right-hand sides.

Location-based recovery returns flat point-major arrays. model.quadrature() returns Quadrature; pass quadrature.locations to recovery methods, and use quadrature.weights_area, quadrature.weights_volume, and quadrature.points_per_element for integrating quantities over elements. model.locate_points(...) performs a mesh query for arbitrary physical points, while model.locate_points_in_elements(...) is the cheaper path when element ownership is already known. Existing QuadMeshQuery results can be converted with query.point_locations() and passed to the same recovery methods. Sparse recovery exports use the same locations: model.interpolation_operator(locations), model.strain_operator(locations), and model.stress_operator(locations).

Formulation Notes

The axisymmetric and plane-strain solvers share the same 2D quadrilateral mesh, two displacement unknowns per node, and four-component strain/stress storage. The difference is how that 2D mesh represents a 3D body.

For formulation="axisymmetric", the coordinates are interpreted as (r, z). Each quadrature area sample represents a full ring, so stiffness and load integrals use the volume scale 2*pi*r*dA. The strain vector is [rr, zz, tt, rz]; the out-of-plane hoop strain is not an independent displacement derivative, but is recovered from the radial displacement as epsilon_tt = u_r / r.

For formulation="plane_strain", the coordinates are interpreted as (x, y). Each area sample represents a prismatic slice with user-supplied thickness, so integrals use the volume scale thickness*dA. The strain vector is [xx, yy, zz, xy]; the out-of-plane strain is constrained to epsilon_zz = 0, while sigma_zz can still be nonzero through the constitutive matrix.

Plane strain and plane stress are different 2D reductions. Plane strain models a body that is long, periodic, or otherwise constrained in the out-of-plane direction, with zero out-of-plane strain and generally nonzero out-of-plane stress. Plane stress models a thin sheet or plate with traction-free faces through the thickness, with zero out-of-plane stress and generally nonzero out-of-plane strain. This FEM path currently implements axisymmetric and plane-strain reductions; it does not implement a plane-stress constitutive reduction.

cfsem.solenoid_stress.fem2d

2D structural elasticity finite-element assembly.

This module provides a small displacement-based quadrilateral FEM solver for axisymmetric and plane-strain structural reductions. The backend stores the reduced stiffness matrix and evaluates loads and quadrature recovery matrix-free unless sparse operators are explicitly exported.

The element formulation follows the standard small-strain Galerkin construction

K_e = integral(B^T D B c dA) where c is 2*pi*r for axisymmetric and the thickness of the planar domain for plane strain.

with consistent body-force, surface-pressure, and surface-traction load vectors. The axisymmetric engineering-strain vector is ordered as [e_rr, e_zz, e_tt, g_rz]. In Bower's terminology, the underlying equations are the strain-displacement equation, the elastic stress-strain law, the equation of static equilibrium for stresses, and the boundary conditions on displacement and stress.

References

[1] Allan F. Bower, Applied Mechanics of Solids, CRC Press, 2009. See especially Section 8.1 and Table 8.3 for the general displacement-based finite-element construction and 2D interpolation functions.

[2] E. L. Wilson, "Structural Analysis of Axisymmetric Solids," AIAA Journal, 3(12), pp. 2269-2274, 1965.

[3] R. A. Mitchell, R. M. Woolley, and C. R. Fisher, "Formulation and experimental verification of an axisymmetric finite-element structural analysis," Journal of Research of the National Bureau of Standards Section C, 75C, 1971.

[4] I. Fried, "Notes on the finite element analysis of the axisymmetric elastic solid," International Journal of Solids and Structures, 10(3), 1974.

ElementMeasures dataclass

Per-element cross-section area and represented volume.

areas and volumes both have shape (nelem,). areas has units [area] and volumes has units [volume].

Source code in cfsem/solenoid_stress/fem2d.py 
125
126
127
128
129
130
131
132
133
134
@dataclass(frozen=True, slots=True)
class ElementMeasures:
    """Per-element cross-section area and represented volume.

    `areas` and `volumes` both have shape `(nelem,)`.
    `areas` has units `[area]` and `volumes` has units `[volume]`.
    """

    areas: Float64Array
    volumes: Float64Array

ElevatedQuad9Mesh dataclass

Explicit 9-node analysis mesh inferred from a corner-only quad4 mesh.

analysis_elements use the local quad9 ordering: - corners 0..3 in counter-clockwise order [bottom-left, bottom-right, top-right, top-left] - midsides 4..7 on faces [bottom, right, top, left] - center node 8

input_nodes and analysis_nodes have units [length].

Source code in cfsem/solenoid_stress/fem2d.py 
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
@dataclass(frozen=True, slots=True)
class ElevatedQuad9Mesh:
    """Explicit 9-node analysis mesh inferred from a corner-only quad4 mesh.

    `analysis_elements` use the local quad9 ordering:
    - corners `0..3` in counter-clockwise order `[bottom-left, bottom-right, top-right, top-left]`
    - midsides `4..7` on faces `[bottom, right, top, left]`
    - center node `8`

    `input_nodes` and `analysis_nodes` have units `[length]`.
    """

    input_nodes: Float64Array
    input_elements: UInt64Array
    analysis_nodes: Float64Array
    analysis_elements: UInt64Array
    corner_node_indices: npt.NDArray[np.int64]
    midside_node_indices: npt.NDArray[np.int64]
    center_node_indices: npt.NDArray[np.int64]

PointLocations dataclass

Element-owned physical and reference point locations.

A location is a physical point together with the element that owns or is nearest to that point and the corresponding element-local reference coordinates. Recovery and sparse operator construction use element_indices and reference_points as the source of truth; points is included for caller inspection, plotting, and compatibility with mesh-query outputs.

points has shape (npoint, 2) and units [length]. element_indices has shape (npoint,) and stores unitless analysis-element indices. reference_points has shape (npoint, 2) and stores unitless coordinates in the element's [-1, 1]^2 reference domain. element_type records the element family that produced the locations so model methods can reject locations from an incompatible mesh.

Source code in cfsem/solenoid_stress/fem2d.py 
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
@dataclass(frozen=True, slots=True)
class PointLocations:
    """Element-owned physical and reference point locations.

    A location is a physical point together with the element that owns or is nearest to that point
    and the corresponding element-local reference coordinates. Recovery and sparse operator
    construction use `element_indices` and `reference_points` as the source of truth; `points` is
    included for caller inspection, plotting, and compatibility with mesh-query outputs.

    `points` has shape `(npoint, 2)` and units `[length]`. `element_indices` has shape
    `(npoint,)` and stores unitless analysis-element indices. `reference_points` has shape
    `(npoint, 2)` and stores unitless coordinates in the element's `[-1, 1]^2` reference domain.
    `element_type` records the element family that produced the locations so model methods can
    reject locations from an incompatible mesh.

    """

    points: Float64Array
    element_indices: npt.NDArray[np.int64]
    reference_points: Float64Array
    element_type: str

QuadMeshInterpolation dataclass

Interpolated nodal values and element-location metadata for query points.

values has shape (npoint, ...), where ... is the trailing shape of the nodal values. element_indices stores the nearest element used for interpolation. inside reports whether the nearest-element distance was within the containment tolerance.

Source code in cfsem/solenoid_stress/fem2d.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
@dataclass(frozen=True, slots=True)
class QuadMeshInterpolation:
    """Interpolated nodal values and element-location metadata for query points.

    `values` has shape `(npoint, ...)`, where `...` is the trailing shape of the nodal values.
    `element_indices` stores the nearest element used for interpolation. `inside` reports whether
    the nearest-element distance was within the containment tolerance.
    """

    values: Float64Array
    element_indices: npt.NDArray[np.int64]
    reference_points: Float64Array
    inside: npt.NDArray[np.bool_]

QuadMeshQuery dataclass

One-pass geometric query results for points in a 2D quadrilateral mesh.

The query stores nearest-node, nearest-element, and nearest-face data for each query point. Interpolation and recovery operators can reuse this object without repeating the mesh search. For contained points, the nearest element is the containing element and nearest_element_distances is zero to numerical tolerance.

Source code in cfsem/solenoid_stress/fem2d.py 
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
@dataclass(frozen=True, slots=True)
class QuadMeshQuery:
    """One-pass geometric query results for points in a 2D quadrilateral mesh.

    The query stores nearest-node, nearest-element, and nearest-face data for each query point.
    Interpolation and recovery operators can reuse this object without repeating the mesh search.
    For contained points, the nearest element is the containing element and
    `nearest_element_distances` is zero to numerical tolerance.
    """

    nodes: Float64Array
    elements: UInt64Array
    points: Float64Array
    element_type: str
    nearest_node_indices: npt.NDArray[np.int64]
    nearest_node_points: Float64Array
    nearest_node_distances: Float64Array
    nearest_element_indices: npt.NDArray[np.int64]
    nearest_element_reference_points: Float64Array
    nearest_element_points: Float64Array
    nearest_element_distances: Float64Array
    nearest_face_element_indices: npt.NDArray[np.int64]
    nearest_face_local_faces: npt.NDArray[np.int64]
    nearest_face_reference_coordinates: Float64Array
    nearest_face_points: Float64Array
    nearest_face_distances: Float64Array

    def point_locations(self) -> PointLocations:
        """Return nearest-element locations for recovery and sparse operators.

        The returned locations reuse the element ownership and reference coordinates already found
        by the mesh query. No additional mesh search or point projection is performed. Query points
        outside the mesh are represented by their nearest projected element points; callers that
        need strict containment should check `nearest_element_distances` before using the locations.

        Returns:
            PointLocations: Element-owned nearest-element locations.
        """

        return PointLocations(
            points=self.nearest_element_points,
            element_indices=self.nearest_element_indices,
            reference_points=self.nearest_element_reference_points,
            element_type=self.element_type,
        )
point_locations 
point_locations() -> PointLocations

Return nearest-element locations for recovery and sparse operators.

The returned locations reuse the element ownership and reference coordinates already found by the mesh query. No additional mesh search or point projection is performed. Query points outside the mesh are represented by their nearest projected element points; callers that need strict containment should check nearest_element_distances before using the locations.

Returns:

Name Type Description
PointLocations PointLocations

Element-owned nearest-element locations.
Source code in cfsem/solenoid_stress/fem2d.py 
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
def point_locations(self) -> PointLocations:
    """Return nearest-element locations for recovery and sparse operators.

    The returned locations reuse the element ownership and reference coordinates already found
    by the mesh query. No additional mesh search or point projection is performed. Query points
    outside the mesh are represented by their nearest projected element points; callers that
    need strict containment should check `nearest_element_distances` before using the locations.

    Returns:
        PointLocations: Element-owned nearest-element locations.
    """

    return PointLocations(
        points=self.nearest_element_points,
        element_indices=self.nearest_element_indices,
        reference_points=self.nearest_element_reference_points,
        element_type=self.element_type,
    )

Quadrature dataclass

Element-major quadrature locations and mapped integration weights.

locations stores the physical points, owning elements, and reference coordinates used by recovery and sparse operator methods. weights_area and weights_volume have shape (npoint,); reshape them as (nelem, points_per_element) for integrating quantities over elements.

Source code in cfsem/solenoid_stress/fem2d.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
@dataclass(frozen=True, slots=True)
class Quadrature:
    """Element-major quadrature locations and mapped integration weights.

    `locations` stores the physical points, owning elements, and reference coordinates used by
    recovery and sparse operator methods. `weights_area` and `weights_volume` have shape
    `(npoint,)`; reshape them as `(nelem, points_per_element)` for integrating quantities over
    elements.
    """

    locations: PointLocations
    weights_area: Float64Array
    weights_volume: Float64Array
    points_per_element: int

Structural2DFEMModel

Reusable 2D structural FEM model with sparse operators and reduced solve state.

Structural FEM numeric arrays are float64; floating inputs must already use float64 arrays.

Load and stiffness operators are exported from the Rust backend on demand: - body_force_to_rhs, pressure_to_rhs, traction_to_rhs, and temperature_to_rhs map load amplitudes to the reduced structural right-hand side, Field recovery uses explicit PointLocations objects. Use quadrature().locations for quadrature locations, locate_points(...) for arbitrary physical points, locate_points_in_elements(...) when element ownership is already known, or QuadMeshQuery.point_locations() to reuse an existing mesh query.

Matrix-free field methods (strain, stress, thermal_strain, thermal_stress) evaluate values directly at supplied locations. Sparse operator methods (interpolation_operator, strain_operator, stress_operator) materialize user-owned SciPy matrices for workflows that apply the same located recovery many times.

Key public array shapes and units: - stiffness has shape (ndof_reduced, ndof_reduced) with entry units [generalized force / displacement] = [energy / distance^2], - body_force_to_rhs has shape (ndof_reduced, 2 * nelem) with entry units [volume], - pressure_to_rhs has shape (ndof_reduced, n_pressure_faces) with entry units [area], - traction_to_rhs has shape (ndof_reduced, 2 * n_traction_faces) with entry units [area], - temperature_to_rhs has shape (ndof_reduced, n_temperature_nodes) with entry units [generalized force / temperature] = [energy / (distance * temperature)].

input_nodes and analysis_nodes have shape (nnode, 2) and units [length]. input_elements and analysis_elements expose the original and analysis connectivity.

Source code in cfsem/solenoid_stress/fem2d.py 
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
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
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
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
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
class Structural2DFEMModel:
    """Reusable 2D structural FEM model with sparse operators and reduced solve state.

    Structural FEM numeric arrays are `float64`; floating inputs must already use `float64` arrays.

    Load and stiffness operators are exported from the Rust backend on demand:
    - `body_force_to_rhs`, `pressure_to_rhs`, `traction_to_rhs`, and `temperature_to_rhs`
      map load amplitudes to the reduced structural right-hand side,
    Field recovery uses explicit `PointLocations` objects. Use `quadrature().locations` for quadrature
    locations, `locate_points(...)` for arbitrary physical points,
    `locate_points_in_elements(...)` when element ownership is already known, or
    `QuadMeshQuery.point_locations()` to reuse an existing mesh query.

    Matrix-free field methods (`strain`, `stress`, `thermal_strain`, `thermal_stress`) evaluate
    values directly at supplied locations. Sparse operator methods (`interpolation_operator`,
    `strain_operator`, `stress_operator`) materialize user-owned SciPy matrices for workflows that
    apply the same located recovery many times.

    Key public array shapes and units:
    - `stiffness` has shape `(ndof_reduced, ndof_reduced)` with entry units
      `[generalized force / displacement] = [energy / distance^2]`,
    - `body_force_to_rhs` has shape `(ndof_reduced, 2 * nelem)` with entry units `[volume]`,
    - `pressure_to_rhs` has shape `(ndof_reduced, n_pressure_faces)` with entry units `[area]`,
    - `traction_to_rhs` has shape `(ndof_reduced, 2 * n_traction_faces)` with entry units
      `[area]`,
    - `temperature_to_rhs` has shape `(ndof_reduced, n_temperature_nodes)` with entry units
      `[generalized force / temperature] = [energy / (distance * temperature)]`.

    `input_nodes` and `analysis_nodes` have shape `(nnode, 2)` and units `[length]`.
    `input_elements` and `analysis_elements` expose the original and analysis connectivity.
    """

    def __init__(
        self,
        *,
        backend: Any,
        input_nodes: Float64Array,
        input_elements: UInt64Array,
        analysis_nodes: Float64Array,
        analysis_elements: UInt64Array,
        elevated: ElevatedQuad9Mesh | None,
        material_ids: UInt64Array,
        material_table: Float64Array,
        thermal_material_table: Float64Array | None,
        material_orientation_angles: Float64Array,
        pressure_faces: UInt64Array,
        traction_faces: UInt64Array,
        formulation: str,
        thickness: float,
        element_type: str,
        free_dofs: npt.NDArray[np.int64],
        fixed_dofs: npt.NDArray[np.int64],
        fixed_values: Float64Array,
        ndof_full: int,
        ndof_reduced: int,
        nelem: int,
        nq_per_element: int,
        n_temperature_nodes: int,
    ) -> None:
        self._backend = backend
        self._input_nodes = input_nodes
        self._input_elements = input_elements
        self._elevated = elevated
        self._material_ids = material_ids
        self._material_table = material_table
        self._thermal_material_table = thermal_material_table
        self._material_orientation_angles = material_orientation_angles
        self.pressure_faces = pressure_faces
        self.traction_faces = traction_faces
        self.analysis_nodes = analysis_nodes
        self.analysis_elements = analysis_elements
        self.free_dofs = free_dofs
        self.fixed_dofs = fixed_dofs
        self.fixed_values = fixed_values
        self.formulation = formulation
        self.thickness = thickness
        self.element_type = element_type
        if formulation == "axisymmetric":
            self.coordinate_labels = ("r", "z")
            self.displacement_labels = ("u_r", "u_z")
            self.tensor_labels = ("rr", "zz", "tt", "rz")
            self.measure_label = "swept_volume"
        else:
            self.coordinate_labels = ("x", "y")
            self.displacement_labels = ("u_x", "u_y")
            self.tensor_labels = ("xx", "yy", "zz", "xy")
            self.measure_label = "volume"
        self.ndof_full = int(ndof_full)
        self.ndof_reduced = int(ndof_reduced)
        self.nelem = int(nelem)
        self.nq_per_element = int(nq_per_element)
        self.n_temperature_nodes = int(n_temperature_nodes)
        self._quadrature_cache: Quadrature | None = None
        self._element_measures_cache: ElementMeasures | None = None

    @property
    def input_nodes(self) -> Float64Array:
        """Corner-node input mesh coordinates with shape `(nnode, 2)` and units `[length]`."""

        return self._input_nodes

    @property
    def input_elements(self) -> npt.NDArray[np.uint64]:
        """Input mesh connectivity with shape `(nelem, 4)`."""

        return self._input_elements

    @cached_property
    def constant_rhs(self) -> Float64Array:
        """Load-independent reduced RHS contribution, exported from Rust on first access."""

        return np.asarray(self._backend.constant_rhs(), dtype=np.float64)

    def quadrature(self) -> Quadrature:
        """Return element-major quadrature locations and mapped integration weights.

        The returned locations are built directly from the model's quadrature rule and element
        geometry, so no global mesh query or point inversion is performed. Results are cached
        because the model mesh, quadrature rule, and geometry are immutable after assembly.

        Returns:
            Quadrature: Flat element-major locations with `nelem * nq_per_element` rows plus
            mapped area and volume integration weights. Pass `quadrature.locations` to recovery
            and sparse operator methods.
        """

        cache = self._quadrature_cache
        if cache is None:
            points, element_indices, reference_points, weights_area, weights_volume, points_per_element = (
                self._backend.quadrature()
            )
            locations = PointLocations(
                points=np.asarray(points, dtype=np.float64).reshape(-1, 2),
                element_indices=np.asarray(element_indices, dtype=np.int64),
                reference_points=np.asarray(reference_points, dtype=np.float64).reshape(-1, 2),
                element_type=self.element_type,
            )
            cache = Quadrature(
                locations=locations,
                weights_area=np.asarray(weights_area, dtype=np.float64),
                weights_volume=np.asarray(weights_volume, dtype=np.float64),
                points_per_element=int(points_per_element),
            )
            self._quadrature_cache = cache
        return cache

    def locate_points(
        self,
        points: ArrayLike,
        *,
        outside: str = "nearest",
        tolerance: float | None = None,
        max_iterations: int = 20,
    ) -> PointLocations:
        """Locate arbitrary physical points in this model's analysis mesh.

        This uses the current brute-force quadrilateral mesh query and returns the nearest element
        plus reference coordinates for each query point. Points outside the mesh are projected to
        the nearest element unless `outside="raise"` or `outside="error"` is supplied.

        Args:
            points: Physical coordinates with shape `(npoint, 2)` and units `[length]`.
            outside: Outside-mesh policy. `"nearest"` returns nearest-element projections;
                `"raise"` and `"error"` raise if any point is outside `tolerance`.
            tolerance: Nonnegative physical distance used to classify contained points. Defaults
                to `1e-10`.
            max_iterations: Maximum local inverse-map iterations per element during the query.

        Returns:
            PointLocations: Located points for recovery and sparse operator construction.

        Raises:
            ValueError: If `outside` requests an error and at least one point is outside the mesh.
        """

        query = query_quad_mesh(
            self.analysis_nodes,
            self.analysis_elements,
            points,
            element_type=self.element_type,
            max_iterations=max_iterations,
        )
        outside_policy = str(outside).strip().lower()
        assert outside_policy in {
            "nearest",
            "raise",
            "error",
        }, f"unsupported outside policy {outside!r}; use 'nearest' or 'raise'"
        tol = _normalize_query_tolerance(tolerance)
        inside = query.nearest_element_distances <= tol
        if outside_policy in {"raise", "error"} and not np.all(inside):
            first = int(np.flatnonzero(~inside)[0])
            raise ValueError(f"query point {first} is outside the quad mesh")
        return PointLocations(
            points=query.nearest_element_points,
            element_indices=query.nearest_element_indices,
            reference_points=query.nearest_element_reference_points,
            element_type=self.element_type,
        )

    def locate_points_in_elements(
        self,
        points: ArrayLike,
        element_indices: ArrayLike,
        *,
        max_iterations: int = 20,
    ) -> PointLocations:
        """Project physical points into caller-supplied owning elements.

        This is the fast path when the caller already knows element ownership, such as when
        reusing element indices returned from `quadrature()` or a previous mesh query. It performs
        one local element projection per point and does not scan the global mesh.

        Args:
            points: Physical coordinates with shape `(npoint, 2)` and units `[length]`.
            element_indices: Owning analysis-element indices with shape `(npoint,)`.
            max_iterations: Maximum inverse-map iterations for each local element projection.

        Returns:
            PointLocations: Projected physical points, caller-supplied element indices, and
            reference coordinates.
        """

        points_arr = _normalize_query_points(points)
        element_indices_arr = np.asarray(element_indices, dtype=np.uint64).reshape(-1)
        projected_points, projected_elements, reference_points = self._backend.locate_points_in_elements(
            points_arr,
            element_indices_arr,
            int(max_iterations),
        )
        return PointLocations(
            points=np.asarray(projected_points, dtype=np.float64).reshape(-1, 2),
            element_indices=np.asarray(projected_elements, dtype=np.int64),
            reference_points=np.asarray(reference_points, dtype=np.float64).reshape(-1, 2),
            element_type=self.element_type,
        )

    @cached_property
    def _temperature_elevation(self) -> sp.csr_matrix:
        """Return the cached input-to-analysis temperature elevation operator.

        Corner-node quad9 inputs are elevated inside the Rust backend before assembly, while the
        Python API still accepts temperatures on the original input nodes.  This operator bridges
        those two spaces for exported scipy operators.  It is cached because the input mesh and
        inferred analysis mesh are immutable for the lifetime of the model.
        """

        elevated = self._elevated
        assert elevated is not None, "temperature elevation is available only for inferred quad9 meshes"
        return _temperature_elevation_operator(elevated)

    @cached_property
    def stiffness(self) -> sp.csc_matrix:
        """Reduced stiffness matrix with shape `(ndof_reduced, ndof_reduced)`.

        Entries have units `[generalized force / displacement] = [energy / distance^2]`.
        The SciPy matrix is exported from the Rust backend on first access and then cached.
        """

        return _csc_matrix_from_binding(self._backend.stiffness_csc())

    @property
    def body_force_to_rhs(self) -> sp.csr_matrix:
        """Operator mapping per-element body-force density to the reduced RHS.

        Shape is `(ndof_reduced, 2 * nelem)`. Entries have units `[volume]`.
        """

        return _csr_matrix_from_binding(self._backend.body_force_to_rhs_csr())

    @property
    def pressure_to_rhs(self) -> sp.csr_matrix:
        """Operator mapping scalar pressure amplitudes to the reduced RHS.

        Shape is `(ndof_reduced, n_pressure_faces)`. Entries have units `[area]`.
        """

        return _csr_matrix_from_binding(self._backend.pressure_to_rhs_csr())

    @property
    def traction_to_rhs(self) -> sp.csr_matrix:
        """Operator mapping vector traction amplitudes to the reduced RHS.

        Shape is `(ndof_reduced, 2 * n_traction_faces)`. Entries have units `[area]`.
        """

        return _csr_matrix_from_binding(self._backend.traction_to_rhs_csr())

    @property
    def temperature_to_rhs(self) -> sp.csr_matrix:
        """Operator mapping input-node temperatures to the reduced RHS.

        Shape is `(ndof_reduced, n_temperature_nodes)`. Entries have units
        `[generalized force / temperature] = [energy / (distance * temperature)]`.
        """

        analysis_operator = _csr_matrix_from_binding(self._backend.temperature_to_rhs_csr())
        return (
            sp.csr_matrix(analysis_operator @ self._temperature_elevation)
            if self._elevated is not None and self.n_temperature_nodes > 0
            else analysis_operator
        )

    def interpolation_operator(self, locations: PointLocations) -> sp.csr_matrix:
        """Build a sparse interpolation operator for located points.

        Args:
            locations: Element-owned point locations from `quadrature().locations`,
                `locate_points(...)`, `locate_points_in_elements(...)`, or
                `QuadMeshQuery.point_locations()`.

        Returns:
            csr_matrix: Sparse operator with shape `(npoint, n_analysis_nodes)`. Multiplying by a
            scalar nodal field with shape `(n_analysis_nodes,)` returns interpolated values with
            shape `(npoint,)`; multiplying by `(n_analysis_nodes, ncomponent)` interpolates each
            component independently. Entries are unitless shape-function values.
        """

        locations = self._validate_locations(locations)
        return _coo_operator_from_binding(
            _quad_mesh_interpolation_operator_f64(
                self.analysis_nodes,
                self.analysis_elements,
                locations.element_indices.astype(np.uint64, copy=False),
                locations.reference_points,
                self.element_type,
            ),
        )

    def strain_operator(self, locations: PointLocations) -> sp.csr_matrix:
        """Build a sparse total-strain recovery operator for located points.

        Args:
            locations: Element-owned point locations from `quadrature().locations`,
                `locate_points(...)`, `locate_points_in_elements(...)`, or
                `QuadMeshQuery.point_locations()`.

        Returns:
            csr_matrix: Sparse operator with shape `(4 * npoint, 2 * n_analysis_nodes)`. Rows are
            grouped by point and tensor component. Multiplying by full analysis displacements with
            shape `(2 * n_analysis_nodes,)` returns flat strain samples with shape `(4 * npoint,)`.
            Entries have units `[1 / length]`.
        """

        locations = self._validate_locations(locations)
        return _coo_operator_from_binding(
            _quad_mesh_strain_operator_f64(
                self.analysis_nodes,
                self.analysis_elements,
                locations.element_indices.astype(np.uint64, copy=False),
                locations.reference_points,
                self.element_type,
                _formulation_code(self.formulation),
                self.thickness,
            ),
        )

    def stress_operator(self, locations: PointLocations) -> sp.csr_matrix:
        """Build a sparse elastic-stress recovery operator for located points.

        Args:
            locations: Element-owned point locations from `quadrature().locations`,
                `locate_points(...)`, `locate_points_in_elements(...)`, or
                `QuadMeshQuery.point_locations()`.

        Returns:
            csr_matrix: Sparse operator with shape `(4 * npoint, 2 * n_analysis_nodes)`. Rows are
            grouped by point and tensor component. Multiplying by full analysis displacements with
            shape `(2 * n_analysis_nodes,)` returns flat stress samples with shape `(4 * npoint,)`.
            Entries have units `[stress / length]`.
        """

        locations = self._validate_locations(locations)
        return _coo_operator_from_binding(
            _quad_mesh_stress_operator_f64(
                self.analysis_nodes,
                self.analysis_elements,
                locations.element_indices.astype(np.uint64, copy=False),
                locations.reference_points,
                self._material_ids,
                self._material_table,
                self._material_orientation_angles,
                self.element_type,
                _formulation_code(self.formulation),
                self.thickness,
            ),
        )

    def element_measures(self) -> ElementMeasures:
        """Return cross-section area and represented volume for each element.

        Returns:
            ElementMeasures: Per-element measures with:
                `areas` of shape `(nelem,)` and units `[area]`,
                `volumes` of shape `(nelem,)` and units `[volume]`.
        """

        cache = self._element_measures_cache
        if cache is not None:
            return cache
        quadrature = self.quadrature()
        weights_area = quadrature.weights_area.reshape(self.nelem, quadrature.points_per_element)
        weights_volume = quadrature.weights_volume.reshape(self.nelem, quadrature.points_per_element)
        cache = ElementMeasures(
            areas=np.asarray(weights_area.sum(axis=1), dtype=np.float64),
            volumes=np.asarray(weights_volume.sum(axis=1), dtype=np.float64),
        )
        self._element_measures_cache = cache
        return cache

    def _normalize_temperature_for_backend(
        self,
        nodal_temperature: ArrayLike | None,
    ) -> npt.NDArray[np.floating[Any]] | None:
        if self.n_temperature_nodes == 0:
            values = (
                np.zeros((0,), dtype=np.float64)
                if nodal_temperature is None
                else np.asarray(nodal_temperature).reshape(-1)
            )
            assert values.size == 0, "nodal_temperature was provided, but this model has no thermal operator"
            return None
        if nodal_temperature is None:
            raise ValueError("nodal_temperature is required because this model includes thermal materials")
        input_temperature = _normalize_nodal_temperature(nodal_temperature, self._input_nodes.shape[0])
        return (
            np.asarray(self._temperature_elevation @ input_temperature, dtype=np.float64)
            if self._elevated is not None
            else input_temperature
        )

    def _validate_locations(self, locations: PointLocations) -> PointLocations:
        assert isinstance(locations, PointLocations), "locations must be a PointLocations object"
        assert (
            locations.element_type == self.element_type
        ), f"locations use element_type {locations.element_type!r}, but model uses {self.element_type!r}"
        assert (
            locations.points.ndim == 2 and locations.points.shape[1] == 2
        ), f"locations.points must have shape (npoint, 2); got {locations.points.shape}"
        assert (
            locations.reference_points.ndim == 2 and locations.reference_points.shape[1] == 2
        ), f"locations.reference_points must have shape (npoint, 2); got {locations.reference_points.shape}"
        assert (
            locations.element_indices.ndim == 1
        ), f"locations.element_indices must have shape (npoint,); got {locations.element_indices.shape}"
        npoint = locations.element_indices.shape[0]
        assert (
            locations.points.shape[0] == npoint
        ), f"locations.points has {locations.points.shape[0]} rows, but element_indices has {npoint}"
        assert locations.reference_points.shape[0] == npoint, (
            "locations.reference_points has "
            f"{locations.reference_points.shape[0]} rows, but element_indices has {npoint}"
        )
        return locations

    def build_rhs(
        self,
        body_force: ArrayLike | None = None,
        pressure_values: ArrayLike | None = None,
        traction_values: ArrayLike | None = None,
        nodal_temperature: ArrayLike | None = None,
    ) -> npt.NDArray[np.floating[Any]]:
        """Build one reduced structural right-hand side.

        Args:
            body_force: Elementwise body-force amplitudes with shape `(2,)` or `(nelem, 2)`.
                Components are `[b_r, b_z]` with units `[force / volume]`.
            pressure_values: Pressure amplitudes with shape `(n_pressure_faces,)` and units
                `[force / area]`. Positive values act in the inward normal direction.
            traction_values: Surface traction amplitudes with shape `(2,)` or
                `(n_traction_faces, 2)`. Components are `[t_r, t_z]` with units
                `[force / area]`.
            nodal_temperature: Input-node temperatures with shape `(n_input_nodes,)` and units
                `[temperature]`. Required only when the model includes thermal materials.

        Returns:
            NDArray: Reduced right-hand side with shape `(ndof_reduced,)` and units
            `[generalized force] = [energy / distance]`.

        Raises:
            ValueError: If thermal materials are present but `nodal_temperature` is omitted.
        """

        body_force_arr = None if body_force is None else _normalize_body_force(body_force, self.nelem)
        npressure = int(self.pressure_faces.shape[0])
        pressure_arr = (
            None if pressure_values is None else _normalize_pressure_values(pressure_values, npressure)
        )
        traction_arr = (
            None
            if traction_values is None
            else _normalize_traction_values(
                traction_values,
                int(self.traction_faces.shape[0]),
            )
        )
        temperature_arr = self._normalize_temperature_for_backend(nodal_temperature)
        rhs = self._backend.build_rhs(
            None if body_force_arr is None else body_force_arr.reshape(-1),
            pressure_arr,
            None if traction_arr is None else traction_arr.reshape(-1),
            temperature_arr,
        )
        return np.asarray(rhs, dtype=np.float64)

    def solve(self, rhs: ArrayLike) -> npt.NDArray[np.floating[Any]]:
        """Solve the reduced system and recover the full displacement field.

        Args:
            rhs: Reduced right-hand side with shape `(ndof_reduced,)` and units
                `[generalized force] = [energy / distance]`.

        Returns:
            NDArray: Full displacement vector with shape `(ndof_full,)` and component ordering
            `[u_r0, u_z0, u_r1, u_z1, ...]`. Units are `[length]`.
        """

        rhs_arr = np.asarray(rhs).reshape(-1)
        assert (
            rhs_arr.shape[0] == self.ndof_reduced
        ), f"rhs must have length {self.ndof_reduced}; got {rhs_arr.shape}"
        return np.asarray(self._backend.solve(rhs_arr), dtype=np.float64)

    def recover_full(self, reduced_solution: ArrayLike) -> npt.NDArray[np.floating[Any]]:
        """Reinsert prescribed Dirichlet values into a reduced displacement vector.

        Args:
            reduced_solution: Reduced displacement vector with shape `(ndof_reduced,)` and units
                `[length]`.

        Returns:
            NDArray: Full displacement vector with shape `(ndof_full,)` and component ordering
            `[u_r0, u_z0, u_r1, u_z1, ...]`. Units are `[length]`.
        """

        reduced_arr = np.asarray(reduced_solution).reshape(-1)
        assert (
            reduced_arr.shape[0] == self.ndof_reduced
        ), f"reduced_solution must have length {self.ndof_reduced}; got {reduced_arr.shape}"
        full = np.zeros((self.ndof_full,), dtype=np.float64)
        full[self.fixed_dofs] = self.fixed_values
        full[self.free_dofs] = reduced_arr
        return full

    def _full_displacement_for_backend(self, displacements: ArrayLike) -> npt.NDArray[np.floating[Any]]:
        """Normalize reduced or full displacements to the backend's full flat vector."""

        arr = np.asarray(displacements)
        if arr.ndim == 1 and arr.shape == (self.ndof_reduced,):
            return self.recover_full(arr)
        return _normalize_displacements(displacements, self.analysis_nodes.shape[0]).reshape(-1)

    def strain(
        self,
        locations: PointLocations,
        displacements: ArrayLike,
    ) -> npt.NDArray[np.floating[Any]]:
        """Evaluate total strain at located points without materializing recovery matrices.

        Args:
            locations: Element-owned point locations from `quadrature().locations`,
                `locate_points(...)`, `locate_points_in_elements(...)`, or
                `QuadMeshQuery.point_locations()`.
            displacements: Either the reduced displacement solution with shape
                `(ndof_reduced,)`, or the full analysis displacement field with shape
                `(2 * n_analysis_nodes,)` or `(n_analysis_nodes, 2)`. Displacement units are
                `[length]`.

        Returns:
            NDArray: Total strain with shape `(npoint, 4)` and component ordering `[rr, zz, tt,
            rz]` for axisymmetric models or `[xx, yy, zz, xy]` for plane strain. Strain is
            unitless. For quadrature locations, reshape as `(nelem, nq_per_element, 4)` when an
            element-major view is needed.
        """

        locations = self._validate_locations(locations)
        displacements_full = self._full_displacement_for_backend(displacements)
        strain_flat = self._backend.strain(
            locations.element_indices.astype(np.uint64, copy=False),
            locations.reference_points,
            displacements_full,
        )
        return np.asarray(strain_flat, dtype=np.float64).reshape(-1, 4)

    def stress(
        self,
        locations: PointLocations,
        displacements: ArrayLike,
    ) -> npt.NDArray[np.floating[Any]]:
        """Evaluate stress at located points without materializing recovery matrices.

        Args:
            locations: Element-owned point locations from `quadrature().locations`,
                `locate_points(...)`, `locate_points_in_elements(...)`, or
                `QuadMeshQuery.point_locations()`.
            displacements: Either the reduced displacement solution with shape
                `(ndof_reduced,)`, or the full analysis displacement field with shape
                `(2 * n_analysis_nodes,)` or `(n_analysis_nodes, 2)`. Displacement units are
                `[length]`.

        Returns:
            NDArray: Stress with shape `(npoint, 4)` and component ordering `[rr, zz, tt, rz]` for
            axisymmetric models or `[xx, yy, zz, xy]` for plane strain. Units are `[stress]`. For
            quadrature locations, reshape as `(nelem, nq_per_element, 4)` when an element-major
            view is needed.
        """

        locations = self._validate_locations(locations)
        displacements_full = self._full_displacement_for_backend(displacements)
        stress_flat = self._backend.stress(
            locations.element_indices.astype(np.uint64, copy=False),
            locations.reference_points,
            displacements_full,
        )
        return np.asarray(stress_flat, dtype=np.float64).reshape(-1, 4)

    def thermal_strain(
        self,
        locations: PointLocations,
        nodal_temperature: ArrayLike | None = None,
    ) -> npt.NDArray[np.floating[Any]]:
        """Evaluate thermal strain at located points without materializing recovery matrices.

        Args:
            locations: Element-owned point locations from `quadrature().locations`,
                `locate_points(...)`, `locate_points_in_elements(...)`, or
                `QuadMeshQuery.point_locations()`.
            nodal_temperature: Input-node temperatures with shape `(n_input_nodes,)` and units
                `[temperature]`. Required only when the model includes thermal materials.

        Returns:
            NDArray: Thermal strain with shape `(npoint, 4)` and component ordering `[rr, zz, tt,
            rz]` for axisymmetric models or `[xx, yy, zz, xy]` for plane strain. Strain is
            unitless. Models without thermal materials return zeros and do not require
            `nodal_temperature`.
        """

        locations = self._validate_locations(locations)
        temperature_arr = self._normalize_temperature_for_backend(nodal_temperature)
        thermal_strain_flat = self._backend.thermal_strain(
            locations.element_indices.astype(np.uint64, copy=False),
            locations.reference_points,
            temperature_arr,
        )
        return np.asarray(thermal_strain_flat, dtype=np.float64).reshape(-1, 4)

    def thermal_stress(
        self,
        locations: PointLocations,
        nodal_temperature: ArrayLike | None = None,
    ) -> npt.NDArray[np.floating[Any]]:
        """Evaluate thermal stress at located points without materializing recovery matrices.

        Args:
            locations: Element-owned point locations from `quadrature().locations`,
                `locate_points(...)`, `locate_points_in_elements(...)`, or
                `QuadMeshQuery.point_locations()`.
            nodal_temperature: Input-node temperatures with shape `(n_input_nodes,)` and units
                `[temperature]`. Required only when the model includes thermal materials.

        Returns:
            NDArray: Thermal stress with shape `(npoint, 4)` and component ordering `[rr, zz, tt,
            rz]` for axisymmetric models or `[xx, yy, zz, xy]` for plane strain. Units are
            `[stress]`. Models without thermal materials return zeros and do not require
            `nodal_temperature`.
        """

        locations = self._validate_locations(locations)
        temperature_arr = self._normalize_temperature_for_backend(nodal_temperature)
        thermal_stress_flat = self._backend.thermal_stress(
            locations.element_indices.astype(np.uint64, copy=False),
            locations.reference_points,
            temperature_arr,
        )
        return np.asarray(thermal_stress_flat, dtype=np.float64).reshape(-1, 4)
body_force_to_rhs property 
body_force_to_rhs: csr_matrix

Operator mapping per-element body-force density to the reduced RHS.

Shape is (ndof_reduced, 2 * nelem). Entries have units [volume].

constant_rhs cached property 
constant_rhs: Float64Array

Load-independent reduced RHS contribution, exported from Rust on first access.

input_elements property 
input_elements: NDArray[uint64]

Input mesh connectivity with shape (nelem, 4).

input_nodes property 
input_nodes: Float64Array

Corner-node input mesh coordinates with shape (nnode, 2) and units [length].

pressure_to_rhs property 
pressure_to_rhs: csr_matrix

Operator mapping scalar pressure amplitudes to the reduced RHS.

Shape is (ndof_reduced, n_pressure_faces). Entries have units [area].

stiffness cached property 
stiffness: csc_matrix

Reduced stiffness matrix with shape (ndof_reduced, ndof_reduced).

Entries have units [generalized force / displacement] = [energy / distance^2]. The SciPy matrix is exported from the Rust backend on first access and then cached.

temperature_to_rhs property 
temperature_to_rhs: csr_matrix

Operator mapping input-node temperatures to the reduced RHS.

Shape is (ndof_reduced, n_temperature_nodes). Entries have units [generalized force / temperature] = [energy / (distance * temperature)].

traction_to_rhs property 
traction_to_rhs: csr_matrix

Operator mapping vector traction amplitudes to the reduced RHS.

Shape is (ndof_reduced, 2 * n_traction_faces). Entries have units [area].

build_rhs 
build_rhs(
    body_force: ArrayLike | None = None,
    pressure_values: ArrayLike | None = None,
    traction_values: ArrayLike | None = None,
    nodal_temperature: ArrayLike | None = None,
) -> npt.NDArray[np.floating[Any]]

Build one reduced structural right-hand side.

Parameters:

Name Type Description Default
body_force ArrayLike | None

Elementwise body-force amplitudes with shape (2,) or (nelem, 2). Components are [b_r, b_z] with units [force / volume].

 None
pressure_values ArrayLike | None

Pressure amplitudes with shape (n_pressure_faces,) and units [force / area]. Positive values act in the inward normal direction.

 None
traction_values ArrayLike | None

Surface traction amplitudes with shape (2,) or (n_traction_faces, 2). Components are [t_r, t_z] with units [force / area].

 None
nodal_temperature ArrayLike | None

Input-node temperatures with shape (n_input_nodes,) and units [temperature]. Required only when the model includes thermal materials.

 None

Returns:

Name Type Description
NDArray NDArray[floating[Any]]

Reduced right-hand side with shape (ndof_reduced,) and units
NDArray[floating[Any]]

[generalized force] = [energy / distance].

Raises:

Type Description
ValueError

If thermal materials are present but nodal_temperature is omitted.
Source code in cfsem/solenoid_stress/fem2d.py 
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
def build_rhs(
    self,
    body_force: ArrayLike | None = None,
    pressure_values: ArrayLike | None = None,
    traction_values: ArrayLike | None = None,
    nodal_temperature: ArrayLike | None = None,
) -> npt.NDArray[np.floating[Any]]:
    """Build one reduced structural right-hand side.

    Args:
        body_force: Elementwise body-force amplitudes with shape `(2,)` or `(nelem, 2)`.
            Components are `[b_r, b_z]` with units `[force / volume]`.
        pressure_values: Pressure amplitudes with shape `(n_pressure_faces,)` and units
            `[force / area]`. Positive values act in the inward normal direction.
        traction_values: Surface traction amplitudes with shape `(2,)` or
            `(n_traction_faces, 2)`. Components are `[t_r, t_z]` with units
            `[force / area]`.
        nodal_temperature: Input-node temperatures with shape `(n_input_nodes,)` and units
            `[temperature]`. Required only when the model includes thermal materials.

    Returns:
        NDArray: Reduced right-hand side with shape `(ndof_reduced,)` and units
        `[generalized force] = [energy / distance]`.

    Raises:
        ValueError: If thermal materials are present but `nodal_temperature` is omitted.
    """

    body_force_arr = None if body_force is None else _normalize_body_force(body_force, self.nelem)
    npressure = int(self.pressure_faces.shape[0])
    pressure_arr = (
        None if pressure_values is None else _normalize_pressure_values(pressure_values, npressure)
    )
    traction_arr = (
        None
        if traction_values is None
        else _normalize_traction_values(
            traction_values,
            int(self.traction_faces.shape[0]),
        )
    )
    temperature_arr = self._normalize_temperature_for_backend(nodal_temperature)
    rhs = self._backend.build_rhs(
        None if body_force_arr is None else body_force_arr.reshape(-1),
        pressure_arr,
        None if traction_arr is None else traction_arr.reshape(-1),
        temperature_arr,
    )
    return np.asarray(rhs, dtype=np.float64)
element_measures 
element_measures() -> ElementMeasures

Return cross-section area and represented volume for each element.

Returns:

Name Type Description
ElementMeasures ElementMeasures

Per-element measures with: areas of shape (nelem,) and units [area], volumes of shape (nelem,) and units [volume].
Source code in cfsem/solenoid_stress/fem2d.py 
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
def element_measures(self) -> ElementMeasures:
    """Return cross-section area and represented volume for each element.

    Returns:
        ElementMeasures: Per-element measures with:
            `areas` of shape `(nelem,)` and units `[area]`,
            `volumes` of shape `(nelem,)` and units `[volume]`.
    """

    cache = self._element_measures_cache
    if cache is not None:
        return cache
    quadrature = self.quadrature()
    weights_area = quadrature.weights_area.reshape(self.nelem, quadrature.points_per_element)
    weights_volume = quadrature.weights_volume.reshape(self.nelem, quadrature.points_per_element)
    cache = ElementMeasures(
        areas=np.asarray(weights_area.sum(axis=1), dtype=np.float64),
        volumes=np.asarray(weights_volume.sum(axis=1), dtype=np.float64),
    )
    self._element_measures_cache = cache
    return cache
interpolation_operator 
interpolation_operator(
    locations: PointLocations,
) -> sp.csr_matrix

Build a sparse interpolation operator for located points.

Parameters:

Name Type Description Default
locations PointLocations

Element-owned point locations from quadrature().locations, locate_points(...), locate_points_in_elements(...), or QuadMeshQuery.point_locations().

 required

Returns:

Name Type Description
csr_matrix csr_matrix

Sparse operator with shape (npoint, n_analysis_nodes). Multiplying by a
csr_matrix

scalar nodal field with shape (n_analysis_nodes,) returns interpolated values with
csr_matrix

shape (npoint,); multiplying by (n_analysis_nodes, ncomponent) interpolates each
csr_matrix

component independently. Entries are unitless shape-function values.
Source code in cfsem/solenoid_stress/fem2d.py 
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
def interpolation_operator(self, locations: PointLocations) -> sp.csr_matrix:
    """Build a sparse interpolation operator for located points.

    Args:
        locations: Element-owned point locations from `quadrature().locations`,
            `locate_points(...)`, `locate_points_in_elements(...)`, or
            `QuadMeshQuery.point_locations()`.

    Returns:
        csr_matrix: Sparse operator with shape `(npoint, n_analysis_nodes)`. Multiplying by a
        scalar nodal field with shape `(n_analysis_nodes,)` returns interpolated values with
        shape `(npoint,)`; multiplying by `(n_analysis_nodes, ncomponent)` interpolates each
        component independently. Entries are unitless shape-function values.
    """

    locations = self._validate_locations(locations)
    return _coo_operator_from_binding(
        _quad_mesh_interpolation_operator_f64(
            self.analysis_nodes,
            self.analysis_elements,
            locations.element_indices.astype(np.uint64, copy=False),
            locations.reference_points,
            self.element_type,
        ),
    )
locate_points 
locate_points(
    points: ArrayLike,
    *,
    outside: str = "nearest",
    tolerance: float | None = None,
    max_iterations: int = 20,
) -> PointLocations

Locate arbitrary physical points in this model's analysis mesh.

This uses the current brute-force quadrilateral mesh query and returns the nearest element plus reference coordinates for each query point. Points outside the mesh are projected to the nearest element unless outside="raise" or outside="error" is supplied.

Parameters:

Name Type Description Default
points ArrayLike

Physical coordinates with shape (npoint, 2) and units [length].

 required
outside str

Outside-mesh policy. "nearest" returns nearest-element projections; "raise" and "error" raise if any point is outside tolerance.

 'nearest'
tolerance float | None

Nonnegative physical distance used to classify contained points. Defaults to 1e-10.

 None
max_iterations int

Maximum local inverse-map iterations per element during the query.

 20

Returns:

Name Type Description
PointLocations PointLocations

Located points for recovery and sparse operator construction.

Raises:

Type Description
ValueError

If outside requests an error and at least one point is outside the mesh.
Source code in cfsem/solenoid_stress/fem2d.py 
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
def locate_points(
    self,
    points: ArrayLike,
    *,
    outside: str = "nearest",
    tolerance: float | None = None,
    max_iterations: int = 20,
) -> PointLocations:
    """Locate arbitrary physical points in this model's analysis mesh.

    This uses the current brute-force quadrilateral mesh query and returns the nearest element
    plus reference coordinates for each query point. Points outside the mesh are projected to
    the nearest element unless `outside="raise"` or `outside="error"` is supplied.

    Args:
        points: Physical coordinates with shape `(npoint, 2)` and units `[length]`.
        outside: Outside-mesh policy. `"nearest"` returns nearest-element projections;
            `"raise"` and `"error"` raise if any point is outside `tolerance`.
        tolerance: Nonnegative physical distance used to classify contained points. Defaults
            to `1e-10`.
        max_iterations: Maximum local inverse-map iterations per element during the query.

    Returns:
        PointLocations: Located points for recovery and sparse operator construction.

    Raises:
        ValueError: If `outside` requests an error and at least one point is outside the mesh.
    """

    query = query_quad_mesh(
        self.analysis_nodes,
        self.analysis_elements,
        points,
        element_type=self.element_type,
        max_iterations=max_iterations,
    )
    outside_policy = str(outside).strip().lower()
    assert outside_policy in {
        "nearest",
        "raise",
        "error",
    }, f"unsupported outside policy {outside!r}; use 'nearest' or 'raise'"
    tol = _normalize_query_tolerance(tolerance)
    inside = query.nearest_element_distances <= tol
    if outside_policy in {"raise", "error"} and not np.all(inside):
        first = int(np.flatnonzero(~inside)[0])
        raise ValueError(f"query point {first} is outside the quad mesh")
    return PointLocations(
        points=query.nearest_element_points,
        element_indices=query.nearest_element_indices,
        reference_points=query.nearest_element_reference_points,
        element_type=self.element_type,
    )
locate_points_in_elements 
locate_points_in_elements(
    points: ArrayLike,
    element_indices: ArrayLike,
    *,
    max_iterations: int = 20,
) -> PointLocations

Project physical points into caller-supplied owning elements.

This is the fast path when the caller already knows element ownership, such as when reusing element indices returned from quadrature() or a previous mesh query. It performs one local element projection per point and does not scan the global mesh.

Parameters:

Name Type Description Default
points ArrayLike

Physical coordinates with shape (npoint, 2) and units [length].

 required
element_indices ArrayLike

Owning analysis-element indices with shape (npoint,).

 required
max_iterations int

Maximum inverse-map iterations for each local element projection.

 20

Returns:

Name Type Description
PointLocations PointLocations

Projected physical points, caller-supplied element indices, and
PointLocations

reference coordinates.
Source code in cfsem/solenoid_stress/fem2d.py 
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
def locate_points_in_elements(
    self,
    points: ArrayLike,
    element_indices: ArrayLike,
    *,
    max_iterations: int = 20,
) -> PointLocations:
    """Project physical points into caller-supplied owning elements.

    This is the fast path when the caller already knows element ownership, such as when
    reusing element indices returned from `quadrature()` or a previous mesh query. It performs
    one local element projection per point and does not scan the global mesh.

    Args:
        points: Physical coordinates with shape `(npoint, 2)` and units `[length]`.
        element_indices: Owning analysis-element indices with shape `(npoint,)`.
        max_iterations: Maximum inverse-map iterations for each local element projection.

    Returns:
        PointLocations: Projected physical points, caller-supplied element indices, and
        reference coordinates.
    """

    points_arr = _normalize_query_points(points)
    element_indices_arr = np.asarray(element_indices, dtype=np.uint64).reshape(-1)
    projected_points, projected_elements, reference_points = self._backend.locate_points_in_elements(
        points_arr,
        element_indices_arr,
        int(max_iterations),
    )
    return PointLocations(
        points=np.asarray(projected_points, dtype=np.float64).reshape(-1, 2),
        element_indices=np.asarray(projected_elements, dtype=np.int64),
        reference_points=np.asarray(reference_points, dtype=np.float64).reshape(-1, 2),
        element_type=self.element_type,
    )
quadrature 
quadrature() -> Quadrature

Return element-major quadrature locations and mapped integration weights.

The returned locations are built directly from the model's quadrature rule and element geometry, so no global mesh query or point inversion is performed. Results are cached because the model mesh, quadrature rule, and geometry are immutable after assembly.

Returns:

Name Type Description
Quadrature Quadrature

Flat element-major locations with nelem * nq_per_element rows plus
Quadrature

mapped area and volume integration weights. Pass quadrature.locations to recovery
Quadrature

and sparse operator methods.
Source code in cfsem/solenoid_stress/fem2d.py 
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
def quadrature(self) -> Quadrature:
    """Return element-major quadrature locations and mapped integration weights.

    The returned locations are built directly from the model's quadrature rule and element
    geometry, so no global mesh query or point inversion is performed. Results are cached
    because the model mesh, quadrature rule, and geometry are immutable after assembly.

    Returns:
        Quadrature: Flat element-major locations with `nelem * nq_per_element` rows plus
        mapped area and volume integration weights. Pass `quadrature.locations` to recovery
        and sparse operator methods.
    """

    cache = self._quadrature_cache
    if cache is None:
        points, element_indices, reference_points, weights_area, weights_volume, points_per_element = (
            self._backend.quadrature()
        )
        locations = PointLocations(
            points=np.asarray(points, dtype=np.float64).reshape(-1, 2),
            element_indices=np.asarray(element_indices, dtype=np.int64),
            reference_points=np.asarray(reference_points, dtype=np.float64).reshape(-1, 2),
            element_type=self.element_type,
        )
        cache = Quadrature(
            locations=locations,
            weights_area=np.asarray(weights_area, dtype=np.float64),
            weights_volume=np.asarray(weights_volume, dtype=np.float64),
            points_per_element=int(points_per_element),
        )
        self._quadrature_cache = cache
    return cache
recover_full 
recover_full(
    reduced_solution: ArrayLike,
) -> npt.NDArray[np.floating[Any]]

Reinsert prescribed Dirichlet values into a reduced displacement vector.

Parameters:

Name Type Description Default
reduced_solution ArrayLike

Reduced displacement vector with shape (ndof_reduced,) and units [length].

 required

Returns:

Name Type Description
NDArray NDArray[floating[Any]]

Full displacement vector with shape (ndof_full,) and component ordering
NDArray[floating[Any]]

[u_r0, u_z0, u_r1, u_z1, ...]. Units are [length].
Source code in cfsem/solenoid_stress/fem2d.py 
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
def recover_full(self, reduced_solution: ArrayLike) -> npt.NDArray[np.floating[Any]]:
    """Reinsert prescribed Dirichlet values into a reduced displacement vector.

    Args:
        reduced_solution: Reduced displacement vector with shape `(ndof_reduced,)` and units
            `[length]`.

    Returns:
        NDArray: Full displacement vector with shape `(ndof_full,)` and component ordering
        `[u_r0, u_z0, u_r1, u_z1, ...]`. Units are `[length]`.
    """

    reduced_arr = np.asarray(reduced_solution).reshape(-1)
    assert (
        reduced_arr.shape[0] == self.ndof_reduced
    ), f"reduced_solution must have length {self.ndof_reduced}; got {reduced_arr.shape}"
    full = np.zeros((self.ndof_full,), dtype=np.float64)
    full[self.fixed_dofs] = self.fixed_values
    full[self.free_dofs] = reduced_arr
    return full
solve 
solve(rhs: ArrayLike) -> npt.NDArray[np.floating[Any]]

Solve the reduced system and recover the full displacement field.

Parameters:

Name Type Description Default
rhs ArrayLike

Reduced right-hand side with shape (ndof_reduced,) and units [generalized force] = [energy / distance].

 required

Returns:

Name Type Description
NDArray NDArray[floating[Any]]

Full displacement vector with shape (ndof_full,) and component ordering
NDArray[floating[Any]]

[u_r0, u_z0, u_r1, u_z1, ...]. Units are [length].
Source code in cfsem/solenoid_stress/fem2d.py 
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
def solve(self, rhs: ArrayLike) -> npt.NDArray[np.floating[Any]]:
    """Solve the reduced system and recover the full displacement field.

    Args:
        rhs: Reduced right-hand side with shape `(ndof_reduced,)` and units
            `[generalized force] = [energy / distance]`.

    Returns:
        NDArray: Full displacement vector with shape `(ndof_full,)` and component ordering
        `[u_r0, u_z0, u_r1, u_z1, ...]`. Units are `[length]`.
    """

    rhs_arr = np.asarray(rhs).reshape(-1)
    assert (
        rhs_arr.shape[0] == self.ndof_reduced
    ), f"rhs must have length {self.ndof_reduced}; got {rhs_arr.shape}"
    return np.asarray(self._backend.solve(rhs_arr), dtype=np.float64)
strain 
strain(
    locations: PointLocations, displacements: ArrayLike
) -> npt.NDArray[np.floating[Any]]

Evaluate total strain at located points without materializing recovery matrices.

Parameters:

Name Type Description Default
locations PointLocations

Element-owned point locations from quadrature().locations, locate_points(...), locate_points_in_elements(...), or QuadMeshQuery.point_locations().

 required
displacements ArrayLike

Either the reduced displacement solution with shape (ndof_reduced,), or the full analysis displacement field with shape (2 * n_analysis_nodes,) or (n_analysis_nodes, 2). Displacement units are [length].

 required

Returns:

Name Type Description
NDArray NDArray[floating[Any]]

Total strain with shape (npoint, 4) and component ordering `[rr, zz, tt,
NDArray[floating[Any]]

rz]for axisymmetric models or[xx, yy, zz, xy]` for plane strain. Strain is
NDArray[floating[Any]]

unitless. For quadrature locations, reshape as (nelem, nq_per_element, 4) when an
NDArray[floating[Any]]

element-major view is needed.
Source code in cfsem/solenoid_stress/fem2d.py 
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
def strain(
    self,
    locations: PointLocations,
    displacements: ArrayLike,
) -> npt.NDArray[np.floating[Any]]:
    """Evaluate total strain at located points without materializing recovery matrices.

    Args:
        locations: Element-owned point locations from `quadrature().locations`,
            `locate_points(...)`, `locate_points_in_elements(...)`, or
            `QuadMeshQuery.point_locations()`.
        displacements: Either the reduced displacement solution with shape
            `(ndof_reduced,)`, or the full analysis displacement field with shape
            `(2 * n_analysis_nodes,)` or `(n_analysis_nodes, 2)`. Displacement units are
            `[length]`.

    Returns:
        NDArray: Total strain with shape `(npoint, 4)` and component ordering `[rr, zz, tt,
        rz]` for axisymmetric models or `[xx, yy, zz, xy]` for plane strain. Strain is
        unitless. For quadrature locations, reshape as `(nelem, nq_per_element, 4)` when an
        element-major view is needed.
    """

    locations = self._validate_locations(locations)
    displacements_full = self._full_displacement_for_backend(displacements)
    strain_flat = self._backend.strain(
        locations.element_indices.astype(np.uint64, copy=False),
        locations.reference_points,
        displacements_full,
    )
    return np.asarray(strain_flat, dtype=np.float64).reshape(-1, 4)
strain_operator 
strain_operator(locations: PointLocations) -> sp.csr_matrix

Build a sparse total-strain recovery operator for located points.

Parameters:

Name Type Description Default
locations PointLocations

Element-owned point locations from quadrature().locations, locate_points(...), locate_points_in_elements(...), or QuadMeshQuery.point_locations().

 required

Returns:

Name Type Description
csr_matrix csr_matrix

Sparse operator with shape (4 * npoint, 2 * n_analysis_nodes). Rows are
csr_matrix

grouped by point and tensor component. Multiplying by full analysis displacements with
csr_matrix

shape (2 * n_analysis_nodes,) returns flat strain samples with shape (4 * npoint,).
csr_matrix

Entries have units [1 / length].
Source code in cfsem/solenoid_stress/fem2d.py 
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
def strain_operator(self, locations: PointLocations) -> sp.csr_matrix:
    """Build a sparse total-strain recovery operator for located points.

    Args:
        locations: Element-owned point locations from `quadrature().locations`,
            `locate_points(...)`, `locate_points_in_elements(...)`, or
            `QuadMeshQuery.point_locations()`.

    Returns:
        csr_matrix: Sparse operator with shape `(4 * npoint, 2 * n_analysis_nodes)`. Rows are
        grouped by point and tensor component. Multiplying by full analysis displacements with
        shape `(2 * n_analysis_nodes,)` returns flat strain samples with shape `(4 * npoint,)`.
        Entries have units `[1 / length]`.
    """

    locations = self._validate_locations(locations)
    return _coo_operator_from_binding(
        _quad_mesh_strain_operator_f64(
            self.analysis_nodes,
            self.analysis_elements,
            locations.element_indices.astype(np.uint64, copy=False),
            locations.reference_points,
            self.element_type,
            _formulation_code(self.formulation),
            self.thickness,
        ),
    )
stress 
stress(
    locations: PointLocations, displacements: ArrayLike
) -> npt.NDArray[np.floating[Any]]

Evaluate stress at located points without materializing recovery matrices.

Parameters:

Name Type Description Default
locations PointLocations

Element-owned point locations from quadrature().locations, locate_points(...), locate_points_in_elements(...), or QuadMeshQuery.point_locations().

 required
displacements ArrayLike

Either the reduced displacement solution with shape (ndof_reduced,), or the full analysis displacement field with shape (2 * n_analysis_nodes,) or (n_analysis_nodes, 2). Displacement units are [length].

 required

Returns:

Name Type Description
NDArray NDArray[floating[Any]]

Stress with shape (npoint, 4) and component ordering [rr, zz, tt, rz] for
NDArray[floating[Any]]

axisymmetric models or [xx, yy, zz, xy] for plane strain. Units are [stress]. For
NDArray[floating[Any]]

quadrature locations, reshape as (nelem, nq_per_element, 4) when an element-major
NDArray[floating[Any]]

view is needed.
Source code in cfsem/solenoid_stress/fem2d.py 
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
def stress(
    self,
    locations: PointLocations,
    displacements: ArrayLike,
) -> npt.NDArray[np.floating[Any]]:
    """Evaluate stress at located points without materializing recovery matrices.

    Args:
        locations: Element-owned point locations from `quadrature().locations`,
            `locate_points(...)`, `locate_points_in_elements(...)`, or
            `QuadMeshQuery.point_locations()`.
        displacements: Either the reduced displacement solution with shape
            `(ndof_reduced,)`, or the full analysis displacement field with shape
            `(2 * n_analysis_nodes,)` or `(n_analysis_nodes, 2)`. Displacement units are
            `[length]`.

    Returns:
        NDArray: Stress with shape `(npoint, 4)` and component ordering `[rr, zz, tt, rz]` for
        axisymmetric models or `[xx, yy, zz, xy]` for plane strain. Units are `[stress]`. For
        quadrature locations, reshape as `(nelem, nq_per_element, 4)` when an element-major
        view is needed.
    """

    locations = self._validate_locations(locations)
    displacements_full = self._full_displacement_for_backend(displacements)
    stress_flat = self._backend.stress(
        locations.element_indices.astype(np.uint64, copy=False),
        locations.reference_points,
        displacements_full,
    )
    return np.asarray(stress_flat, dtype=np.float64).reshape(-1, 4)
stress_operator 
stress_operator(locations: PointLocations) -> sp.csr_matrix

Build a sparse elastic-stress recovery operator for located points.

Parameters:

Name Type Description Default
locations PointLocations

Element-owned point locations from quadrature().locations, locate_points(...), locate_points_in_elements(...), or QuadMeshQuery.point_locations().

 required

Returns:

Name Type Description
csr_matrix csr_matrix

Sparse operator with shape (4 * npoint, 2 * n_analysis_nodes). Rows are
csr_matrix

grouped by point and tensor component. Multiplying by full analysis displacements with
csr_matrix

shape (2 * n_analysis_nodes,) returns flat stress samples with shape (4 * npoint,).
csr_matrix

Entries have units [stress / length].
Source code in cfsem/solenoid_stress/fem2d.py 
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
def stress_operator(self, locations: PointLocations) -> sp.csr_matrix:
    """Build a sparse elastic-stress recovery operator for located points.

    Args:
        locations: Element-owned point locations from `quadrature().locations`,
            `locate_points(...)`, `locate_points_in_elements(...)`, or
            `QuadMeshQuery.point_locations()`.

    Returns:
        csr_matrix: Sparse operator with shape `(4 * npoint, 2 * n_analysis_nodes)`. Rows are
        grouped by point and tensor component. Multiplying by full analysis displacements with
        shape `(2 * n_analysis_nodes,)` returns flat stress samples with shape `(4 * npoint,)`.
        Entries have units `[stress / length]`.
    """

    locations = self._validate_locations(locations)
    return _coo_operator_from_binding(
        _quad_mesh_stress_operator_f64(
            self.analysis_nodes,
            self.analysis_elements,
            locations.element_indices.astype(np.uint64, copy=False),
            locations.reference_points,
            self._material_ids,
            self._material_table,
            self._material_orientation_angles,
            self.element_type,
            _formulation_code(self.formulation),
            self.thickness,
        ),
    )
thermal_strain 
thermal_strain(
    locations: PointLocations,
    nodal_temperature: ArrayLike | None = None,
) -> npt.NDArray[np.floating[Any]]

Evaluate thermal strain at located points without materializing recovery matrices.

Parameters:

Name Type Description Default
locations PointLocations

Element-owned point locations from quadrature().locations, locate_points(...), locate_points_in_elements(...), or QuadMeshQuery.point_locations().

 required
nodal_temperature ArrayLike | None

Input-node temperatures with shape (n_input_nodes,) and units [temperature]. Required only when the model includes thermal materials.

 None

Returns:

Name Type Description
NDArray NDArray[floating[Any]]

Thermal strain with shape (npoint, 4) and component ordering `[rr, zz, tt,
NDArray[floating[Any]]

rz]for axisymmetric models or[xx, yy, zz, xy]` for plane strain. Strain is
NDArray[floating[Any]]

unitless. Models without thermal materials return zeros and do not require
NDArray[floating[Any]]

nodal_temperature.
Source code in cfsem/solenoid_stress/fem2d.py 
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
def thermal_strain(
    self,
    locations: PointLocations,
    nodal_temperature: ArrayLike | None = None,
) -> npt.NDArray[np.floating[Any]]:
    """Evaluate thermal strain at located points without materializing recovery matrices.

    Args:
        locations: Element-owned point locations from `quadrature().locations`,
            `locate_points(...)`, `locate_points_in_elements(...)`, or
            `QuadMeshQuery.point_locations()`.
        nodal_temperature: Input-node temperatures with shape `(n_input_nodes,)` and units
            `[temperature]`. Required only when the model includes thermal materials.

    Returns:
        NDArray: Thermal strain with shape `(npoint, 4)` and component ordering `[rr, zz, tt,
        rz]` for axisymmetric models or `[xx, yy, zz, xy]` for plane strain. Strain is
        unitless. Models without thermal materials return zeros and do not require
        `nodal_temperature`.
    """

    locations = self._validate_locations(locations)
    temperature_arr = self._normalize_temperature_for_backend(nodal_temperature)
    thermal_strain_flat = self._backend.thermal_strain(
        locations.element_indices.astype(np.uint64, copy=False),
        locations.reference_points,
        temperature_arr,
    )
    return np.asarray(thermal_strain_flat, dtype=np.float64).reshape(-1, 4)
thermal_stress 
thermal_stress(
    locations: PointLocations,
    nodal_temperature: ArrayLike | None = None,
) -> npt.NDArray[np.floating[Any]]

Evaluate thermal stress at located points without materializing recovery matrices.

Parameters:

Name Type Description Default
locations PointLocations

Element-owned point locations from quadrature().locations, locate_points(...), locate_points_in_elements(...), or QuadMeshQuery.point_locations().

 required
nodal_temperature ArrayLike | None

Input-node temperatures with shape (n_input_nodes,) and units [temperature]. Required only when the model includes thermal materials.

 None

Returns:

Name Type Description
NDArray NDArray[floating[Any]]

Thermal stress with shape (npoint, 4) and component ordering `[rr, zz, tt,
NDArray[floating[Any]]

rz]for axisymmetric models or[xx, yy, zz, xy]` for plane strain. Units are
NDArray[floating[Any]]

[stress]. Models without thermal materials return zeros and do not require
NDArray[floating[Any]]

nodal_temperature.
Source code in cfsem/solenoid_stress/fem2d.py 
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
def thermal_stress(
    self,
    locations: PointLocations,
    nodal_temperature: ArrayLike | None = None,
) -> npt.NDArray[np.floating[Any]]:
    """Evaluate thermal stress at located points without materializing recovery matrices.

    Args:
        locations: Element-owned point locations from `quadrature().locations`,
            `locate_points(...)`, `locate_points_in_elements(...)`, or
            `QuadMeshQuery.point_locations()`.
        nodal_temperature: Input-node temperatures with shape `(n_input_nodes,)` and units
            `[temperature]`. Required only when the model includes thermal materials.

    Returns:
        NDArray: Thermal stress with shape `(npoint, 4)` and component ordering `[rr, zz, tt,
        rz]` for axisymmetric models or `[xx, yy, zz, xy]` for plane strain. Units are
        `[stress]`. Models without thermal materials return zeros and do not require
        `nodal_temperature`.
    """

    locations = self._validate_locations(locations)
    temperature_arr = self._normalize_temperature_for_backend(nodal_temperature)
    thermal_stress_flat = self._backend.thermal_stress(
        locations.element_indices.astype(np.uint64, copy=False),
        locations.reference_points,
        temperature_arr,
    )
    return np.asarray(thermal_stress_flat, dtype=np.float64).reshape(-1, 4)

assemble_structural_2d 

assemble_structural_2d(
    nodes: ArrayLike,
    elements: ArrayLike,
    material_ids: ArrayLike,
    material_table: ArrayLike,
    *,
    formulation: str = "axisymmetric",
    thickness: float | None = None,
    material_orientation_angles: ArrayLike | None = None,
    pressure_faces: ArrayLike | None = None,
    traction_faces: ArrayLike | None = None,
    thermal_material_table: ArrayLike | None = None,
    prescribed: Mapping[int, float] | None = None,
    quadrature: str | int = "gl3",
    element_type: str = "quad4",
    par: bool = True,
) -> Structural2DFEMModel

Assemble the reusable 2D structural FEM model.

Parameters:

Name Type Description Default
nodes ArrayLike

Corner-node coordinates with shape (nnode, 2). Coordinates are (r, z) for formulation="axisymmetric" and (x, y) for formulation="plane_strain". Units are [length]. Floating input arrays must have dtype float64.

 required
elements ArrayLike

Connectivity with shape (nelem, 4) for element_type="quad4". For element_type="quad9", pass either corner-only (nelem, 4) connectivity to infer a straight-sided quad9 mesh, or explicit (nelem, 9) connectivity in local order [corner0, corner1, corner2, corner3, face0_mid, face1_mid, face2_mid, face3_mid, center]. Corner nodes must be ordered counter-clockwise in the 2D analysis plane.

 required
material_ids ArrayLike

Dense material row indices with shape (nelem,).

 required
material_table ArrayLike

Elastic stress-strain matrices with shape (nmat, 4, 4). Matrix units are [stress / strain] = [pressure].

 required
formulation str

Symmetry reduction, either "axisymmetric" or "plane_strain".

 'axisymmetric'
thickness float | None

Plane-strain out-of-plane thickness. Required only for formulation="plane_strain".

 None
material_orientation_angles ArrayLike | None

Optional scalar or per-element angles, in radians, rotating local material axes into the global 2D frame before assembly.

 None
pressure_faces ArrayLike | None

Optional pressure-load topology with shape (n_pressure_faces, 2). Each row is [element_index, local_face].

 None
traction_faces ArrayLike | None

Optional traction-load topology with shape (n_traction_faces, 2). Each row is [element_index, local_face].

 None
thermal_material_table ArrayLike | None

Optional thermal material rows with shape (nmat, 5) storing [alpha_r, alpha_z, alpha_t, alpha_rz, T_ref]. Thermal expansion coefficients have units [strain / temperature] and T_ref has units [temperature].

 None
prescribed Mapping[int, float] | None

Optional mapping from full displacement DOF index to prescribed displacement value. Displacement units are [length].

 None
quadrature str | int

Quadrature rule selector, either gl3, gl4, 3, or 4.

 'gl3'
element_type str

Analysis element family, either quad4 or quad9.

 'quad4'
par bool

Whether to assemble stiffness and computed-on-call sparse exports using threaded element batches.

 True

Returns:

Name Type Description
Structural2DFEMModel Structural2DFEMModel

Reusable model with backend solve state and user-owned sparse
Structural2DFEMModel

operator exports.

Raises:

Type Description
ValueError

If quadrature is unsupported.
AssertionError

If array shapes are invalid or if mapping-style material inputs are passed instead of dense arrays.
Source code in cfsem/solenoid_stress/fem2d.py 
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
def assemble_structural_2d(
    nodes: ArrayLike,
    elements: ArrayLike,
    material_ids: ArrayLike,
    material_table: ArrayLike,
    *,
    formulation: str = "axisymmetric",
    thickness: float | None = None,
    material_orientation_angles: ArrayLike | None = None,
    pressure_faces: ArrayLike | None = None,
    traction_faces: ArrayLike | None = None,
    thermal_material_table: ArrayLike | None = None,
    prescribed: Mapping[int, float] | None = None,
    quadrature: str | int = "gl3",
    element_type: str = "quad4",
    par: bool = True,
) -> Structural2DFEMModel:
    """Assemble the reusable 2D structural FEM model.

    Args:
        nodes: Corner-node coordinates with shape `(nnode, 2)`. Coordinates are `(r, z)` for
            `formulation="axisymmetric"` and `(x, y)` for `formulation="plane_strain"`.
            Units are `[length]`. Floating input arrays must have dtype `float64`.
        elements: Connectivity with shape `(nelem, 4)` for `element_type="quad4"`. For
            `element_type="quad9"`, pass either corner-only `(nelem, 4)` connectivity to infer a
            straight-sided quad9 mesh, or explicit `(nelem, 9)` connectivity in local order
            `[corner0, corner1, corner2, corner3, face0_mid, face1_mid, face2_mid, face3_mid,
            center]`. Corner nodes must be ordered counter-clockwise in the 2D analysis plane.
        material_ids: Dense material row indices with shape `(nelem,)`.
        material_table: Elastic stress-strain matrices with shape `(nmat, 4, 4)`. Matrix units
            are `[stress / strain] = [pressure]`.
        formulation: Symmetry reduction, either `"axisymmetric"` or `"plane_strain"`.
        thickness: Plane-strain out-of-plane thickness. Required only for
            `formulation="plane_strain"`.
        material_orientation_angles: Optional scalar or per-element angles, in radians, rotating
            local material axes into the global 2D frame before assembly.
        pressure_faces: Optional pressure-load topology with shape `(n_pressure_faces, 2)`. Each
            row is `[element_index, local_face]`.
        traction_faces: Optional traction-load topology with shape `(n_traction_faces, 2)`. Each
            row is `[element_index, local_face]`.
        thermal_material_table: Optional thermal material rows with shape `(nmat, 5)` storing
            `[alpha_r, alpha_z, alpha_t, alpha_rz, T_ref]`. Thermal expansion coefficients have
            units `[strain / temperature]` and `T_ref` has units `[temperature]`.
        prescribed: Optional mapping from full displacement DOF index to prescribed displacement
            value. Displacement units are `[length]`.
        quadrature: Quadrature rule selector, either `gl3`, `gl4`, `3`, or `4`.
        element_type: Analysis element family, either `quad4` or `quad9`.
        par: Whether to assemble stiffness and computed-on-call sparse exports using threaded
            element batches.

    Returns:
        Structural2DFEMModel: Reusable model with backend solve state and user-owned sparse
        operator exports.

    Raises:
        ValueError: If `quadrature` is unsupported.
        AssertionError: If array shapes are invalid or if mapping-style material inputs are passed
            instead of dense arrays.
    """

    nodes_arr = _normalize_nodes(nodes)
    material_ids_arr, material_table_arr = _normalize_materials(material_ids, material_table)
    thermal_material_table_arr = _normalize_thermal_material_table(
        thermal_material_table,
    )
    pressure_faces_arr = _normalize_face_pairs("pressure_faces", pressure_faces)
    traction_faces_arr = _normalize_face_pairs("traction_faces", traction_faces)
    prescribed_dofs, prescribed_values = _normalize_prescribed_dirichlet(prescribed)
    quadrature_code = _quadrature_code(quadrature)
    normalized_formulation = _normalize_formulation(formulation)
    thickness_value = _normalize_thickness(normalized_formulation, thickness)
    normalized_element_type = _normalize_element_type(element_type)
    if normalized_element_type == "quad4":
        elements_arr = _normalize_elements(elements, 4)
        analysis_nodes, analysis_elements, elevated = nodes_arr, elements_arr, None
    else:
        elements_arr = _normalize_elements(elements, (4, 9))
        if elements_arr.shape[1] == 4:
            elevated = infer_quad9_mesh(nodes_arr, elements_arr)
            analysis_nodes, analysis_elements = elevated.analysis_nodes, elevated.analysis_elements
        else:
            analysis_nodes, analysis_elements, elevated = nodes_arr, elements_arr, None
    material_orientation_angles_arr = _normalize_material_orientation_angles(
        material_orientation_angles,
        int(analysis_elements.shape[0]),
    )
    backend = _assemble_model_2d_f64(
        analysis_nodes,
        analysis_elements,
        material_ids_arr,
        material_table_arr,
        pressure_faces_arr,
        traction_faces_arr,
        (
            thermal_material_table_arr
            if thermal_material_table_arr is not None
            else np.zeros((0, 5), dtype=np.float64)
        ),
        material_orientation_angles_arr,
        prescribed_dofs,
        prescribed_values,
        4 if normalized_element_type == "quad4" else 9,
        _formulation_code(normalized_formulation),
        thickness_value,
        quadrature_code,
        par,
    )
    n_temperature_nodes = 0 if thermal_material_table_arr is None else nodes_arr.shape[0]

    model = Structural2DFEMModel(
        backend=backend,
        input_nodes=nodes_arr,
        input_elements=elements_arr,
        analysis_nodes=analysis_nodes,
        analysis_elements=analysis_elements,
        elevated=elevated,
        material_ids=material_ids_arr,
        material_table=material_table_arr,
        thermal_material_table=thermal_material_table_arr,
        material_orientation_angles=material_orientation_angles_arr,
        pressure_faces=pressure_faces_arr,
        traction_faces=traction_faces_arr,
        formulation=normalized_formulation,
        thickness=thickness_value,
        element_type=normalized_element_type,
        free_dofs=np.asarray(backend.free_dofs(), dtype=np.int64),
        fixed_dofs=np.asarray(backend.fixed_dofs(), dtype=np.int64),
        fixed_values=np.asarray(backend.fixed_values(), dtype=np.float64),
        ndof_full=int(backend.ndof_full),
        ndof_reduced=int(backend.ndof_reduced),
        nelem=elements_arr.shape[0],
        nq_per_element=int(backend.nq_per_element),
        n_temperature_nodes=n_temperature_nodes,
    )
    return model

cfsem_radial_material 

cfsem_radial_material(
    youngs_modulus: float, poisson_ratio: float
) -> npt.NDArray[np.floating[Any]]

Construct the reduced elastic matrix used by the 1D radial solver.

Parameters:

Name Type Description Default
youngs_modulus float

Young's modulus with units [pressure].

 required
poisson_ratio float

Poisson ratio with units [dimensionless].

 required

Returns:

Name Type Description
NDArray NDArray[floating[Any]]

Elastic stress-strain matrix with shape (4, 4) in component order
NDArray[floating[Any]]

[rr, zz, tt, rz]. Units are [stress / strain] = [pressure].
Source code in cfsem/solenoid_stress/fem2d.py 
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
def cfsem_radial_material(
    youngs_modulus: float,
    poisson_ratio: float,
) -> npt.NDArray[np.floating[Any]]:
    """Construct the reduced elastic matrix used by the 1D radial solver.

    Args:
        youngs_modulus: Young's modulus with units `[pressure]`.
        poisson_ratio: Poisson ratio with units `[dimensionless]`.

    Returns:
        NDArray: Elastic stress-strain matrix with shape `(4, 4)` in component order
        `[rr, zz, tt, rz]`. Units are `[stress / strain] = [pressure]`.
    """

    return _as_float64_array(_cfsem_radial_material_f64(youngs_modulus, poisson_ratio)).reshape(4, 4)

infer_quad9_mesh 

infer_quad9_mesh(
    nodes: ArrayLike, elements: ArrayLike
) -> ElevatedQuad9Mesh

Elevate a corner-only quad mesh to an explicit 9-node Lagrange mesh.

Parameters:

Name Type Description Default
nodes ArrayLike

Corner-node coordinates with shape (nnode, 2) in (r, z) order. Units are [length].

 required
elements ArrayLike

Quad4 connectivity with shape (nelem, 4). Corner nodes must be ordered counter-clockwise in the (r, z) plane.

 required

Returns:

Name Type Description
ElevatedQuad9Mesh ElevatedQuad9Mesh

Elevated analysis mesh with: analysis_nodes of shape (n_analysis_nodes, 2) and units [length], analysis_elements of shape (nelem, 9), corner_node_indices, midside_node_indices, and center_node_indices as one-dimensional index arrays.
Source code in cfsem/solenoid_stress/fem2d.py 
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
def infer_quad9_mesh(nodes: ArrayLike, elements: ArrayLike) -> ElevatedQuad9Mesh:
    """Elevate a corner-only quad mesh to an explicit 9-node Lagrange mesh.

    Args:
        nodes: Corner-node coordinates with shape `(nnode, 2)` in `(r, z)` order. Units are
            `[length]`.
        elements: Quad4 connectivity with shape `(nelem, 4)`. Corner nodes must be ordered
            counter-clockwise in the `(r, z)` plane.

    Returns:
        ElevatedQuad9Mesh: Elevated analysis mesh with:
            `analysis_nodes` of shape `(n_analysis_nodes, 2)` and units `[length]`,
            `analysis_elements` of shape `(nelem, 9)`,
            `corner_node_indices`, `midside_node_indices`, and `center_node_indices` as
            one-dimensional index arrays.
    """

    nodes_arr = _normalize_nodes(nodes)
    elements_arr = _normalize_elements(elements)
    (
        analysis_nodes_flat,
        analysis_elements_flat,
        corner_node_indices,
        midside_node_indices,
        center_node_indices,
    ) = _infer_quad9_mesh_f64(nodes_arr, elements_arr)

    return ElevatedQuad9Mesh(
        input_nodes=nodes_arr,
        input_elements=elements_arr,
        analysis_nodes=np.asarray(analysis_nodes_flat, dtype=np.float64).reshape(-1, 2),
        analysis_elements=np.asarray(analysis_elements_flat, dtype=np.uint64).reshape(-1, 9),
        corner_node_indices=np.asarray(corner_node_indices, dtype=np.int64),
        midside_node_indices=np.asarray(midside_node_indices, dtype=np.int64),
        center_node_indices=np.asarray(center_node_indices, dtype=np.int64),
    )

interpolate_quad_mesh_values 

interpolate_quad_mesh_values(
    nodes: ArrayLike,
    elements: ArrayLike,
    nodal_values: ArrayLike,
    points: ArrayLike,
    *,
    element_type: str = "quad4",
    outside: str = "raise",
    tolerance: float | None = None,
    max_iterations: int = 20,
) -> QuadMeshInterpolation

Interpolate nodal values at arbitrary physical points in a quadrilateral mesh.

The interpolation uses the element's actual shape functions. nodal_values may have shape (nnode,) or (nnode, ...); the returned values have shape (npoint,) or (npoint, ...).

Point location is Rust-backed but brute-force and scans all elements once per query point. Outside policies are applied from the nearest-element distance: "raise" errors, "nan" masks outside values, and "nearest" returns the nearest-element interpolation. Complexity is O(npoint * nelem * max_iterations) for point location plus O(npoint * nodes_per_element * ncomponent) for interpolation.

Parameters:

Name Type Description Default
nodes ArrayLike

Mesh node coordinates with shape (nnode, 2).

 required
elements ArrayLike

Quad connectivity with shape (nelem, 4) or (nelem, 9).

 required
nodal_values ArrayLike

Values at mesh nodes with shape (nnode,) or (nnode, ...).

 required
points ArrayLike

Query point coordinates with shape (npoint, 2).

 required
element_type str

Element family, either "quad4" or "quad9".

 'quad4'
outside str

Outside-mesh policy: "raise"/"error", "nan", or "nearest".

 'raise'
tolerance float | None

Physical and reference-space tolerance for point containment.

 None
max_iterations int

Maximum Newton/projection iterations per element.

 20

Returns:

Type Description
QuadMeshInterpolation

Interpolated values plus element indices, reference coordinates, and inside flags for the
QuadMeshInterpolation

query points. values has shape (npoint,) or (npoint, ...) and the same units as
QuadMeshInterpolation

nodal_values; element_indices is unitless with shape (npoint,); reference_points
QuadMeshInterpolation

is unitless with shape (npoint, 2); inside has shape (npoint,).
Source code in cfsem/solenoid_stress/fem2d.py 
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
def interpolate_quad_mesh_values(
    nodes: ArrayLike,
    elements: ArrayLike,
    nodal_values: ArrayLike,
    points: ArrayLike,
    *,
    element_type: str = "quad4",
    outside: str = "raise",
    tolerance: float | None = None,
    max_iterations: int = 20,
) -> QuadMeshInterpolation:
    """Interpolate nodal values at arbitrary physical points in a quadrilateral mesh.

    The interpolation uses the element's actual shape functions. `nodal_values` may have shape
    `(nnode,)` or `(nnode, ...)`; the returned values have shape `(npoint,)` or `(npoint, ...)`.

    Point location is Rust-backed but brute-force and scans all elements once per query point.
    Outside policies are applied from the nearest-element distance: `"raise"` errors,
    `"nan"` masks outside values, and `"nearest"` returns the nearest-element interpolation.
    Complexity is `O(npoint * nelem * max_iterations)` for point location plus
    `O(npoint * nodes_per_element * ncomponent)` for interpolation.

    Args:
        nodes: Mesh node coordinates with shape `(nnode, 2)`.
        elements: Quad connectivity with shape `(nelem, 4)` or `(nelem, 9)`.
        nodal_values: Values at mesh nodes with shape `(nnode,)` or `(nnode, ...)`.
        points: Query point coordinates with shape `(npoint, 2)`.
        element_type: Element family, either `"quad4"` or `"quad9"`.
        outside: Outside-mesh policy: `"raise"`/`"error"`, `"nan"`, or `"nearest"`.
        tolerance: Physical and reference-space tolerance for point containment.
        max_iterations: Maximum Newton/projection iterations per element.

    Returns:
        Interpolated values plus element indices, reference coordinates, and inside flags for the
        query points. `values` has shape `(npoint,)` or `(npoint, ...)` and the same units as
        `nodal_values`; `element_indices` is unitless with shape `(npoint,)`; `reference_points`
        is unitless with shape `(npoint, 2)`; `inside` has shape `(npoint,)`.
    """

    query = query_quad_mesh(
        nodes,
        elements,
        points,
        element_type=element_type,
        max_iterations=max_iterations,
    )
    values_arr = np.asarray(nodal_values)
    assert (
        values_arr.ndim >= 1 and values_arr.shape[0] == query.nodes.shape[0]
    ), f"nodal_values must have shape (nnode,) or (nnode, ...); got {values_arr.shape}"
    values_shape = values_arr.shape[1:]
    values_2d = values_arr.reshape(query.nodes.shape[0], -1)
    outside_policy = str(outside).strip().lower()
    assert outside_policy in {
        "nearest",
        "nan",
        "raise",
        "error",
    }, f"unsupported outside policy {outside!r}; use 'raise', 'nan', or 'nearest'"
    tol = _normalize_query_tolerance(tolerance)
    inside = query.nearest_element_distances <= tol
    if outside_policy in {"raise", "error"} and not np.all(inside):
        first = int(np.flatnonzero(~inside)[0])
        raise ValueError(f"query point {first} is outside the quad mesh")
    operator = quad_mesh_interpolation_operator(query)
    values = np.asarray(operator @ values_2d).reshape((query.points.shape[0], *values_shape))
    if outside_policy == "nan" and np.any(~inside):
        values[~inside] = np.nan
    return QuadMeshInterpolation(
        values=values,
        element_indices=query.nearest_element_indices,
        reference_points=query.nearest_element_reference_points,
        inside=inside,
    )

isotropic_axisymmetric_material 

isotropic_axisymmetric_material(
    youngs_modulus: float, poisson_ratio: float
) -> npt.NDArray[np.floating[Any]]

Construct the isotropic axisymmetric elastic stress-strain matrix.

Parameters:

Name Type Description Default
youngs_modulus float

Young's modulus with units [pressure].

 required
poisson_ratio float

Poisson ratio with units [dimensionless].

 required

Returns:

Name Type Description
NDArray NDArray[floating[Any]]

Elastic stress-strain matrix with shape (4, 4) in component order
NDArray[floating[Any]]

[rr, zz, tt, rz]. Units are [stress / strain] = [pressure].
Source code in cfsem/solenoid_stress/fem2d.py 
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
def isotropic_axisymmetric_material(
    youngs_modulus: float,
    poisson_ratio: float,
) -> npt.NDArray[np.floating[Any]]:
    """Construct the isotropic axisymmetric elastic stress-strain matrix.

    Args:
        youngs_modulus: Young's modulus with units `[pressure]`.
        poisson_ratio: Poisson ratio with units `[dimensionless]`.

    Returns:
        NDArray: Elastic stress-strain matrix with shape `(4, 4)` in component order
        `[rr, zz, tt, rz]`. Units are `[stress / strain] = [pressure]`.
    """

    return _as_float64_array(_isotropic_axisymmetric_material_f64(youngs_modulus, poisson_ratio)).reshape(
        4, 4
    )

isotropic_axisymmetric_thermal_material 

isotropic_axisymmetric_thermal_material(
    alpha: float, reference_temperature: float = 0.0
) -> npt.NDArray[np.floating[Any]]

Construct isotropic thermal-expansion data.

Parameters:

Name Type Description Default
alpha float

Isotropic thermal expansion coefficient with units [strain / temperature].

 required
reference_temperature float

Stress-free reference temperature with units [temperature].

 0.0

Returns:

Name Type Description
NDArray NDArray[floating[Any]]

Thermal material row with shape (5,) storing
NDArray[floating[Any]]

[alpha_r, alpha_z, alpha_t, alpha_rz, T_ref]. The first four entries have units
NDArray[floating[Any]]

[strain / temperature]; T_ref has units [temperature].
Source code in cfsem/solenoid_stress/fem2d.py 
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
def isotropic_axisymmetric_thermal_material(
    alpha: float,
    reference_temperature: float = 0.0,
) -> npt.NDArray[np.floating[Any]]:
    """Construct isotropic thermal-expansion data.

    Args:
        alpha: Isotropic thermal expansion coefficient with units `[strain / temperature]`.
        reference_temperature: Stress-free reference temperature with units `[temperature]`.

    Returns:
        NDArray: Thermal material row with shape `(5,)` storing
        `[alpha_r, alpha_z, alpha_t, alpha_rz, T_ref]`. The first four entries have units
        `[strain / temperature]`; `T_ref` has units `[temperature]`.
    """

    return _as_float64_array(_isotropic_axisymmetric_thermal_material_f64(alpha, reference_temperature))

isotropic_plane_strain_material 

isotropic_plane_strain_material(
    youngs_modulus: float, poisson_ratio: float
) -> npt.NDArray[np.floating[Any]]

Construct the isotropic plane-strain elastic stress-strain matrix.

Returns a dense (4, 4) constitutive matrix in [xx, yy, zz, xy] order. The plane-strain solver sets epsilon_zz = 0, but this matrix still recovers the nonzero sigma_zz implied by the in-plane strains.

Source code in cfsem/solenoid_stress/fem2d.py 
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
def isotropic_plane_strain_material(
    youngs_modulus: float,
    poisson_ratio: float,
) -> npt.NDArray[np.floating[Any]]:
    """Construct the isotropic plane-strain elastic stress-strain matrix.

    Returns a dense `(4, 4)` constitutive matrix in `[xx, yy, zz, xy]` order. The plane-strain
    solver sets `epsilon_zz = 0`, but this matrix still recovers the nonzero `sigma_zz` implied
    by the in-plane strains.
    """

    return _as_float64_array(_isotropic_plane_strain_material_f64(youngs_modulus, poisson_ratio)).reshape(
        4, 4
    )

isotropic_plane_strain_thermal_material 

isotropic_plane_strain_thermal_material(
    alpha: float, reference_temperature: float = 0.0
) -> npt.NDArray[np.floating[Any]]

Construct isotropic plane-strain thermal-expansion data.

Returns a row [alpha_x, alpha_y, alpha_z, alpha_xy, T_ref] with equal normal expansion coefficients and zero engineering shear expansion.

Source code in cfsem/solenoid_stress/fem2d.py 
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
def isotropic_plane_strain_thermal_material(
    alpha: float,
    reference_temperature: float = 0.0,
) -> npt.NDArray[np.floating[Any]]:
    """Construct isotropic plane-strain thermal-expansion data.

    Returns a row `[alpha_x, alpha_y, alpha_z, alpha_xy, T_ref]` with equal normal expansion
    coefficients and zero engineering shear expansion.
    """

    return _as_float64_array(_isotropic_plane_strain_thermal_material_f64(alpha, reference_temperature))

orthotropic_axisymmetric_thermal_material 

orthotropic_axisymmetric_thermal_material(
    alpha_r: float,
    alpha_z: float,
    alpha_t: float,
    reference_temperature: float = 0.0,
) -> npt.NDArray[np.floating[Any]]

Construct orthotropic thermal-expansion data.

Parameters:

Name Type Description Default
alpha_r float

Radial thermal expansion coefficient with units [strain / temperature].

 required
alpha_z float

Axial thermal expansion coefficient with units [strain / temperature].

 required
alpha_t float

Hoop thermal expansion coefficient with units [strain / temperature].

 required
reference_temperature float

Stress-free reference temperature with units [temperature].

 0.0

Returns:

Name Type Description
NDArray NDArray[floating[Any]]

Thermal material row with shape (5,) storing
NDArray[floating[Any]]

[alpha_r, alpha_z, alpha_t, alpha_rz, T_ref]. The first four entries have units
NDArray[floating[Any]]

[strain / temperature]; T_ref has units [temperature].
Source code in cfsem/solenoid_stress/fem2d.py 
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
def orthotropic_axisymmetric_thermal_material(
    alpha_r: float,
    alpha_z: float,
    alpha_t: float,
    reference_temperature: float = 0.0,
) -> npt.NDArray[np.floating[Any]]:
    """Construct orthotropic thermal-expansion data.

    Args:
        alpha_r: Radial thermal expansion coefficient with units `[strain / temperature]`.
        alpha_z: Axial thermal expansion coefficient with units `[strain / temperature]`.
        alpha_t: Hoop thermal expansion coefficient with units `[strain / temperature]`.
        reference_temperature: Stress-free reference temperature with units `[temperature]`.

    Returns:
        NDArray: Thermal material row with shape `(5,)` storing
        `[alpha_r, alpha_z, alpha_t, alpha_rz, T_ref]`. The first four entries have units
        `[strain / temperature]`; `T_ref` has units `[temperature]`.
    """

    return _as_float64_array(
        _orthotropic_axisymmetric_thermal_material_f64(
            alpha_r,
            alpha_z,
            alpha_t,
            reference_temperature,
        )
    )

orthotropic_plane_strain_thermal_material 

orthotropic_plane_strain_thermal_material(
    alpha_x: float,
    alpha_y: float,
    alpha_z: float,
    reference_temperature: float = 0.0,
) -> npt.NDArray[np.floating[Any]]

Construct orthotropic plane-strain thermal-expansion data.

Returns a row [alpha_x, alpha_y, alpha_z, alpha_xy, T_ref] with zero engineering shear expansion. Use material_orientation_angles during assembly to rotate local orthotropic axes.

Source code in cfsem/solenoid_stress/fem2d.py 
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
def orthotropic_plane_strain_thermal_material(
    alpha_x: float,
    alpha_y: float,
    alpha_z: float,
    reference_temperature: float = 0.0,
) -> npt.NDArray[np.floating[Any]]:
    """Construct orthotropic plane-strain thermal-expansion data.

    Returns a row `[alpha_x, alpha_y, alpha_z, alpha_xy, T_ref]` with zero engineering shear
    expansion. Use `material_orientation_angles` during assembly to rotate local orthotropic axes.
    """

    return orthotropic_axisymmetric_thermal_material(
        alpha_x,
        alpha_y,
        alpha_z,
        reference_temperature,
    )

pack_material_tables_from_tags 

pack_material_tables_from_tags(
    material_ids: ArrayLike,
    material_table_by_tag: Mapping[int, ArrayLike],
    thermal_material_table_by_tag: Mapping[int, ArrayLike]
    | None = None,
) -> tuple[
    npt.NDArray[np.uint64],
    npt.NDArray[np.floating[Any]],
    npt.NDArray[np.floating[Any]] | None,
]

Pack tagged material definitions into the dense FEM input format.

Parameters:

Name Type Description Default
material_ids ArrayLike

Element material tags with shape (nelem,).

 required
material_table_by_tag Mapping[int, ArrayLike]

Mapping from external material tag to elastic stress-strain matrix with shape (4, 4). Matrix units are [stress / strain] = [pressure].

 required
thermal_material_table_by_tag Mapping[int, ArrayLike] | None

Optional mapping from external material tag to thermal row with shape (5,) storing [alpha_r, alpha_z, alpha_t, alpha_rz, T_ref]. Thermal expansion coefficients have units [strain / temperature] and T_ref has units [temperature].

 None

Returns: tuple: (packed_material_ids, packed_material_table, packed_thermal_material_table) where: packed_material_ids has shape (nelem,), packed_material_table has shape (nmat, 4, 4), packed_thermal_material_table has shape (nmat, 5) when provided, otherwise None.

Raises:

Type Description
ValueError

If an element tag is missing from material_table_by_tag, or if thermal tags do not match the elastic tags exactly.
Source code in cfsem/solenoid_stress/fem2d.py 
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
def pack_material_tables_from_tags(
    material_ids: ArrayLike,
    material_table_by_tag: Mapping[int, ArrayLike],
    thermal_material_table_by_tag: Mapping[int, ArrayLike] | None = None,
) -> tuple[
    npt.NDArray[np.uint64],
    npt.NDArray[np.floating[Any]],
    npt.NDArray[np.floating[Any]] | None,
]:
    """Pack tagged material definitions into the dense FEM input format.

    Args:
        material_ids: Element material tags with shape `(nelem,)`.
        material_table_by_tag: Mapping from external material tag to elastic stress-strain matrix
            with shape `(4, 4)`. Matrix units are `[stress / strain] = [pressure]`.
        thermal_material_table_by_tag: Optional mapping from external material tag to thermal row
            with shape `(5,)` storing `[alpha_r, alpha_z, alpha_t, alpha_rz, T_ref]`. Thermal
            expansion coefficients have units `[strain / temperature]` and `T_ref` has units
            `[temperature]`.
    Returns:
        tuple: `(packed_material_ids, packed_material_table, packed_thermal_material_table)` where:
            `packed_material_ids` has shape `(nelem,)`,
            `packed_material_table` has shape `(nmat, 4, 4)`,
            `packed_thermal_material_table` has shape `(nmat, 5)` when provided, otherwise `None`.

    Raises:
        ValueError: If an element tag is missing from `material_table_by_tag`, or if thermal tags
            do not match the elastic tags exactly.
    """

    assert material_table_by_tag, "material_table_by_tag cannot be empty"
    resolved_dtype = np.dtype(np.float64)
    ids = np.asarray(material_ids, dtype=np.uint64)
    assert ids.ndim == 1, f"material_ids must have shape (nelem,); got {ids.shape}"

    material_tags = sorted(int(tag) for tag in material_table_by_tag)
    tag_to_index = {tag: index for index, tag in enumerate(material_tags)}
    try:
        packed_ids = np.asarray([tag_to_index[int(tag)] for tag in ids], dtype=np.uint64)
    except KeyError as exc:
        raise ValueError(
            f"material_ids contains tag {exc.args[0]} that is missing from material_table_by_tag"
        ) from exc

    material_rows = []
    for tag in material_tags:
        matrix = np.asarray(material_table_by_tag[tag], dtype=resolved_dtype)
        assert matrix.shape == (
            4,
            4,
        ), f"material_table_by_tag[{tag}] must have shape (4, 4); got {matrix.shape}"
        material_rows.append(matrix)
    packed_material_table = np.ascontiguousarray(np.stack(material_rows, axis=0), dtype=resolved_dtype)

    packed_thermal_table: npt.NDArray[np.floating[Any]] | None
    if thermal_material_table_by_tag is None:
        packed_thermal_table = None
    else:
        thermal_tags = {int(tag) for tag in thermal_material_table_by_tag}
        if thermal_tags != set(material_tags):
            raise ValueError(
                "thermal_material_table_by_tag must have exactly the same keys as material_table_by_tag"
            )
        thermal_rows = []
        for tag in material_tags:
            row = np.asarray(thermal_material_table_by_tag[tag], dtype=resolved_dtype)
            assert row.shape == (
                5,
            ), f"thermal_material_table_by_tag[{tag}] must have shape (5,); got {row.shape}"
            thermal_rows.append(row)
        packed_thermal_table = np.ascontiguousarray(np.stack(thermal_rows, axis=0), dtype=resolved_dtype)
        assert not np.any(
            packed_thermal_table[:, 3] != 0.0
        ), "shear thermal expansion (alpha_rz) is not yet supported"

    return packed_ids, packed_material_table, packed_thermal_table

quad_mesh_interpolation_operator 

quad_mesh_interpolation_operator(
    query: QuadMeshQuery,
) -> sp.csr_matrix

Return a reusable sparse operator mapping nodal scalar values to query-point values.

The returned matrix has shape (npoint, nnode). Applying it to a dense (nnode,) vector gives scalar values at the query points; applying it to (nnode, ncomponent) interpolates multiple nodal fields with the same operator. The operator always uses the query's nearest element.

Parameters:

Name Type Description Default
query QuadMeshQuery

Mesh query data from query_quad_mesh(...).

 required

Returns:

Type Description
csr_matrix

Sparse interpolation operator with shape (npoint, nnode). Entries are unitless shape
csr_matrix

function values, so output values have the same units as the nodal values supplied during
csr_matrix

matrix multiplication.
Source code in cfsem/solenoid_stress/fem2d.py 
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
def quad_mesh_interpolation_operator(
    query: QuadMeshQuery,
) -> sp.csr_matrix:
    """Return a reusable sparse operator mapping nodal scalar values to query-point values.

    The returned matrix has shape `(npoint, nnode)`. Applying it to a dense `(nnode,)` vector gives
    scalar values at the query points; applying it to `(nnode, ncomponent)` interpolates multiple
    nodal fields with the same operator. The operator always uses the query's nearest element.

    Args:
        query: Mesh query data from `query_quad_mesh(...)`.

    Returns:
        Sparse interpolation operator with shape `(npoint, nnode)`. Entries are unitless shape
        function values, so output values have the same units as the nodal values supplied during
        matrix multiplication.
    """

    return _coo_operator_from_binding(
        _quad_mesh_interpolation_operator_f64(
            query.nodes,
            query.elements,
            np.asarray(query.nearest_element_indices, dtype=np.uint64),
            query.nearest_element_reference_points,
            query.element_type,
        ),
    )

query_quad_mesh 

query_quad_mesh(
    nodes: ArrayLike,
    elements: ArrayLike,
    points: ArrayLike,
    *,
    element_type: str = "quad4",
    max_iterations: int = 20,
) -> QuadMeshQuery

Query nearest node, nearest element, and nearest face in one pass.

The Rust backend scans all nodes once and all elements once per query point. The element scan computes nearest-element and nearest-face metadata together so downstream interpolation and recovery operators do not repeat point location. Complexity is O(npoint * (nnode + nelem * max_iterations)). A point is contained when its nearest-element distance is zero to the caller's tolerance.

Parameters:

Name Type Description Default
nodes ArrayLike

Mesh node coordinates with shape (nnode, 2) and units [length].

 required
elements ArrayLike

Quad connectivity with shape (nelem, 4) for quad4 or (nelem, 9) for quad9. Entries are unitless node indices.

 required
points ArrayLike

Query point coordinates with shape (npoint, 2) and units [length].

 required
element_type str

Element family, either "quad4" or "quad9".

 'quad4'
max_iterations int

Maximum Newton/projection iterations per element. Unitless.

 20

Returns:

Type Description
QuadMeshQuery

Query data with nearest-node, nearest-element, and nearest-face arrays. Coordinate arrays
QuadMeshQuery

have units [length], distances have units [length], reference coordinates are unitless,
QuadMeshQuery

and index arrays are unitless.
Source code in cfsem/solenoid_stress/fem2d.py 
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
def query_quad_mesh(
    nodes: ArrayLike,
    elements: ArrayLike,
    points: ArrayLike,
    *,
    element_type: str = "quad4",
    max_iterations: int = 20,
) -> QuadMeshQuery:
    """Query nearest node, nearest element, and nearest face in one pass.

    The Rust backend scans all nodes once and all elements once per query point. The element scan
    computes nearest-element and nearest-face metadata together so downstream interpolation and
    recovery operators do not repeat point location. Complexity is
    `O(npoint * (nnode + nelem * max_iterations))`. A point is contained when its
    nearest-element distance is zero to the caller's tolerance.

    Args:
        nodes: Mesh node coordinates with shape `(nnode, 2)` and units `[length]`.
        elements: Quad connectivity with shape `(nelem, 4)` for `quad4` or `(nelem, 9)` for
            `quad9`. Entries are unitless node indices.
        points: Query point coordinates with shape `(npoint, 2)` and units `[length]`.
        element_type: Element family, either `"quad4"` or `"quad9"`.
        max_iterations: Maximum Newton/projection iterations per element. Unitless.

    Returns:
        Query data with nearest-node, nearest-element, and nearest-face arrays. Coordinate arrays
        have units `[length]`, distances have units `[length]`, reference coordinates are unitless,
        and index arrays are unitless.
    """

    normalized_element_type = _normalize_element_type(element_type)
    nodes_arr = _normalize_nodes(nodes)
    elements_arr = _normalize_elements(
        elements,
        4 if normalized_element_type == "quad4" else 9,
    )
    points_arr = _normalize_query_points(points)
    data = _quad_mesh_query_f64(
        nodes_arr,
        elements_arr,
        points_arr,
        normalized_element_type,
        int(max_iterations),
    )

    return QuadMeshQuery(
        nodes=nodes_arr,
        elements=elements_arr,
        points=points_arr,
        element_type=normalized_element_type,
        nearest_node_indices=np.asarray(data["nearest_node_indices"], dtype=np.int64),
        nearest_node_points=np.asarray(data["nearest_node_points"], dtype=np.float64).reshape(-1, 2),
        nearest_node_distances=np.asarray(data["nearest_node_distances"], dtype=np.float64),
        nearest_element_indices=np.asarray(data["nearest_element_indices"], dtype=np.int64),
        nearest_element_reference_points=np.asarray(
            data["nearest_element_reference_points"],
            dtype=np.float64,
        ).reshape(-1, 2),
        nearest_element_points=np.asarray(data["nearest_element_points"], dtype=np.float64).reshape(-1, 2),
        nearest_element_distances=np.asarray(data["nearest_element_distances"], dtype=np.float64),
        nearest_face_element_indices=np.asarray(data["nearest_face_element_indices"], dtype=np.int64),
        nearest_face_local_faces=np.asarray(data["nearest_face_local_faces"], dtype=np.int64),
        nearest_face_reference_coordinates=np.asarray(
            data["nearest_face_reference_coordinates"],
            dtype=np.float64,
        ),
        nearest_face_points=np.asarray(data["nearest_face_points"], dtype=np.float64).reshape(-1, 2),
        nearest_face_distances=np.asarray(data["nearest_face_distances"], dtype=np.float64),
    )

Analytic Reference Formulas

cfsem.solenoid_stress.s_long_solenoid 

s_long_solenoid(
    r: NDArray,
    ri: float,
    ro: float,
    j: float,
    bzi: float,
    bzo: float,
    poisson_ratio: float,
) -> tuple[NDArray, NDArray]

Radial and hoop stress in an infinitely long solenoid under linearly-varying self field. The "infinite length" assumption is equivalent to assuming zero R-Z shear ("deck of cards") and assuming no B-field in the R-direction (no Z-load or stress).

The linearly varying B-field from inside to outside allows slightly extending this to partially account for the finite length of a real solenoid, which produces a region of negative Bz near the outer radius (as opposed to the true infinite solenoid, for which Bz trends to exactly zero at the outer radius).

Iwasa 2e pg 101 eqns 3.77a,b .

Assumes * Infinitely long solenoid (no R-field or Z-stress). * Linear B-field fall-off between inner and outer radius * "Very long" solenoid - allows some negative field at the OD, but always linearly varying * Uniform current density; no bulk regions of non-conducting structure * Radial stress at inner and outer radius is zero (BC due to no support) * Isotropic material * No thermal stress

Can acommodate a uniform or linearly-varying background field, but not general fields.

Parameters:

Name Type Description Default
r NDArray

[m] (n x 1) array of radius points at which to evaluate the stress

 required
ri float

[m] inner radius

 required
ro float

[m] outer radius

 required
j float

[A/m^2] current density

 required
bzi float

[T] axial B-field at inner radius

 required
bzo float

[T] axial B-field at outer radius

 required
poisson_ratio float

[dimensionless] Material property; off-axis stress coupling term

 required

Returns:

Type Description
tuple[NDArray, NDArray]

s_radial, s_hoop - each (n x 1) with units of [Pa]
Source code in cfsem/solenoid_stress/solenoid_handcalc.py 
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
def s_long_solenoid(
    r: NDArray,
    ri: float,
    ro: float,
    j: float,
    bzi: float,
    bzo: float,
    poisson_ratio: float,
) -> tuple[NDArray, NDArray]:
    """
    Radial and hoop stress in an infinitely long solenoid under linearly-varying self field.
    The "infinite length" assumption is equivalent to assuming zero R-Z shear ("deck of cards")
    and assuming no B-field in the R-direction (no Z-load or stress).

    The linearly varying B-field from inside to outside allows slightly extending this
    to partially account for the finite length of a real solenoid, which produces
    a region of negative Bz near the outer radius (as opposed to the true infinite solenoid,
    for which Bz trends to exactly zero at the outer radius).

    Iwasa 2e pg 101 eqns 3.77a,b .

    Assumes
    * Infinitely long solenoid (no R-field or Z-stress).
    * Linear B-field fall-off between inner and outer radius
      * "Very long" solenoid - allows some negative field at the OD, but always linearly varying
    * Uniform current density; no bulk regions of non-conducting structure
    * Radial stress at inner and outer radius is zero (BC due to no support)
    * Isotropic material
    * No thermal stress

    Can acommodate a uniform or linearly-varying background field, but not general fields.

    Args:
        r: [m] (n x 1) array of radius points at which to evaluate the stress
        ri: [m] inner radius
        ro: [m] outer radius
        j: [A/m^2] current density
        bzi: [T] axial B-field at inner radius
        bzo: [T] axial B-field at outer radius
        poisson_ratio: [dimensionless] Material property; off-axis stress coupling term

    Returns:
        s_radial, s_hoop - each (n x 1) with units of [Pa]
    """

    # Terms shared between s_radial and s_hoop
    nu = poisson_ratio
    rho = r / ri
    alpha = ro / ri
    kappa = bzo / bzi

    jbr = j * bzi * ri  # [Pa]
    term1 = jbr / (alpha - 1.0)  # [Pa]
    term2 = (2.0 + nu) / 3.0
    term3 = (3.0 + nu) / 8.0
    term4 = alpha - kappa
    term5 = 1.0 - kappa

    # s_radial
    term6 = ((alpha**2 + alpha + 1.0 - alpha**2 / rho**2) / (alpha + 1.0)) - rho
    term7 = term3 * term5 * (alpha**2 + 1.0 - alpha**2 / rho**2 - rho**2)
    s_radial = term1 * (term2 * term4 * term6 - term7)  # [Pa]

    # s_hoop
    term8 = term2 * (alpha**2 + alpha + 1.0 + alpha**2 / rho**2) / (alpha + 1.0)
    term9 = rho * (1.0 + 2.0 * nu) / 3.0
    term10 = term4 * (term8 - term9)

    term11 = term3 * (alpha**2 + 1.0 + alpha**2 / rho**2)
    term12 = rho**2 * (1.0 + 3.0 * nu) / 8.0
    term13 = term5 * (term11 - term12)

    s_hoop = term1 * (term10 - term13)  # [Pa]

    return s_radial, s_hoop  # [Pa]

cfsem.solenoid_stress.s_radial_thick_wall_cylinder 

s_radial_thick_wall_cylinder(
    r: NDArray,
    ri: float,
    ro: float,
    pin: float,
    pout: float,
) -> NDArray

Radial stress at a location in a thick walled cylinder under pressure load with ends "capped", although the capped constraint does not affect the hoop or radial stress compared to an infinite-length constraint.

https://www.engineeringtoolbox.com/stress-thick-walled-tube-d_949.html https://www.suncam.com/miva/downloads/docs/303.pdf

Parameters:

Name Type Description Default
r NDArray

[m] radius at which to evaluate

 required
ri float

[m] inner radius

 required
ro float

[m] outer radius

 required
pin float

[Pa] inside pressure

 required
pout float

[Pa] outside pressure

 required

Returns:

Type Description
NDArray

[Pa] radial stress
Source code in cfsem/solenoid_stress/thick_wall_cylinder_handcalc.py 
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
def s_radial_thick_wall_cylinder(r: NDArray, ri: float, ro: float, pin: float, pout: float) -> NDArray:
    """
    Radial stress at a location in a thick walled cylinder under pressure load
    with ends "capped", although the capped constraint does not affect the hoop or radial stress
    compared to an infinite-length constraint.

    https://www.engineeringtoolbox.com/stress-thick-walled-tube-d_949.html
    https://www.suncam.com/miva/downloads/docs/303.pdf

    Args:
        r: [m] radius at which to evaluate
        ri: [m] inner radius
        ro: [m] outer radius
        pin: [Pa] inside pressure
        pout: [Pa] outside pressure

    Returns:
        [Pa] radial stress
    """
    # Factors of pi cancel out
    # fmt: off
    s_radial = (
        ((pin * ri**2 - pout * ro**2) / (ro**2 - ri**2)) + \
        (ri**2 * ro**2 * (pout - pin) / (r**2 * (ro**2 - ri**2)))
    )
    # fmt: on

    return s_radial  # [Pa] radial stress

cfsem.solenoid_stress.s_hoop_thick_wall_cylinder 

s_hoop_thick_wall_cylinder(
    r: NDArray,
    ri: float,
    ro: float,
    pin: float,
    pout: float,
) -> NDArray

Hoop stress at a location in a thick walled cylinder under pressure load with ends "capped", although the capped constraint does not affect the hoop or radial stress compared to an infinite-length constraint.

https://www.engineeringtoolbox.com/stress-thick-walled-tube-d_949.html https://www.suncam.com/miva/downloads/docs/303.pdf

Parameters:

Name Type Description Default
r NDArray

[m] radius at which to evaluate

 required
ri float

[m] inner radius

 required
ro float

[m] outer radius

 required
pin float

[Pa] inside pressure

 required
pout float

[Pa] outside pressure

 required

Returns:

Type Description
NDArray

[Pa] hoop stress
Source code in cfsem/solenoid_stress/thick_wall_cylinder_handcalc.py 
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
def s_hoop_thick_wall_cylinder(r: NDArray, ri: float, ro: float, pin: float, pout: float) -> NDArray:
    """
    Hoop stress at a location in a thick walled cylinder under pressure load
    with ends "capped", although the capped constraint does not affect the hoop or radial stress
    compared to an infinite-length constraint.

    https://www.engineeringtoolbox.com/stress-thick-walled-tube-d_949.html
    https://www.suncam.com/miva/downloads/docs/303.pdf

    Args:
        r: [m] radius at which to evaluate
        ri: [m] inner radius
        ro: [m] outer radius
        pin: [Pa] inside pressure
        pout: [Pa] outside pressure

    Returns:
        [Pa] hoop stress
    """
    # Factors of pi cancel out
    # fmt: off
    s_hoop = (
        (pin * ri ** 2 - pout * ro ** 2) / (ro ** 2 - ri ** 2) -
        (ri ** 2 * ro ** 2 * (pout - pin) / (r ** 2 * (ro ** 2 - ri ** 2)))
    )
    # fmt: on

    return s_hoop  # [Pa] hoop stress