Skip to content

models.grain.fmm

Fast Marching Method (FMM) base classes for grain geometries whose cross-section cannot be described analytically. The burning surface at any web distance is read off a single distance map of the initial port, solved with scikit-fmm.

  • FMMGrainSegment — Shared regression logic: builds the distance map and derives web thickness, the regressed face map, and port and burn area from it.
  • FMMGrainSegment2D — Constant cross-section geometries (Star, D-grain, wagon-wheel, multi-port, rod-and-tube). Burn area is the core perimeter (a single contour of the distance map) times the grain length, plus any exposed end faces.
  • FMMGrainSegment3D — Axially varying cross-section geometries (Conical, Finocyl). The port is a stack of axial slices through a 3D distance map; burn area is the marching-cubes area of the regressing iso-surface, and port area is read from a chosen slice.
  • FMMSTLGrainSegment — Arbitrary grain geometry imported from an STL mesh file.

Used internally by the concrete geometry classes in geometries.

machwave.models.grain.fmm

Fast Marching Method (FMM) grain geometry models.

Supports both 2D and 3D geometries, as well as STL meshes.

References

https://math.berkeley.edu/~sethian/2006/Explanations/fast_marching_explain.html

FMMGrainSegment

Bases: GrainSegment, ABC

Fast Marching Method (FMM) implementation of a grain segment.

Source code in machwave/models/grain/fmm/base.py 
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
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
class FMMGrainSegment(grain.GrainSegment, ABC):
    """Fast Marching Method (FMM) implementation of a grain segment."""

    def __init__(
        self,
        length: float,
        outer_diameter: float,
        inhibited_surfaces: grain_base.InhibitedSurfaces | None = None,
        grid_resolution: int = DEFAULT_GRID_RESOLUTION,
        density_ratio: float = 1.0,
    ) -> None:
        """
        Initialize an FMM grain segment.

        Args:
            length: Segment length [m].
            outer_diameter: Outer diameter [m].
            inhibited_surfaces: Surfaces inhibited from burning.
            grid_resolution: Resolution of the face map, x and y axes.
            density_ratio: Ratio of real to ideal propellant density.
        """
        self.grid_resolution = grid_resolution

        # Cache variables:
        self.coordinate_grids = None
        self.outer_diameter_mask = None
        self.masked_face = None
        self.regression_map = None
        self.web_thickness = None

        # Cache variables that only work for a specific web distance:
        self._solid_mask_web = None
        self._solid_mask = None
        self._solid_indices_web = None
        self._solid_indices = None

        super().__init__(
            length=length,
            outer_diameter=outer_diameter,
            inhibited_surfaces=inhibited_surfaces,
            density_ratio=density_ratio,
        )

    def validate(self) -> None:
        """
        Validate the internal geometry of the grain.

        Raises:
            GrainGeometryError: If the grid resolution is below the valid
                threshold.
        """
        super().validate()

        if not self.grid_resolution >= MINIMUM_GRID_RESOLUTION:
            raise grain.GrainGeometryError(
                f"Grid resolution must be at least {MINIMUM_GRID_RESOLUTION}, "
                f"got {self.grid_resolution}"
            )

    @abstractmethod
    def get_coordinate_grids(self) -> tuple:
        """Return the coordinate grids for the grain, one per spatial axis."""
        pass

    @abstractmethod
    def get_outer_diameter_mask(self) -> np.ndarray:
        """Return a boolean mask of cells outside the outer-diameter boundary."""
        pass

    @abstractmethod
    def generate_initial_face_map(self) -> np.typing.NDArray[np.int_]:
        """Generate the initial face map for the geometry (1 = propellant, 0 = void)."""
        pass

    def get_empty_face_map(self) -> np.ndarray:
        """Return a face map of all ones, shaped like the first stored map."""
        return np.ones_like(self.get_coordinate_grids()[0])

    def normalize(self, value: int | float) -> float:
        """
        Convert a dimensional value to a fraction of the outer radius.

        Args:
            value: Dimensional value (e.g., length) to normalize [m].

        Returns:
            Value expressed as a fraction of the outer radius.
        """
        return value / (0.5 * self.outer_diameter)

    def denormalize(
        self, value: int | float | NDArray[np.float64]
    ) -> float | NDArray[np.float64]:
        """Convert a normalized input back into a dimensional value [m]."""
        return (value / 2) * (self.outer_diameter)

    def get_radial_grid_spacing(self) -> float:
        """
        Return the distance between adjacent cross-section cells [m].

        The grid spans the outer diameter over `grid_resolution` samples, so it
        holds `grid_resolution - 1` intervals.
        """
        return self.outer_diameter / (self.grid_resolution - 1)

    def cells_to_meters(
        self, value: float | NDArray[np.float64]
    ) -> float | NDArray[np.float64]:
        """
        Convert a cell-distance value to [m].

        Args:
            value: Distance in cell units.

        Returns:
            Distance [m].
        """
        return value * self.get_radial_grid_spacing()

    def cells_to_square_meters(
        self, value: float | NDArray[np.float64]
    ) -> float | NDArray[np.float64]:
        """
        Convert a cell-area value to [m^2].

        Args:
            value: Area in cell units.

        Returns:
            Area [m^2].
        """
        return value * self.get_radial_grid_spacing() ** 2

    def _apply_surface_inhibition(
        self,
        face_map: NDArray[np.int_],
        excluded_mask: NDArray[np.bool_],
    ) -> tuple[NDArray[np.int_], NDArray[np.bool_]]:
        """
        Apply surface inhibition maps.

        Base implementation handles outer surface inhibition. 2D and 3D subclasses add
        their own logic for other inhibited surfaces.

        Args:
            face_map: Mutable copy of the initial face map.
            excluded_mask: Mutable copy of the outer diameter mask; cells excluded from
                the burn domain, extended here with inhibited surfaces.

        Returns:
            `(face_map, excluded_mask)` with inhibition applied.
        """
        if not self.inhibited_surfaces.outer_surface:
            inside = ~excluded_mask  # Invert the mask
            eroded = binary_erosion(inside)  # Erode the inside to find the boundary
            outer_surface_boundary = inside & np.logical_not(eroded)
            face_map[outer_surface_boundary] = 0

        return face_map, excluded_mask

    def get_masked_face(self) -> np.ndarray:
        """Return a masked representation of the face map."""
        if self.masked_face is None:
            face_map = self.generate_initial_face_map().copy()
            excluded_mask = self.get_outer_diameter_mask().copy()

            face_map, excluded_mask = self._apply_surface_inhibition(
                face_map, excluded_mask
            )

            self.masked_face = np.ma.MaskedArray(face_map, excluded_mask)
        return self.masked_face

    @property
    def has_cross_section_regression(self) -> bool:
        """Return whether the cross-section has any burning surface."""
        masked_face = self.get_masked_face()
        unmasked = ~np.ma.getmaskarray(masked_face)
        return bool(np.any(masked_face.data[np.asarray(unmasked, dtype=bool)] == 0))

    def _compute_regression_distance(self, masked_face: np.ndarray) -> np.ndarray:
        """Return the regression distance from the burning surface in normalized web units."""
        distance = skfmm.distance(masked_face, dx=self.get_radial_grid_spacing())
        return distance * (2.0 / self.outer_diameter)

    def get_regression_map(self):
        """
        Return the distance map for grain regression.

        Uses the fast marching method (scikit-fmm) on the masked face. Each
        value is the distance from the initial face along the cross-section.
        With no burning surface (end burner), returns a static map set to the
        axial burnout web, so the grain regresses along its length, not a cell.
        """
        if self.regression_map is None:
            masked_face = self.get_masked_face()

            if self.has_cross_section_regression:
                self.regression_map = self._compute_regression_distance(masked_face)
            else:
                # End burner: web = length split across exposed ends. Fill the
                # whole array with one constant (not 0 outside the mask) so the
                # perimeter tracer finds no spurious contour.
                exposed_end_count = (not self.inhibited_surfaces.upper_end) + (
                    not self.inhibited_surfaces.lower_end
                )
                normalized_axial_web_distance = (
                    self.normalize(self.length / exposed_end_count)
                    if exposed_end_count
                    else 0.0
                )
                self.regression_map = np.ma.MaskedArray(
                    np.full(masked_face.shape, normalized_axial_web_distance),
                    mask=np.ma.getmaskarray(masked_face),
                )
        return self.regression_map

    def get_web_thickness(self) -> float:
        """
        Return the maximum web thickness of the grain [m].

        Derived from the distance map as the largest normalized distance,
        converted back into real units.
        """
        if self.web_thickness is None:
            self.web_thickness = float(
                self.denormalize(float(np.amax(self.get_regression_map())))
            )
        return float(self.web_thickness)

    def get_face_map(self, web_distance: float) -> np.typing.NDArray[np.int64]:
        """
        Returns a matrix representing the grain face based on the given web distance.

        The returned array can contain:
        -1 for masked or invalid points,
        0 for points below the threshold,
        1 for points above the threshold.

        Args:
            web_distance: The distance traveled into the grain web.

        Returns:
            A NumPy array with -1, 0, or 1 indicating the grain face at the specified web distance.
        """
        web_distance_normalized = self.normalize(web_distance)
        regression_map = self.get_regression_map()
        excluded_mask = np.ma.getmaskarray(regression_map)

        # Create a masked array, where excluded cells are masked out
        occupancy_state = np.ma.MaskedArray(
            (regression_map > web_distance_normalized).astype(np.int64),
            mask=excluded_mask,
        )

        # Fill masked entries with -1, valid/true entries remain 1 or 0
        return occupancy_state.filled(-1)

    def _get_solid_mask(self, web_distance: float) -> NDArray[np.bool_]:
        """
        Boolean mask of solid cells at a web distance.

        Memoized per web distance, since it is used for volume, center of gravity
        and moment of inertia calculations.
        """
        if self._solid_mask is None or self._solid_mask_web != web_distance:
            regression_map = self.get_regression_map()
            web_distance_normalized = self.normalize(web_distance)
            excluded_mask = np.ma.getmaskarray(regression_map)
            self._solid_mask = (
                np.ma.getdata(regression_map) > web_distance_normalized
            ) & ~excluded_mask
            self._solid_mask_web = web_distance
        return self._solid_mask

    def _get_solid_indices(self, web_distance: float) -> tuple[NDArray[np.int_], ...]:
        """Indices of solid cells, memoized."""
        if self._solid_indices is None or self._solid_indices_web != web_distance:
            self._solid_indices = np.where(self._get_solid_mask(web_distance))
            self._solid_indices_web = web_distance
        return self._solid_indices

    @abstractmethod
    def get_contours(
        self, web_distance: float, *args: Any, **kwargs: Any
    ) -> list[NDArray[np.float64]]:
        """
        Return the contours of the regression map for a given web distance.

        Must be implemented by a subclass to compute the contour data based on
        the given web distance and any additional parameters.

        Args:
            web_distance: Depth of regression into the grain web [m].
            *args: Additional positional arguments.
            **kwargs: Additional keyword arguments.

        Returns:
            An array representing the computed contours of the grain regression.
        """
        pass

has_cross_section_regression property

Return whether the cross-section has any burning surface.

__init__(length, outer_diameter, inhibited_surfaces=None, grid_resolution=DEFAULT_GRID_RESOLUTION, density_ratio=1.0)

Initialize an FMM grain segment.

Parameters:

Name Type Description Default
length float

Segment length [m].

 required
outer_diameter float

Outer diameter [m].

 required
inhibited_surfaces InhibitedSurfaces | None

Surfaces inhibited from burning.

 None
grid_resolution int

Resolution of the face map, x and y axes.

 DEFAULT_GRID_RESOLUTION
density_ratio float

Ratio of real to ideal propellant density.

 1.0
Source code in machwave/models/grain/fmm/base.py 
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
def __init__(
    self,
    length: float,
    outer_diameter: float,
    inhibited_surfaces: grain_base.InhibitedSurfaces | None = None,
    grid_resolution: int = DEFAULT_GRID_RESOLUTION,
    density_ratio: float = 1.0,
) -> None:
    """
    Initialize an FMM grain segment.

    Args:
        length: Segment length [m].
        outer_diameter: Outer diameter [m].
        inhibited_surfaces: Surfaces inhibited from burning.
        grid_resolution: Resolution of the face map, x and y axes.
        density_ratio: Ratio of real to ideal propellant density.
    """
    self.grid_resolution = grid_resolution

    # Cache variables:
    self.coordinate_grids = None
    self.outer_diameter_mask = None
    self.masked_face = None
    self.regression_map = None
    self.web_thickness = None

    # Cache variables that only work for a specific web distance:
    self._solid_mask_web = None
    self._solid_mask = None
    self._solid_indices_web = None
    self._solid_indices = None

    super().__init__(
        length=length,
        outer_diameter=outer_diameter,
        inhibited_surfaces=inhibited_surfaces,
        density_ratio=density_ratio,
    )

cells_to_meters(value)

Convert a cell-distance value to [m].

Parameters:

Name Type Description Default
value float | NDArray[float64]

Distance in cell units.

 required

Returns:

Type Description
float | NDArray[float64]

Distance [m].
Source code in machwave/models/grain/fmm/base.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
def cells_to_meters(
    self, value: float | NDArray[np.float64]
) -> float | NDArray[np.float64]:
    """
    Convert a cell-distance value to [m].

    Args:
        value: Distance in cell units.

    Returns:
        Distance [m].
    """
    return value * self.get_radial_grid_spacing()

cells_to_square_meters(value)

Convert a cell-area value to [m^2].

Parameters:

Name Type Description Default
value float | NDArray[float64]

Area in cell units.

 required

Returns:

Type Description
float | NDArray[float64]

Area [m^2].
Source code in machwave/models/grain/fmm/base.py 
135
136
137
138
139
140
141
142
143
144
145
146
147
def cells_to_square_meters(
    self, value: float | NDArray[np.float64]
) -> float | NDArray[np.float64]:
    """
    Convert a cell-area value to [m^2].

    Args:
        value: Area in cell units.

    Returns:
        Area [m^2].
    """
    return value * self.get_radial_grid_spacing() ** 2

denormalize(value)

Convert a normalized input back into a dimensional value [m].

Source code in machwave/models/grain/fmm/base.py 
106
107
108
109
110
def denormalize(
    self, value: int | float | NDArray[np.float64]
) -> float | NDArray[np.float64]:
    """Convert a normalized input back into a dimensional value [m]."""
    return (value / 2) * (self.outer_diameter)

generate_initial_face_map() abstractmethod

Generate the initial face map for the geometry (1 = propellant, 0 = void).

Source code in machwave/models/grain/fmm/base.py 
85
86
87
88
@abstractmethod
def generate_initial_face_map(self) -> np.typing.NDArray[np.int_]:
    """Generate the initial face map for the geometry (1 = propellant, 0 = void)."""
    pass

get_contours(web_distance, *args, **kwargs) abstractmethod

Return the contours of the regression map for a given web distance.

Must be implemented by a subclass to compute the contour data based on the given web distance and any additional parameters.

Parameters:

Name Type Description Default
web_distance float

Depth of regression into the grain web [m].

 required
*args Any

Additional positional arguments.

 ()
**kwargs Any

Additional keyword arguments.

 {}

Returns:

Type Description
list[NDArray[float64]]

An array representing the computed contours of the grain regression.
Source code in machwave/models/grain/fmm/base.py 
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
@abstractmethod
def get_contours(
    self, web_distance: float, *args: Any, **kwargs: Any
) -> list[NDArray[np.float64]]:
    """
    Return the contours of the regression map for a given web distance.

    Must be implemented by a subclass to compute the contour data based on
    the given web distance and any additional parameters.

    Args:
        web_distance: Depth of regression into the grain web [m].
        *args: Additional positional arguments.
        **kwargs: Additional keyword arguments.

    Returns:
        An array representing the computed contours of the grain regression.
    """
    pass

get_coordinate_grids() abstractmethod

Return the coordinate grids for the grain, one per spatial axis.

Source code in machwave/models/grain/fmm/base.py 
75
76
77
78
@abstractmethod
def get_coordinate_grids(self) -> tuple:
    """Return the coordinate grids for the grain, one per spatial axis."""
    pass

get_empty_face_map()

Return a face map of all ones, shaped like the first stored map.

Source code in machwave/models/grain/fmm/base.py 
90
91
92
def get_empty_face_map(self) -> np.ndarray:
    """Return a face map of all ones, shaped like the first stored map."""
    return np.ones_like(self.get_coordinate_grids()[0])

get_face_map(web_distance)

Returns a matrix representing the grain face based on the given web distance.

The returned array can contain: -1 for masked or invalid points, 0 for points below the threshold, 1 for points above the threshold.

Parameters:

Name Type Description Default
web_distance float

The distance traveled into the grain web.

 required

Returns:

Type Description
NDArray[int64]

A NumPy array with -1, 0, or 1 indicating the grain face at the specified web distance.
Source code in machwave/models/grain/fmm/base.py 
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
def get_face_map(self, web_distance: float) -> np.typing.NDArray[np.int64]:
    """
    Returns a matrix representing the grain face based on the given web distance.

    The returned array can contain:
    -1 for masked or invalid points,
    0 for points below the threshold,
    1 for points above the threshold.

    Args:
        web_distance: The distance traveled into the grain web.

    Returns:
        A NumPy array with -1, 0, or 1 indicating the grain face at the specified web distance.
    """
    web_distance_normalized = self.normalize(web_distance)
    regression_map = self.get_regression_map()
    excluded_mask = np.ma.getmaskarray(regression_map)

    # Create a masked array, where excluded cells are masked out
    occupancy_state = np.ma.MaskedArray(
        (regression_map > web_distance_normalized).astype(np.int64),
        mask=excluded_mask,
    )

    # Fill masked entries with -1, valid/true entries remain 1 or 0
    return occupancy_state.filled(-1)

get_masked_face()

Return a masked representation of the face map.

Source code in machwave/models/grain/fmm/base.py 
176
177
178
179
180
181
182
183
184
185
186
187
def get_masked_face(self) -> np.ndarray:
    """Return a masked representation of the face map."""
    if self.masked_face is None:
        face_map = self.generate_initial_face_map().copy()
        excluded_mask = self.get_outer_diameter_mask().copy()

        face_map, excluded_mask = self._apply_surface_inhibition(
            face_map, excluded_mask
        )

        self.masked_face = np.ma.MaskedArray(face_map, excluded_mask)
    return self.masked_face

get_outer_diameter_mask() abstractmethod

Return a boolean mask of cells outside the outer-diameter boundary.

Source code in machwave/models/grain/fmm/base.py 
80
81
82
83
@abstractmethod
def get_outer_diameter_mask(self) -> np.ndarray:
    """Return a boolean mask of cells outside the outer-diameter boundary."""
    pass

get_radial_grid_spacing()

Return the distance between adjacent cross-section cells [m].

The grid spans the outer diameter over grid_resolution samples, so it holds grid_resolution - 1 intervals.

Source code in machwave/models/grain/fmm/base.py 
112
113
114
115
116
117
118
119
def get_radial_grid_spacing(self) -> float:
    """
    Return the distance between adjacent cross-section cells [m].

    The grid spans the outer diameter over `grid_resolution` samples, so it
    holds `grid_resolution - 1` intervals.
    """
    return self.outer_diameter / (self.grid_resolution - 1)

get_regression_map()

Return the distance map for grain regression.

Uses the fast marching method (scikit-fmm) on the masked face. Each value is the distance from the initial face along the cross-section. With no burning surface (end burner), returns a static map set to the axial burnout web, so the grain regresses along its length, not a cell.

Source code in machwave/models/grain/fmm/base.py 
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
def get_regression_map(self):
    """
    Return the distance map for grain regression.

    Uses the fast marching method (scikit-fmm) on the masked face. Each
    value is the distance from the initial face along the cross-section.
    With no burning surface (end burner), returns a static map set to the
    axial burnout web, so the grain regresses along its length, not a cell.
    """
    if self.regression_map is None:
        masked_face = self.get_masked_face()

        if self.has_cross_section_regression:
            self.regression_map = self._compute_regression_distance(masked_face)
        else:
            # End burner: web = length split across exposed ends. Fill the
            # whole array with one constant (not 0 outside the mask) so the
            # perimeter tracer finds no spurious contour.
            exposed_end_count = (not self.inhibited_surfaces.upper_end) + (
                not self.inhibited_surfaces.lower_end
            )
            normalized_axial_web_distance = (
                self.normalize(self.length / exposed_end_count)
                if exposed_end_count
                else 0.0
            )
            self.regression_map = np.ma.MaskedArray(
                np.full(masked_face.shape, normalized_axial_web_distance),
                mask=np.ma.getmaskarray(masked_face),
            )
    return self.regression_map

get_web_thickness()

Return the maximum web thickness of the grain [m].

Derived from the distance map as the largest normalized distance, converted back into real units.

Source code in machwave/models/grain/fmm/base.py 
233
234
235
236
237
238
239
240
241
242
243
244
def get_web_thickness(self) -> float:
    """
    Return the maximum web thickness of the grain [m].

    Derived from the distance map as the largest normalized distance,
    converted back into real units.
    """
    if self.web_thickness is None:
        self.web_thickness = float(
            self.denormalize(float(np.amax(self.get_regression_map())))
        )
    return float(self.web_thickness)

normalize(value)

Convert a dimensional value to a fraction of the outer radius.

Parameters:

Name Type Description Default
value int | float

Dimensional value (e.g., length) to normalize [m].

 required

Returns:

Type Description
float

Value expressed as a fraction of the outer radius.
Source code in machwave/models/grain/fmm/base.py 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
def normalize(self, value: int | float) -> float:
    """
    Convert a dimensional value to a fraction of the outer radius.

    Args:
        value: Dimensional value (e.g., length) to normalize [m].

    Returns:
        Value expressed as a fraction of the outer radius.
    """
    return value / (0.5 * self.outer_diameter)

validate()

Validate the internal geometry of the grain.

Raises:

Type Description
GrainGeometryError

If the grid resolution is below the valid threshold.
Source code in machwave/models/grain/fmm/base.py 
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
def validate(self) -> None:
    """
    Validate the internal geometry of the grain.

    Raises:
        GrainGeometryError: If the grid resolution is below the valid
            threshold.
    """
    super().validate()

    if not self.grid_resolution >= MINIMUM_GRID_RESOLUTION:
        raise grain.GrainGeometryError(
            f"Grid resolution must be at least {MINIMUM_GRID_RESOLUTION}, "
            f"got {self.grid_resolution}"
        )

FMMGrainSegment2D

Bases: FMMGrainSegment, GrainSegment2D, ABC

Fast Marching Method (FMM) implementation for 2D grain segment.

Source code in machwave/models/grain/fmm/_2d.py 
 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
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
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
class FMMGrainSegment2D(fmm_base.FMMGrainSegment, grain.GrainSegment2D, ABC):
    """Fast Marching Method (FMM) implementation for 2D grain segment."""

    def __init__(
        self,
        length: float,
        outer_diameter: float,
        inhibited_surfaces: grain_base.InhibitedSurfaces | None = None,
        grid_resolution: int = fmm_base.DEFAULT_GRID_RESOLUTION,
        density_ratio: float = 1.0,
    ) -> None:
        self.face_area_interpolator: Callable[[float], float] | None = None
        self.burn_area_interpolator: Callable[[float], float] | None = None
        super().__init__(
            length=length,
            outer_diameter=outer_diameter,
            inhibited_surfaces=inhibited_surfaces,
            grid_resolution=grid_resolution,
            density_ratio=density_ratio,
        )

    def get_coordinate_grids(self) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
        """
        Return the coordinate grids `(map_x, map_y)`.

        Each array has shape `(grid_resolution, grid_resolution)` and ranges from -1 to 1.
        """
        if self.coordinate_grids is None:
            map_x, map_y = np.meshgrid(
                np.linspace(-1, 1, self.grid_resolution, dtype=np.float64),
                np.linspace(-1, 1, self.grid_resolution, dtype=np.float64),
            )
            self.coordinate_grids = (map_x, map_y)
        return self.coordinate_grids

    def get_outer_diameter_mask(self) -> NDArray[np.bool_]:
        """Return a boolean mask indicating which points lie outside the unit circle."""
        if self.outer_diameter_mask is None:
            map_x, map_y = self.get_coordinate_grids()
            self.outer_diameter_mask = (map_x**2 + map_y**2) > 1
        return self.outer_diameter_mask

    def _apply_surface_inhibition(
        self,
        face_map: NDArray[np.int_],
        excluded_mask: NDArray[np.bool_],
    ) -> tuple[NDArray[np.int_], NDArray[np.bool_]]:
        inner_surface_inhibited_cells = (
            (face_map == 0) if self.inhibited_surfaces.inner_surface else None
        )
        face_map, excluded_mask = super()._apply_surface_inhibition(
            face_map, excluded_mask
        )
        if inner_surface_inhibited_cells is not None:
            excluded_mask = excluded_mask | inner_surface_inhibited_cells
        return face_map, excluded_mask

    def get_contours(self, web_distance: float) -> list[NDArray[np.float64]]:
        """
        Return contour arrays for the given web distance.

        Each contour is typically an `(N, 2)` array of `(row, col)` points.
        """
        iso_level = self.normalize(web_distance)
        return fmm_contours.get_iso_contours(self.get_regression_map(), iso_level)

    def get_face_area_interpolator(self) -> Callable[[float], float]:
        """Return an interpolator mapping normalized web distance to face area [m^2]."""
        if self.face_area_interpolator is None:
            regression_map = self.get_regression_map()
            valid = np.logical_not(self.get_outer_diameter_mask())

            # Build face-area curve without per-step full-map scans
            regression_distances_sorted = np.asarray(
                regression_map[valid], dtype=np.float64
            ).ravel()
            regression_distances_sorted.sort()
            max_regression_distance = (
                float(regression_distances_sorted[-1])
                if regression_distances_sorted.size
                else 0.0
            )

            iso_level_count = int(max_regression_distance * self.grid_resolution) + 2
            iso_levels_normalized = (
                np.arange(iso_level_count, dtype=np.float64) / self.grid_resolution
            )

            count_at_or_below_level = np.searchsorted(
                regression_distances_sorted, iso_levels_normalized, side="right"
            )
            solid_cell_count = float(
                regression_distances_sorted.size
            ) - count_at_or_below_level.astype(np.float64)
            face_area_per_iso_level = np.asarray(
                self.cells_to_square_meters(solid_cell_count), dtype=np.float64
            )

            face_area_smoothed = filters.smooth_savitzky_golay(face_area_per_iso_level)
            self.face_area_interpolator = interp1d(
                iso_levels_normalized,
                face_area_smoothed,
                bounds_error=False,
                fill_value=(
                    float(face_area_smoothed[0]),
                    float(face_area_smoothed[-1]),
                ),  # type: ignore[arg-type]
                assume_sorted=True,
            )

        return self.face_area_interpolator

    def get_face_area(self, web_distance: float) -> float:
        """Return the face area at the given web distance."""
        web_distance_normalized = self.normalize(web_distance)
        return float(self.get_face_area_interpolator()(web_distance_normalized))

    def get_port_area(self, web_distance: float) -> float:
        """Return the open cross-sectional port area [m^2]."""
        face_area = self.get_face_area(web_distance)
        return geometric.get_circle_area(self.outer_diameter) - face_area

    def get_burn_area_interpolator(self) -> Callable[[float], float]:
        """Return a cached interpolator for burn area [m^2] vs normalized web."""
        if self.burn_area_interpolator is None:
            regression_map = self.get_regression_map()
            valid = np.logical_not(self.get_outer_diameter_mask())
            regression_distances_sorted = np.asarray(
                regression_map[valid], dtype=np.float64
            ).ravel()
            if regression_distances_sorted.size == 0:
                self.burn_area_interpolator = interp1d(
                    np.asarray([0.0], dtype=np.float64),
                    np.asarray([0.0], dtype=np.float64),
                    bounds_error=False,
                    fill_value=0.0,
                    assume_sorted=True,
                )
                return self.burn_area_interpolator

            regression_distances_sorted.sort()
            max_regression_distance = float(regression_distances_sorted[-1])
            iso_level_count = int(max_regression_distance * self.grid_resolution) + 2
            iso_levels_normalized = (
                np.arange(iso_level_count, dtype=np.float64) / self.grid_resolution
            )

            count_at_or_below_level = np.searchsorted(
                regression_distances_sorted, iso_levels_normalized, side="right"
            )
            solid_cell_count = float(
                regression_distances_sorted.size
            ) - count_at_or_below_level.astype(np.float64)
            face_area_per_iso_level = np.asarray(
                self.cells_to_square_meters(solid_cell_count), dtype=np.float64
            )

            core_perimeter_per_iso_level = np.empty_like(
                iso_levels_normalized, dtype=np.float64
            )
            for i, dist in enumerate(iso_levels_normalized):
                contours = fmm_contours.get_iso_contours(regression_map, float(dist))
                core_perimeter_per_iso_level[i] = float(
                    sum(
                        self.cells_to_meters(fmm_contours.get_length(contour))
                        for contour in contours
                    )
                )

            web_distances_denormalized = np.asarray(
                self.denormalize(iso_levels_normalized), dtype=np.float64
            )
            grain_length_per_web = np.asarray(
                [self.get_length(float(wd)) for wd in web_distances_denormalized],
                dtype=np.float64,
            )
            core_area_per_iso_level = (
                core_perimeter_per_iso_level * grain_length_per_web
            )
            exposed_end_count = (not self.inhibited_surfaces.upper_end) + (
                not self.inhibited_surfaces.lower_end
            )
            exposed_end_area_per_iso_level = exposed_end_count * face_area_per_iso_level
            burn_area_per_iso_level = (
                core_area_per_iso_level + exposed_end_area_per_iso_level
            )

            burn_area_smoothed = filters.smooth_savitzky_golay(burn_area_per_iso_level)
            self.burn_area_interpolator = interp1d(
                iso_levels_normalized,
                burn_area_smoothed,
                bounds_error=False,
                fill_value=(
                    float(burn_area_smoothed[0]),
                    float(burn_area_smoothed[-1]),
                ),  # type: ignore[arg-type]
                assume_sorted=True,
            )

        return self.burn_area_interpolator

    def get_burn_area(self, web_distance: float) -> float:
        """Return burn area [m^2] at a given web distance."""
        if web_distance > self.get_web_thickness():
            return 0.0
        web_distance_normalized = self.normalize(web_distance)
        value = float(self.get_burn_area_interpolator()(web_distance_normalized))
        return max(0.0, value)

    def get_core_perimeter(self, web_distance: float) -> float:
        """Return the perimeter of the open core at the given web distance."""
        contours = self.get_contours(web_distance)
        return float(
            sum(
                self.cells_to_meters(fmm_contours.get_length(contour))
                for contour in contours
            )
        )

    def get_core_area(self, web_distance: float) -> float:
        """
        Return the core (internal) area [m^2] at the given web distance.

        Computed as `perimeter * length` (2D approximation).
        """
        return self.get_core_perimeter(web_distance) * self.get_length(web_distance)

    def _validate_web_distance(self, web_distance: float) -> None:
        """
        Validate that web distance does not exceed web thickness.

        Raises:
            GrainGeometryError: If web distance exceeds web thickness.
        """
        if web_distance > self.get_web_thickness():
            raise grain.GrainGeometryError(
                "The web distance traveled is greater than the grain segment's web thickness."
            )

    def _find_solid_material_indices(
        self, web_distance: float
    ) -> tuple[NDArray[np.int_], NDArray[np.int_]]:
        """
        Get indices of active material at given web distance.

        Args:
            web_distance: Web distance traveled [m].

        Returns:
            Tuple of (y_indices, x_indices) for active material.

        Raises:
            GrainGeometryError: If no active material is found.
        """
        y_indices, x_indices = self._get_solid_indices(web_distance)

        if len(x_indices) == 0 or len(y_indices) == 0:
            raise grain.GrainGeometryError(
                "No active material found at the given web distance."
            )

        return y_indices, x_indices

    def _indices_to_grid_coordinates(
        self, y_indices: NDArray[np.int_], x_indices: NDArray[np.int_]
    ) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
        """
        Convert indices to normalized coordinates centered at origin.

        Args:
            y_indices: Y-axis indices from face map.
            x_indices: X-axis indices from face map.

        Returns:
            Tuple of (x_grid, y_grid) in normalized units.
        """
        grid_center_index = (self.grid_resolution - 1) / 2
        x_grid = (x_indices - grid_center_index).astype(np.float64)
        y_grid = (y_indices - grid_center_index).astype(np.float64)
        return x_grid, y_grid

    def get_center_of_gravity(self, web_distance: float) -> NDArray[np.float64]:
        """
        Return the center of gravity of a 2D FMM grain segment.

        Args:
            web_distance: Web distance traveled [m].

        Returns:
            Center of gravity as `(x, y, z)` [m], measured from the segment port.
        """
        self._validate_web_distance(web_distance)
        y_indices, x_indices = self._find_solid_material_indices(web_distance)
        x_grid, y_grid = self._indices_to_grid_coordinates(y_indices, x_indices)

        # Convert to physical coordinates
        x_denormalized = np.asarray(self.cells_to_meters(x_grid), dtype=np.float64)
        y_denormalized = np.asarray(self.cells_to_meters(y_grid), dtype=np.float64)

        # Axial bounds accounting for end face inhibition
        lower_end_recession = (
            web_distance if not self.inhibited_surfaces.lower_end else 0.0
        )
        upper_end_recession = (
            web_distance if not self.inhibited_surfaces.upper_end else 0.0
        )
        axial_center_of_gravity_position = (
            lower_end_recession + (self.length - upper_end_recession)
        ) / 2.0

        z_denormalized = np.full_like(x_denormalized, axial_center_of_gravity_position)

        return mechanics.get_center_of_gravity(
            x_denormalized, y_denormalized, z_denormalized
        )

    def get_moment_of_inertia(
        self, ideal_density: float, web_distance: float = 0.0
    ) -> NDArray[np.float64]:
        """
        Calculate the moment of inertia tensor at the segment's its center of gravity.

        Args:
            web_distance: Web distance traveled [m].
            ideal_density: Propellant ideal density [kg/m^3].

        Returns:
            A 3x3 inertia tensor [kg-m^2].
        """
        self._validate_web_distance(web_distance)
        y_indices, x_indices = self._find_solid_material_indices(web_distance)
        x_grid, y_grid = self._indices_to_grid_coordinates(y_indices, x_indices)

        current_length = self.get_length(web_distance)

        # Convert to meters
        x_denormalized = np.asarray(self.cells_to_meters(x_grid), dtype=np.float64)
        y_denormalized = np.asarray(self.cells_to_meters(y_grid), dtype=np.float64)

        # Radial CoG from the same coordinates, avoiding an index re-extraction.
        # Axial CoG is unused here: 2D elements have no axial spread.
        center_of_gravity = mechanics.get_center_of_gravity(
            x_denormalized, y_denormalized, np.zeros_like(x_denormalized)
        )

        # Coordinates relative to the center of gravity
        x_relative = x_denormalized - center_of_gravity[1]
        y_relative = y_denormalized - center_of_gravity[2]
        z_relative = np.zeros_like(x_relative)  # 2D grain, no z variation in elements

        # Total mass
        total_volume = self.get_volume(web_distance)
        total_mass = total_volume * ideal_density * self.density_ratio

        # Mass per cross-sectional element
        solid_element_count = len(x_indices)
        element_mass = total_mass / solid_element_count

        # Get base inertia from point masses (radial only)
        inertia_tensor = mechanics.get_moment_of_inertia_tensor(
            x_relative, y_relative, z_relative, element_mass
        )

        # For 2D grain (uniform along axial direction), add axial contribution
        # For a uniform rod of length L with mass M: axial_inertia_contribution = M*L^2/12 (about the center of gravity)
        axial_inertia_contribution = total_mass * current_length**2 / 12

        # Add axial contribution to radial axes (indices 1 and 2 in [z,x,y] system)
        inertia_tensor[1, 1] += axial_inertia_contribution  # Ixx (moment about x-axis)
        inertia_tensor[2, 2] += axial_inertia_contribution  # Iyy (moment about y-axis)

        return inertia_tensor

get_burn_area(web_distance)

Return burn area [m^2] at a given web distance.

Source code in machwave/models/grain/fmm/_2d.py 
219
220
221
222
223
224
225
def get_burn_area(self, web_distance: float) -> float:
    """Return burn area [m^2] at a given web distance."""
    if web_distance > self.get_web_thickness():
        return 0.0
    web_distance_normalized = self.normalize(web_distance)
    value = float(self.get_burn_area_interpolator()(web_distance_normalized))
    return max(0.0, value)

get_burn_area_interpolator()

Return a cached interpolator for burn area [m^2] vs normalized web.

Source code in machwave/models/grain/fmm/_2d.py 
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
def get_burn_area_interpolator(self) -> Callable[[float], float]:
    """Return a cached interpolator for burn area [m^2] vs normalized web."""
    if self.burn_area_interpolator is None:
        regression_map = self.get_regression_map()
        valid = np.logical_not(self.get_outer_diameter_mask())
        regression_distances_sorted = np.asarray(
            regression_map[valid], dtype=np.float64
        ).ravel()
        if regression_distances_sorted.size == 0:
            self.burn_area_interpolator = interp1d(
                np.asarray([0.0], dtype=np.float64),
                np.asarray([0.0], dtype=np.float64),
                bounds_error=False,
                fill_value=0.0,
                assume_sorted=True,
            )
            return self.burn_area_interpolator

        regression_distances_sorted.sort()
        max_regression_distance = float(regression_distances_sorted[-1])
        iso_level_count = int(max_regression_distance * self.grid_resolution) + 2
        iso_levels_normalized = (
            np.arange(iso_level_count, dtype=np.float64) / self.grid_resolution
        )

        count_at_or_below_level = np.searchsorted(
            regression_distances_sorted, iso_levels_normalized, side="right"
        )
        solid_cell_count = float(
            regression_distances_sorted.size
        ) - count_at_or_below_level.astype(np.float64)
        face_area_per_iso_level = np.asarray(
            self.cells_to_square_meters(solid_cell_count), dtype=np.float64
        )

        core_perimeter_per_iso_level = np.empty_like(
            iso_levels_normalized, dtype=np.float64
        )
        for i, dist in enumerate(iso_levels_normalized):
            contours = fmm_contours.get_iso_contours(regression_map, float(dist))
            core_perimeter_per_iso_level[i] = float(
                sum(
                    self.cells_to_meters(fmm_contours.get_length(contour))
                    for contour in contours
                )
            )

        web_distances_denormalized = np.asarray(
            self.denormalize(iso_levels_normalized), dtype=np.float64
        )
        grain_length_per_web = np.asarray(
            [self.get_length(float(wd)) for wd in web_distances_denormalized],
            dtype=np.float64,
        )
        core_area_per_iso_level = (
            core_perimeter_per_iso_level * grain_length_per_web
        )
        exposed_end_count = (not self.inhibited_surfaces.upper_end) + (
            not self.inhibited_surfaces.lower_end
        )
        exposed_end_area_per_iso_level = exposed_end_count * face_area_per_iso_level
        burn_area_per_iso_level = (
            core_area_per_iso_level + exposed_end_area_per_iso_level
        )

        burn_area_smoothed = filters.smooth_savitzky_golay(burn_area_per_iso_level)
        self.burn_area_interpolator = interp1d(
            iso_levels_normalized,
            burn_area_smoothed,
            bounds_error=False,
            fill_value=(
                float(burn_area_smoothed[0]),
                float(burn_area_smoothed[-1]),
            ),  # type: ignore[arg-type]
            assume_sorted=True,
        )

    return self.burn_area_interpolator

get_center_of_gravity(web_distance)

Return the center of gravity of a 2D FMM grain segment.

Parameters:

Name Type Description Default
web_distance float

Web distance traveled [m].

 required

Returns:

Type Description
NDArray[float64]

Center of gravity as (x, y, z) [m], measured from the segment port.
Source code in machwave/models/grain/fmm/_2d.py 
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
def get_center_of_gravity(self, web_distance: float) -> NDArray[np.float64]:
    """
    Return the center of gravity of a 2D FMM grain segment.

    Args:
        web_distance: Web distance traveled [m].

    Returns:
        Center of gravity as `(x, y, z)` [m], measured from the segment port.
    """
    self._validate_web_distance(web_distance)
    y_indices, x_indices = self._find_solid_material_indices(web_distance)
    x_grid, y_grid = self._indices_to_grid_coordinates(y_indices, x_indices)

    # Convert to physical coordinates
    x_denormalized = np.asarray(self.cells_to_meters(x_grid), dtype=np.float64)
    y_denormalized = np.asarray(self.cells_to_meters(y_grid), dtype=np.float64)

    # Axial bounds accounting for end face inhibition
    lower_end_recession = (
        web_distance if not self.inhibited_surfaces.lower_end else 0.0
    )
    upper_end_recession = (
        web_distance if not self.inhibited_surfaces.upper_end else 0.0
    )
    axial_center_of_gravity_position = (
        lower_end_recession + (self.length - upper_end_recession)
    ) / 2.0

    z_denormalized = np.full_like(x_denormalized, axial_center_of_gravity_position)

    return mechanics.get_center_of_gravity(
        x_denormalized, y_denormalized, z_denormalized
    )

get_contours(web_distance)

Return contour arrays for the given web distance.

Each contour is typically an (N, 2) array of (row, col) points.

Source code in machwave/models/grain/fmm/_2d.py 
75
76
77
78
79
80
81
82
def get_contours(self, web_distance: float) -> list[NDArray[np.float64]]:
    """
    Return contour arrays for the given web distance.

    Each contour is typically an `(N, 2)` array of `(row, col)` points.
    """
    iso_level = self.normalize(web_distance)
    return fmm_contours.get_iso_contours(self.get_regression_map(), iso_level)

get_coordinate_grids()

Return the coordinate grids (map_x, map_y).

Each array has shape (grid_resolution, grid_resolution) and ranges from -1 to 1.

Source code in machwave/models/grain/fmm/_2d.py 
39
40
41
42
43
44
45
46
47
48
49
50
51
def get_coordinate_grids(self) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
    """
    Return the coordinate grids `(map_x, map_y)`.

    Each array has shape `(grid_resolution, grid_resolution)` and ranges from -1 to 1.
    """
    if self.coordinate_grids is None:
        map_x, map_y = np.meshgrid(
            np.linspace(-1, 1, self.grid_resolution, dtype=np.float64),
            np.linspace(-1, 1, self.grid_resolution, dtype=np.float64),
        )
        self.coordinate_grids = (map_x, map_y)
    return self.coordinate_grids

get_core_area(web_distance)

Return the core (internal) area [m^2] at the given web distance.

Computed as perimeter * length (2D approximation).

Source code in machwave/models/grain/fmm/_2d.py 
237
238
239
240
241
242
243
def get_core_area(self, web_distance: float) -> float:
    """
    Return the core (internal) area [m^2] at the given web distance.

    Computed as `perimeter * length` (2D approximation).
    """
    return self.get_core_perimeter(web_distance) * self.get_length(web_distance)

get_core_perimeter(web_distance)

Return the perimeter of the open core at the given web distance.

Source code in machwave/models/grain/fmm/_2d.py 
227
228
229
230
231
232
233
234
235
def get_core_perimeter(self, web_distance: float) -> float:
    """Return the perimeter of the open core at the given web distance."""
    contours = self.get_contours(web_distance)
    return float(
        sum(
            self.cells_to_meters(fmm_contours.get_length(contour))
            for contour in contours
        )
    )

get_face_area(web_distance)

Return the face area at the given web distance.

Source code in machwave/models/grain/fmm/_2d.py 
130
131
132
133
def get_face_area(self, web_distance: float) -> float:
    """Return the face area at the given web distance."""
    web_distance_normalized = self.normalize(web_distance)
    return float(self.get_face_area_interpolator()(web_distance_normalized))

get_face_area_interpolator()

Return an interpolator mapping normalized web distance to face area [m^2].

Source code in machwave/models/grain/fmm/_2d.py 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
def get_face_area_interpolator(self) -> Callable[[float], float]:
    """Return an interpolator mapping normalized web distance to face area [m^2]."""
    if self.face_area_interpolator is None:
        regression_map = self.get_regression_map()
        valid = np.logical_not(self.get_outer_diameter_mask())

        # Build face-area curve without per-step full-map scans
        regression_distances_sorted = np.asarray(
            regression_map[valid], dtype=np.float64
        ).ravel()
        regression_distances_sorted.sort()
        max_regression_distance = (
            float(regression_distances_sorted[-1])
            if regression_distances_sorted.size
            else 0.0
        )

        iso_level_count = int(max_regression_distance * self.grid_resolution) + 2
        iso_levels_normalized = (
            np.arange(iso_level_count, dtype=np.float64) / self.grid_resolution
        )

        count_at_or_below_level = np.searchsorted(
            regression_distances_sorted, iso_levels_normalized, side="right"
        )
        solid_cell_count = float(
            regression_distances_sorted.size
        ) - count_at_or_below_level.astype(np.float64)
        face_area_per_iso_level = np.asarray(
            self.cells_to_square_meters(solid_cell_count), dtype=np.float64
        )

        face_area_smoothed = filters.smooth_savitzky_golay(face_area_per_iso_level)
        self.face_area_interpolator = interp1d(
            iso_levels_normalized,
            face_area_smoothed,
            bounds_error=False,
            fill_value=(
                float(face_area_smoothed[0]),
                float(face_area_smoothed[-1]),
            ),  # type: ignore[arg-type]
            assume_sorted=True,
        )

    return self.face_area_interpolator

get_moment_of_inertia(ideal_density, web_distance=0.0)

Calculate the moment of inertia tensor at the segment's its center of gravity.

Parameters:

Name Type Description Default
web_distance float

Web distance traveled [m].

 0.0
ideal_density float

Propellant ideal density [kg/m^3].

 required

Returns:

Type Description
NDArray[float64]

A 3x3 inertia tensor [kg-m^2].
Source code in machwave/models/grain/fmm/_2d.py 
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
def get_moment_of_inertia(
    self, ideal_density: float, web_distance: float = 0.0
) -> NDArray[np.float64]:
    """
    Calculate the moment of inertia tensor at the segment's its center of gravity.

    Args:
        web_distance: Web distance traveled [m].
        ideal_density: Propellant ideal density [kg/m^3].

    Returns:
        A 3x3 inertia tensor [kg-m^2].
    """
    self._validate_web_distance(web_distance)
    y_indices, x_indices = self._find_solid_material_indices(web_distance)
    x_grid, y_grid = self._indices_to_grid_coordinates(y_indices, x_indices)

    current_length = self.get_length(web_distance)

    # Convert to meters
    x_denormalized = np.asarray(self.cells_to_meters(x_grid), dtype=np.float64)
    y_denormalized = np.asarray(self.cells_to_meters(y_grid), dtype=np.float64)

    # Radial CoG from the same coordinates, avoiding an index re-extraction.
    # Axial CoG is unused here: 2D elements have no axial spread.
    center_of_gravity = mechanics.get_center_of_gravity(
        x_denormalized, y_denormalized, np.zeros_like(x_denormalized)
    )

    # Coordinates relative to the center of gravity
    x_relative = x_denormalized - center_of_gravity[1]
    y_relative = y_denormalized - center_of_gravity[2]
    z_relative = np.zeros_like(x_relative)  # 2D grain, no z variation in elements

    # Total mass
    total_volume = self.get_volume(web_distance)
    total_mass = total_volume * ideal_density * self.density_ratio

    # Mass per cross-sectional element
    solid_element_count = len(x_indices)
    element_mass = total_mass / solid_element_count

    # Get base inertia from point masses (radial only)
    inertia_tensor = mechanics.get_moment_of_inertia_tensor(
        x_relative, y_relative, z_relative, element_mass
    )

    # For 2D grain (uniform along axial direction), add axial contribution
    # For a uniform rod of length L with mass M: axial_inertia_contribution = M*L^2/12 (about the center of gravity)
    axial_inertia_contribution = total_mass * current_length**2 / 12

    # Add axial contribution to radial axes (indices 1 and 2 in [z,x,y] system)
    inertia_tensor[1, 1] += axial_inertia_contribution  # Ixx (moment about x-axis)
    inertia_tensor[2, 2] += axial_inertia_contribution  # Iyy (moment about y-axis)

    return inertia_tensor

get_outer_diameter_mask()

Return a boolean mask indicating which points lie outside the unit circle.

Source code in machwave/models/grain/fmm/_2d.py 
53
54
55
56
57
58
def get_outer_diameter_mask(self) -> NDArray[np.bool_]:
    """Return a boolean mask indicating which points lie outside the unit circle."""
    if self.outer_diameter_mask is None:
        map_x, map_y = self.get_coordinate_grids()
        self.outer_diameter_mask = (map_x**2 + map_y**2) > 1
    return self.outer_diameter_mask

get_port_area(web_distance)

Return the open cross-sectional port area [m^2].

Source code in machwave/models/grain/fmm/_2d.py 
135
136
137
138
def get_port_area(self, web_distance: float) -> float:
    """Return the open cross-sectional port area [m^2]."""
    face_area = self.get_face_area(web_distance)
    return geometric.get_circle_area(self.outer_diameter) - face_area

FMMGrainSegment3D

Bases: FMMGrainSegment, GrainSegment3D, ABC

Fast Marching Method (FMM) implementation for 3D grain segment.

Source code in machwave/models/grain/fmm/_3d.py 
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
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
class FMMGrainSegment3D(fmm_base.FMMGrainSegment, grain.GrainSegment3D, ABC):
    """Fast Marching Method (FMM) implementation for 3D grain segment."""

    def __init__(
        self,
        length: float,
        outer_diameter: float,
        inhibited_surfaces: grain_base.InhibitedSurfaces | None = None,
        grid_resolution: int = fmm_base.DEFAULT_GRID_RESOLUTION,
        density_ratio: float = 1.0,
    ) -> None:
        self.burn_area_interpolator: Callable[[float], float] | None = None
        self.volume_interpolator: Callable[[float], float] | None = None
        self.padded_regression_map: np.ndarray | None = None

        # Cache center of gravity and moment of inertia shared moments per web distance
        self._mask_moments_web: float | None = None
        self._mask_moments: tuple[NDArray[np.float64], NDArray[np.float64]] | None = (
            None
        )

        super().__init__(
            length=length,
            outer_diameter=outer_diameter,
            inhibited_surfaces=inhibited_surfaces,
            grid_resolution=grid_resolution,
            density_ratio=density_ratio,
        )

    def validate(self) -> None:
        super().validate()
        self._validate_axial_resolution()

    def _validate_axial_resolution(self) -> None:
        # < 3 slices: get_port_area indexes a size-0 axis and inhibition is skipped
        axial_resolution = self.get_axial_resolution()
        if axial_resolution < 3:
            raise grain.GrainGeometryError(
                f"Axial resolution must be at least 3, got "
                f"{axial_resolution}; increase grid_resolution or the "
                f"length-to-outer-diameter ratio."
            )

    def get_axial_resolution(self) -> int:
        return round(self.grid_resolution * self.length / self.outer_diameter)

    def get_axial_grid_spacing(self) -> float:
        """
        Return the thickness of one axial slice [m].

        The slices tile the length, so each one carries a full share of it.
        """
        return self.length / self.get_axial_resolution()

    def get_coordinate_grids(
        self,
    ) -> tuple[
        NDArray[np.float64],
        NDArray[np.float64],
        NDArray[np.float64],
    ]:
        if self.coordinate_grids is None:
            # Slice centers, half a slice in from each end face.
            axial_resolution = self.get_axial_resolution()
            half_slice = 0.5 / axial_resolution
            map_y, map_z, map_x = np.meshgrid(
                np.linspace(-1, 1, self.grid_resolution),
                np.linspace(1 - half_slice, half_slice, axial_resolution),  # z axis
                np.linspace(-1, 1, self.grid_resolution),
            )

            self.coordinate_grids = (map_x, map_y, map_z)

        return self.coordinate_grids

    def get_outer_diameter_mask(self) -> NDArray[np.bool_]:
        if self.outer_diameter_mask is None:
            map_x, map_y, _ = self.get_coordinate_grids()
            self.outer_diameter_mask = (map_x**2 + map_y**2) > 1

        return self.outer_diameter_mask

    def _apply_surface_inhibition(
        self,
        face_map: NDArray[np.int_],
        excluded_mask: NDArray[np.bool_],
    ) -> tuple[NDArray[np.int_], NDArray[np.bool_]]:
        inner_surface_inhibited_cells = face_map == 0

        if not self.inhibited_surfaces.outer_surface:
            # Erode with the end faces padded as inside, so the boundary the
            # erosion finds is the casing wall and not the end slices.
            inside = ~excluded_mask
            padded_inside = np.pad(
                inside, ((1, 1), (0, 0), (0, 0)), constant_values=True
            )
            eroded = binary_erosion(padded_inside)[1:-1]
            face_map[inside & np.logical_not(eroded)] = 0

        if self.inhibited_surfaces.inner_surface:
            excluded_mask = excluded_mask | inner_surface_inhibited_cells

        return face_map, excluded_mask

    @property
    def has_cross_section_regression(self) -> bool:
        """An exposed end face burns even when the cross section holds no core."""
        return super().has_cross_section_regression or self.exposed_end_count > 0

    def _get_exposed_end_padding(self) -> tuple[int, int]:
        """Slices to add beyond the aft and forward end faces, one if exposed."""
        return (
            int(not self.inhibited_surfaces.lower_end),
            int(not self.inhibited_surfaces.upper_end),
        )

    def _pad_exposed_ends(self, masked_face: np.ndarray) -> np.ndarray:
        """
        Place a void slice beyond each exposed end face.

        The burning surface is the boundary between propellant and void, so an
        end face only burns if there is void past it. Padding puts that void
        outside the segment, which keeps the end slices as propellant.
        """
        lower, upper = self._get_exposed_end_padding()
        if not (lower or upper):
            return masked_face

        pad_width = ((lower, upper), (0, 0), (0, 0))
        data = np.pad(np.ma.getdata(masked_face), pad_width, constant_values=0)
        # The casing runs past the end faces, so the pad carries its mask.
        mask = np.pad(np.ma.getmaskarray(masked_face), pad_width, mode="edge")
        return np.ma.MaskedArray(data, mask)

    def _strip_exposed_ends(self, volume: np.ndarray) -> np.ndarray:
        """Drop the padding slices, leaving the segment itself."""
        lower, upper = self._get_exposed_end_padding()
        if not (lower or upper):
            return volume
        return volume[lower : volume.shape[0] - upper]

    def _compute_regression_distance(self, masked_face: np.ndarray) -> np.ndarray:
        # Regression speed needs to be calibrated for the z axes separately from the x
        # and y axes, because the 3D grid is anisotropic
        axial_grid_spacing = self.get_axial_grid_spacing()
        radial_grid_spacing = self.get_radial_grid_spacing()
        distance = skfmm.distance(
            self._pad_exposed_ends(masked_face),
            dx=[axial_grid_spacing, radial_grid_spacing, radial_grid_spacing],  # type: ignore[arg-type]
        )
        # The padded field keeps the end faces closed for the burn area mesh.
        self.padded_regression_map = distance * (2.0 / self.outer_diameter)
        return self._strip_exposed_ends(self.padded_regression_map)

    def get_padded_regression_map(self) -> np.ndarray:
        """Return the regression map including the void slice past each exposed end."""
        regression_map = self.get_regression_map()
        if self.padded_regression_map is None:
            return regression_map
        return self.padded_regression_map

    def get_axial_index(self, axial_position_normalized: float) -> int:
        """
        Return the slice index for an axial position, clamped to the grid.

        Args:
            axial_position_normalized: Axial position as a fraction of the
                segment length, measured from the aft (nozzle) end.

        Returns:
            Index of the slice holding that position, within
            `[0, axial_resolution - 1]`.
        """
        axial_resolution = self.get_axial_resolution()
        axial_index = int(np.floor(axial_position_normalized * axial_resolution))
        return min(max(axial_index, 0), axial_resolution - 1)

    def get_contours(
        self, web_distance: float, axial_position_normalized: float
    ) -> list[NDArray[np.float64]]:
        """
        Return the contours of one axial slice at a given web distance.

        Args:
            web_distance: Distance traveled into the grain web [m].
            axial_position_normalized: Axial position as a fraction of the
                segment length, measured from the aft (nozzle) end.

        Returns:
            Contour arrays of the slice, each an `(N, 2)` array of points.
        """
        iso_level = self.normalize(web_distance)
        axial_slice = self.get_regression_map()[
            self.get_axial_index(axial_position_normalized)
        ]
        return fmm_contours.get_iso_contours(axial_slice, iso_level)

    def get_burn_area_interpolator(self) -> Callable[[float], float]:
        """Return a cached interpolator for burn area [m^2] vs web distance [m]."""
        if self.burn_area_interpolator is None:
            regression_map = self.get_padded_regression_map()
            regression_values = np.asarray(
                regression_map[~np.ma.getmaskarray(regression_map)], dtype=np.float64
            )
            max_iso_level = (
                float(regression_values.max()) if regression_values.size else 0.0
            )
            if max_iso_level <= 0.0:
                self.burn_area_interpolator = interp1d(
                    np.asarray([0.0], dtype=np.float64),
                    np.asarray([0.0], dtype=np.float64),
                    bounds_error=False,
                    fill_value=0.0,
                    assume_sorted=True,
                )
                return self.burn_area_interpolator

            # Lift the inhibited casing above every iso level so marching cubes
            # meshes only the burning front, not the wall.
            distance_field = np.ma.filled(regression_map, max_iso_level + 1.0)
            axial_grid_spacing = self.get_axial_grid_spacing()
            radial_grid_spacing = self.get_radial_grid_spacing()
            grid_spacing = (
                axial_grid_spacing,
                radial_grid_spacing,
                radial_grid_spacing,
            )

            # The iso-surface at level 0 lies on the voxelized initial face and is
            # degenerate, so sample above it and hold the first area back to web 0.
            iso_levels = np.linspace(
                max_iso_level / ISO_LEVEL_COUNT, max_iso_level, ISO_LEVEL_COUNT
            )
            iso_surface_areas = _compute_iso_surface_areas(
                distance_field, iso_levels, grid_spacing
            )
            web_distances_denormalized = np.concatenate(
                ([0.0], np.asarray(self.denormalize(iso_levels), dtype=np.float64))
            )
            iso_surface_areas = np.concatenate(
                ([iso_surface_areas[0]], iso_surface_areas)
            )

            burn_area_smoothed = filters.smooth_savitzky_golay(
                iso_surface_areas, window_length=9, polyorder=3
            )
            self.burn_area_interpolator = interp1d(
                web_distances_denormalized,
                burn_area_smoothed,
                bounds_error=False,
                fill_value=(
                    float(burn_area_smoothed[0]),
                    float(burn_area_smoothed[-1]),
                ),  # type: ignore[arg-type]
                assume_sorted=True,
            )

        return self.burn_area_interpolator

    def get_burn_area(self, web_distance: float) -> float:
        if web_distance > self.get_web_thickness():
            return 0.0
        return max(0.0, float(self.get_burn_area_interpolator()(web_distance)))

    def get_port_area(self, web_distance: float, z: float = 0.0) -> float:
        """
        Calculates the port area at a given web distance and axial height z.

        This method extracts a single 2D slice from the 3D face map by converting
        the physical height z into an integer index, and then computes the port
        area for that slice.

        Args:
            web_distance: The distance traveled into the grain web.
            z: Axial position along the grain [m], measured from the aft
                (nozzle) end. Defaults to nozzle end.

        Returns:
            A float representing the port area at the specified z slice, in m^2.
        """
        iso_level = self.normalize(web_distance)
        valid = np.logical_not(self.get_outer_diameter_mask())
        solid = np.logical_and(self.get_regression_map() > iso_level, valid)

        axial_index = self.get_axial_index(z / self.length)

        face_area = float(
            self.cells_to_square_meters(float(np.count_nonzero(solid[axial_index])))
        )
        return geometric.get_circle_area(self.outer_diameter) - face_area

    def get_minimum_port_area(self, web_distance: float) -> float:
        """
        Return the smallest port area along the segment [m^2].

        The cross section of a 3D segment varies along its length, so the flow
        passes a bottleneck: the station holding the most solid propellant.

        Args:
            web_distance: The distance traveled into the grain web.

        Returns:
            Smallest port area along the segment [m^2].
        """
        iso_level = self.normalize(web_distance)
        valid = np.logical_not(self.get_outer_diameter_mask())
        solid = np.logical_and(self.get_regression_map() > iso_level, valid)

        solid_cells_per_slice = np.count_nonzero(solid, axis=(1, 2))
        face_area = float(
            self.cells_to_square_meters(float(solid_cells_per_slice.max()))
        )
        return geometric.get_circle_area(self.outer_diameter) - face_area

    def get_voxel_volume(self) -> float:
        """Volume of one cell of the anisotropic grid [m^3]."""
        return self.get_radial_grid_spacing() ** 2 * self.get_axial_grid_spacing()

    def get_volume_interpolator(self) -> Callable[[float], float]:
        """Return a cached interpolator for volume [m^3] for a web distance."""
        if self.volume_interpolator is None:
            regression_map = self.get_regression_map()
            regression_distances_sorted = np.asarray(
                regression_map[~np.ma.getmaskarray(regression_map)], dtype=np.float64
            ).ravel()
            regression_distances_sorted.sort()
            max_regression_distance = (
                float(regression_distances_sorted[-1])
                if regression_distances_sorted.size
                else 0.0
            )
            if max_regression_distance <= 0.0:
                self.volume_interpolator = interp1d(
                    np.asarray([0.0], dtype=np.float64),
                    np.asarray([0.0], dtype=np.float64),
                    bounds_error=False,
                    fill_value=0.0,
                    assume_sorted=True,
                )
                return self.volume_interpolator

            iso_level_count = int(max_regression_distance * self.grid_resolution) + 2
            iso_levels_normalized = (
                np.arange(iso_level_count, dtype=np.float64) / self.grid_resolution
            )

            count_at_or_below_level = np.searchsorted(
                regression_distances_sorted, iso_levels_normalized, side="right"
            )
            solid_voxel_count = float(
                regression_distances_sorted.size
            ) - count_at_or_below_level.astype(np.float64)
            volume_per_iso_level = solid_voxel_count * self.get_voxel_volume()

            self.volume_interpolator = interp1d(
                iso_levels_normalized,
                volume_per_iso_level,
                bounds_error=False,
                fill_value=(
                    float(volume_per_iso_level[0]),
                    float(volume_per_iso_level[-1]),
                ),  # type: ignore[arg-type]
                assume_sorted=True,
            )

        return self.volume_interpolator

    def get_volume(self, web_distance: float) -> float:
        web_distance_normalized = self.normalize(web_distance)
        return max(0.0, float(self.get_volume_interpolator()(web_distance_normalized)))

    def _validate_web_distance(self, web_distance: float) -> None:
        """
        Validate that web distance does not exceed web thickness.

        Raises:
            GrainGeometryError: If web distance exceeds web thickness.
        """
        if web_distance > self.get_web_thickness():
            raise grain.GrainGeometryError(
                "The web distance traveled is greater than the grain "
                "segment's web thickness."
            )

    def _solid_mask_moments(
        self, web_distance: float
    ) -> tuple[NDArray[np.float64], NDArray[np.float64]] | None:
        """Centroid and second moments of the solid voxels, memoized per web."""
        if self._mask_moments_web != web_distance:
            self._mask_moments = self._compute_solid_mask_moments(web_distance)
            self._mask_moments_web = web_distance
        return self._mask_moments

    def _compute_solid_mask_moments(
        self, web_distance: float
    ) -> tuple[NDArray[np.float64], NDArray[np.float64]] | None:
        """
        Centroid and centered second moments of the solid voxels.

        Reduces the cached boolean solid mask onto each coordinate plane and
        takes the moments as dot products with the per-axis coordinate vectors,
        so the per-voxel index arrays and conversions are never built.
        Coordinates match the index path: x and y are centered on the grid, z is
        measured from the port end.

        Returns:
            Tuple of the centroid [z, x, y] [m] and the 3x3 symmetric matrix of
            summed (r_i - cog_i)(r_j - cog_j) over the voxels, axes ordered
            [x, y, z] [m^2]. None when no solid voxels remain.
        """
        mask = self._get_solid_mask(web_distance)
        count = int(np.count_nonzero(mask))
        if count == 0:
            return None

        axial_count, y_count, x_count = mask.shape
        grid_center = (self.grid_resolution - 1) / 2
        radial_grid_spacing = self.get_radial_grid_spacing()
        x = (np.arange(x_count, dtype=np.float64) - grid_center) * radial_grid_spacing
        y = (np.arange(y_count, dtype=np.float64) - grid_center) * radial_grid_spacing
        z = (
            np.arange(axial_count, dtype=np.float64) + 0.5
        ) * self.get_axial_grid_spacing()

        # Voxel counts projected onto each coordinate plane.
        projection_yx = mask.sum(axis=0)  # over z -> (y, x)
        projection_zx = mask.sum(axis=1)  # over y -> (z, x)
        projection_zy = mask.sum(axis=2)  # over x -> (z, y)

        count_x = projection_yx.sum(axis=0)
        count_y = projection_yx.sum(axis=1)
        count_z = projection_zx.sum(axis=1)

        # First moments: sum of each physical coordinate over the voxels.
        sum_x = count_x @ x
        sum_y = count_y @ y
        sum_z = count_z @ z
        centroid = np.array([sum_z, sum_x, sum_y], dtype=np.float64) / count

        # Second moments about the origin; the cross terms need the planar
        # projections because their two axes are coupled.
        sum_xx = count_x @ (x * x)
        sum_yy = count_y @ (y * y)
        sum_zz = count_z @ (z * z)
        sum_xy = y @ projection_yx @ x
        sum_xz = z @ projection_zx @ x
        sum_yz = z @ projection_zy @ y

        # Shift to the centroid (parallel-axis theorem on the summed moments).
        central = np.array(
            [
                [sum_xx, sum_xy, sum_xz],
                [sum_xy, sum_yy, sum_yz],
                [sum_xz, sum_yz, sum_zz],
            ],
            dtype=np.float64,
        )
        first = np.array([sum_x, sum_y, sum_z], dtype=np.float64)
        central -= np.outer(first, first) / count
        return centroid, central

    def get_center_of_gravity(self, web_distance: float) -> NDArray[np.float64]:
        """
        Return the center of gravity of a 3D FMM grain segment.

        Args:
            web_distance: Web distance traveled [m].

        Returns:
            Center of gravity [z, x, y] in meters from the port of the segment.

        Raises:
            GrainGeometryError: If web distance exceeds the web thickness, or
                if no active material is found.
        """
        self._validate_web_distance(web_distance)
        moments = self._solid_mask_moments(web_distance)
        if moments is None:
            raise grain.GrainGeometryError(
                "No active material found at the given web distance."
            )
        centroid, _ = moments
        # Copy so callers cannot mutate the per-web cached array.
        return centroid.copy()

    def get_moment_of_inertia(
        self, ideal_density: float, web_distance: float = 0.0
    ) -> NDArray[np.float64]:
        """
        Calculate the moment of inertia tensor at the segment's its center of gravity.

        Args:
            ideal_density: Propellant ideal density [kg/m^3].
            web_distance: Web distance traveled [m].

        Returns:
            A 3x3 inertia tensor [kg-m^2].

        Raises:
            GrainGeometryError: If web distance exceeds web thickness or
                if no active material is found.
        """
        self._validate_web_distance(web_distance)
        moments = self._solid_mask_moments(web_distance)
        if moments is None:
            raise grain.GrainGeometryError(
                "No active material found at the given web distance."
            )
        _, central_second_moments = moments

        element_mass = self.get_voxel_volume() * ideal_density * self.density_ratio
        return mechanics.get_moment_of_inertia_tensor_from_central_moments(
            central_second_moments, element_mass
        )

has_cross_section_regression property

An exposed end face burns even when the cross section holds no core.

get_axial_grid_spacing()

Return the thickness of one axial slice [m].

The slices tile the length, so each one carries a full share of it.

Source code in machwave/models/grain/fmm/_3d.py 
201
202
203
204
205
206
207
def get_axial_grid_spacing(self) -> float:
    """
    Return the thickness of one axial slice [m].

    The slices tile the length, so each one carries a full share of it.
    """
    return self.length / self.get_axial_resolution()

get_axial_index(axial_position_normalized)

Return the slice index for an axial position, clamped to the grid.

Parameters:

Name Type Description Default
axial_position_normalized float

Axial position as a fraction of the segment length, measured from the aft (nozzle) end.

 required

Returns:

Type Description
int

Index of the slice holding that position, within
int

[0, axial_resolution - 1].
Source code in machwave/models/grain/fmm/_3d.py 
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
def get_axial_index(self, axial_position_normalized: float) -> int:
    """
    Return the slice index for an axial position, clamped to the grid.

    Args:
        axial_position_normalized: Axial position as a fraction of the
            segment length, measured from the aft (nozzle) end.

    Returns:
        Index of the slice holding that position, within
        `[0, axial_resolution - 1]`.
    """
    axial_resolution = self.get_axial_resolution()
    axial_index = int(np.floor(axial_position_normalized * axial_resolution))
    return min(max(axial_index, 0), axial_resolution - 1)

get_burn_area_interpolator()

Return a cached interpolator for burn area [m^2] vs web distance [m].

Source code in machwave/models/grain/fmm/_3d.py 
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
def get_burn_area_interpolator(self) -> Callable[[float], float]:
    """Return a cached interpolator for burn area [m^2] vs web distance [m]."""
    if self.burn_area_interpolator is None:
        regression_map = self.get_padded_regression_map()
        regression_values = np.asarray(
            regression_map[~np.ma.getmaskarray(regression_map)], dtype=np.float64
        )
        max_iso_level = (
            float(regression_values.max()) if regression_values.size else 0.0
        )
        if max_iso_level <= 0.0:
            self.burn_area_interpolator = interp1d(
                np.asarray([0.0], dtype=np.float64),
                np.asarray([0.0], dtype=np.float64),
                bounds_error=False,
                fill_value=0.0,
                assume_sorted=True,
            )
            return self.burn_area_interpolator

        # Lift the inhibited casing above every iso level so marching cubes
        # meshes only the burning front, not the wall.
        distance_field = np.ma.filled(regression_map, max_iso_level + 1.0)
        axial_grid_spacing = self.get_axial_grid_spacing()
        radial_grid_spacing = self.get_radial_grid_spacing()
        grid_spacing = (
            axial_grid_spacing,
            radial_grid_spacing,
            radial_grid_spacing,
        )

        # The iso-surface at level 0 lies on the voxelized initial face and is
        # degenerate, so sample above it and hold the first area back to web 0.
        iso_levels = np.linspace(
            max_iso_level / ISO_LEVEL_COUNT, max_iso_level, ISO_LEVEL_COUNT
        )
        iso_surface_areas = _compute_iso_surface_areas(
            distance_field, iso_levels, grid_spacing
        )
        web_distances_denormalized = np.concatenate(
            ([0.0], np.asarray(self.denormalize(iso_levels), dtype=np.float64))
        )
        iso_surface_areas = np.concatenate(
            ([iso_surface_areas[0]], iso_surface_areas)
        )

        burn_area_smoothed = filters.smooth_savitzky_golay(
            iso_surface_areas, window_length=9, polyorder=3
        )
        self.burn_area_interpolator = interp1d(
            web_distances_denormalized,
            burn_area_smoothed,
            bounds_error=False,
            fill_value=(
                float(burn_area_smoothed[0]),
                float(burn_area_smoothed[-1]),
            ),  # type: ignore[arg-type]
            assume_sorted=True,
        )

    return self.burn_area_interpolator

get_center_of_gravity(web_distance)

Return the center of gravity of a 3D FMM grain segment.

Parameters:

Name Type Description Default
web_distance float

Web distance traveled [m].

 required

Returns:

Type Description
NDArray[float64]

Center of gravity [z, x, y] in meters from the port of the segment.

Raises:

Type Description
GrainGeometryError

If web distance exceeds the web thickness, or if no active material is found.
Source code in machwave/models/grain/fmm/_3d.py 
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
def get_center_of_gravity(self, web_distance: float) -> NDArray[np.float64]:
    """
    Return the center of gravity of a 3D FMM grain segment.

    Args:
        web_distance: Web distance traveled [m].

    Returns:
        Center of gravity [z, x, y] in meters from the port of the segment.

    Raises:
        GrainGeometryError: If web distance exceeds the web thickness, or
            if no active material is found.
    """
    self._validate_web_distance(web_distance)
    moments = self._solid_mask_moments(web_distance)
    if moments is None:
        raise grain.GrainGeometryError(
            "No active material found at the given web distance."
        )
    centroid, _ = moments
    # Copy so callers cannot mutate the per-web cached array.
    return centroid.copy()

get_contours(web_distance, axial_position_normalized)

Return the contours of one axial slice at a given web distance.

Parameters:

Name Type Description Default
web_distance float

Distance traveled into the grain web [m].

 required
axial_position_normalized float

Axial position as a fraction of the segment length, measured from the aft (nozzle) end.

 required

Returns:

Type Description
list[NDArray[float64]]

Contour arrays of the slice, each an (N, 2) array of points.
Source code in machwave/models/grain/fmm/_3d.py 
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
def get_contours(
    self, web_distance: float, axial_position_normalized: float
) -> list[NDArray[np.float64]]:
    """
    Return the contours of one axial slice at a given web distance.

    Args:
        web_distance: Distance traveled into the grain web [m].
        axial_position_normalized: Axial position as a fraction of the
            segment length, measured from the aft (nozzle) end.

    Returns:
        Contour arrays of the slice, each an `(N, 2)` array of points.
    """
    iso_level = self.normalize(web_distance)
    axial_slice = self.get_regression_map()[
        self.get_axial_index(axial_position_normalized)
    ]
    return fmm_contours.get_iso_contours(axial_slice, iso_level)

get_minimum_port_area(web_distance)

Return the smallest port area along the segment [m^2].

The cross section of a 3D segment varies along its length, so the flow passes a bottleneck: the station holding the most solid propellant.

Parameters:

Name Type Description Default
web_distance float

The distance traveled into the grain web.

 required

Returns:

Type Description
float

Smallest port area along the segment [m^2].
Source code in machwave/models/grain/fmm/_3d.py 
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
def get_minimum_port_area(self, web_distance: float) -> float:
    """
    Return the smallest port area along the segment [m^2].

    The cross section of a 3D segment varies along its length, so the flow
    passes a bottleneck: the station holding the most solid propellant.

    Args:
        web_distance: The distance traveled into the grain web.

    Returns:
        Smallest port area along the segment [m^2].
    """
    iso_level = self.normalize(web_distance)
    valid = np.logical_not(self.get_outer_diameter_mask())
    solid = np.logical_and(self.get_regression_map() > iso_level, valid)

    solid_cells_per_slice = np.count_nonzero(solid, axis=(1, 2))
    face_area = float(
        self.cells_to_square_meters(float(solid_cells_per_slice.max()))
    )
    return geometric.get_circle_area(self.outer_diameter) - face_area

get_moment_of_inertia(ideal_density, web_distance=0.0)

Calculate the moment of inertia tensor at the segment's its center of gravity.

Parameters:

Name Type Description Default
ideal_density float

Propellant ideal density [kg/m^3].

 required
web_distance float

Web distance traveled [m].

 0.0

Returns:

Type Description
NDArray[float64]

A 3x3 inertia tensor [kg-m^2].

Raises:

Type Description
GrainGeometryError

If web distance exceeds web thickness or if no active material is found.
Source code in machwave/models/grain/fmm/_3d.py 
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
def get_moment_of_inertia(
    self, ideal_density: float, web_distance: float = 0.0
) -> NDArray[np.float64]:
    """
    Calculate the moment of inertia tensor at the segment's its center of gravity.

    Args:
        ideal_density: Propellant ideal density [kg/m^3].
        web_distance: Web distance traveled [m].

    Returns:
        A 3x3 inertia tensor [kg-m^2].

    Raises:
        GrainGeometryError: If web distance exceeds web thickness or
            if no active material is found.
    """
    self._validate_web_distance(web_distance)
    moments = self._solid_mask_moments(web_distance)
    if moments is None:
        raise grain.GrainGeometryError(
            "No active material found at the given web distance."
        )
    _, central_second_moments = moments

    element_mass = self.get_voxel_volume() * ideal_density * self.density_ratio
    return mechanics.get_moment_of_inertia_tensor_from_central_moments(
        central_second_moments, element_mass
    )

get_padded_regression_map()

Return the regression map including the void slice past each exposed end.

Source code in machwave/models/grain/fmm/_3d.py 
309
310
311
312
313
314
def get_padded_regression_map(self) -> np.ndarray:
    """Return the regression map including the void slice past each exposed end."""
    regression_map = self.get_regression_map()
    if self.padded_regression_map is None:
        return regression_map
    return self.padded_regression_map

get_port_area(web_distance, z=0.0)

Calculates the port area at a given web distance and axial height z.

This method extracts a single 2D slice from the 3D face map by converting the physical height z into an integer index, and then computes the port area for that slice.

Parameters:

Name Type Description Default
web_distance float

The distance traveled into the grain web.

 required
z float

Axial position along the grain [m], measured from the aft (nozzle) end. Defaults to nozzle end.

 0.0

Returns:

Type Description
float

A float representing the port area at the specified z slice, in m^2.
Source code in machwave/models/grain/fmm/_3d.py 
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
def get_port_area(self, web_distance: float, z: float = 0.0) -> float:
    """
    Calculates the port area at a given web distance and axial height z.

    This method extracts a single 2D slice from the 3D face map by converting
    the physical height z into an integer index, and then computes the port
    area for that slice.

    Args:
        web_distance: The distance traveled into the grain web.
        z: Axial position along the grain [m], measured from the aft
            (nozzle) end. Defaults to nozzle end.

    Returns:
        A float representing the port area at the specified z slice, in m^2.
    """
    iso_level = self.normalize(web_distance)
    valid = np.logical_not(self.get_outer_diameter_mask())
    solid = np.logical_and(self.get_regression_map() > iso_level, valid)

    axial_index = self.get_axial_index(z / self.length)

    face_area = float(
        self.cells_to_square_meters(float(np.count_nonzero(solid[axial_index])))
    )
    return geometric.get_circle_area(self.outer_diameter) - face_area

get_volume_interpolator()

Return a cached interpolator for volume [m^3] for a web distance.

Source code in machwave/models/grain/fmm/_3d.py 
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
def get_volume_interpolator(self) -> Callable[[float], float]:
    """Return a cached interpolator for volume [m^3] for a web distance."""
    if self.volume_interpolator is None:
        regression_map = self.get_regression_map()
        regression_distances_sorted = np.asarray(
            regression_map[~np.ma.getmaskarray(regression_map)], dtype=np.float64
        ).ravel()
        regression_distances_sorted.sort()
        max_regression_distance = (
            float(regression_distances_sorted[-1])
            if regression_distances_sorted.size
            else 0.0
        )
        if max_regression_distance <= 0.0:
            self.volume_interpolator = interp1d(
                np.asarray([0.0], dtype=np.float64),
                np.asarray([0.0], dtype=np.float64),
                bounds_error=False,
                fill_value=0.0,
                assume_sorted=True,
            )
            return self.volume_interpolator

        iso_level_count = int(max_regression_distance * self.grid_resolution) + 2
        iso_levels_normalized = (
            np.arange(iso_level_count, dtype=np.float64) / self.grid_resolution
        )

        count_at_or_below_level = np.searchsorted(
            regression_distances_sorted, iso_levels_normalized, side="right"
        )
        solid_voxel_count = float(
            regression_distances_sorted.size
        ) - count_at_or_below_level.astype(np.float64)
        volume_per_iso_level = solid_voxel_count * self.get_voxel_volume()

        self.volume_interpolator = interp1d(
            iso_levels_normalized,
            volume_per_iso_level,
            bounds_error=False,
            fill_value=(
                float(volume_per_iso_level[0]),
                float(volume_per_iso_level[-1]),
            ),  # type: ignore[arg-type]
            assume_sorted=True,
        )

    return self.volume_interpolator

get_voxel_volume()

Volume of one cell of the anisotropic grid [m^3].

Source code in machwave/models/grain/fmm/_3d.py 
469
470
471
def get_voxel_volume(self) -> float:
    """Volume of one cell of the anisotropic grid [m^3]."""
    return self.get_radial_grid_spacing() ** 2 * self.get_axial_grid_spacing()

FMMSTLGrainSegment

Bases: FMMGrainSegment3D, ABC

FMM grain segment loaded from an STL mesh.

Source code in machwave/models/grain/fmm/stl.py 
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
class FMMSTLGrainSegment(grain_fmm.FMMGrainSegment3D, ABC):
    """FMM grain segment loaded from an STL mesh."""

    def __init__(
        self,
        file_path: str,
        outer_diameter: float,
        length: float,
        inhibited_surfaces: grain_base.InhibitedSurfaces | None = None,
        grid_resolution: int = grain_fmm.DEFAULT_GRID_RESOLUTION,
    ) -> None:
        """
        Initialize an STL-backed FMM grain segment.

        Args:
            file_path: Path to a watertight STL mesh of the grain.
            outer_diameter: Outer diameter [m].
            length: Segment length [m].
            inhibited_surfaces: Surfaces inhibited from burning.
            grid_resolution: Grid points per axis of the cross-section.
        """
        self.file_path = file_path
        self.outer_diameter = outer_diameter
        self.length = length

        # "Cache" variables:
        self.face_area_interpolator = None

        super().__init__(
            length=length,
            outer_diameter=outer_diameter,
            inhibited_surfaces=inhibited_surfaces,
            grid_resolution=grid_resolution,
        )

    def validate(self) -> None:
        """Validate STL grain segment geometry."""
        grain.GrainSegment.validate(self)
        if not self.grid_resolution >= 20:
            raise grain.GrainGeometryError(
                f"Grid resolution must be at least 20 for STL grains, got {self.grid_resolution}"
            )
        self._validate_axial_resolution()

    def get_grid_spacing(self) -> float:
        """Return the grid spacing [m] used to voxelize the mesh."""
        return self.outer_diameter / (self.grid_resolution - 1)

    def generate_initial_face_map(self) -> np.typing.NDArray[np.int_]:
        """
        Generate the initial face map by voxelizing an STL mesh.

        The trimesh voxel grid does not line up with the FMM grid, so it is
        resampled onto the canonical `(normalized_length, grid_resolution, grid_resolution)`
        shape the rest of the 3D machinery expects.
        """
        mesh = trimesh.load_mesh(self.file_path)
        assert isinstance(mesh, trimesh.Trimesh), "Expected a single Trimesh"
        assert mesh.is_watertight, "Mesh must be watertight"

        voxels = mesh.voxelized(pitch=self.get_grid_spacing())
        assert voxels is not None, "Voxelization failed"
        voxel_map = voxels.fill().matrix.view(np.ndarray).transpose().astype(np.int_)

        return _resample_nearest(voxel_map, self.get_coordinate_grids()[0].shape)

__init__(file_path, outer_diameter, length, inhibited_surfaces=None, grid_resolution=grain_fmm.DEFAULT_GRID_RESOLUTION)

Initialize an STL-backed FMM grain segment.

Parameters:

Name Type Description Default
file_path str

Path to a watertight STL mesh of the grain.

 required
outer_diameter float

Outer diameter [m].

 required
length float

Segment length [m].

 required
inhibited_surfaces InhibitedSurfaces | None

Surfaces inhibited from burning.

 None
grid_resolution int

Grid points per axis of the cross-section.

 DEFAULT_GRID_RESOLUTION
Source code in machwave/models/grain/fmm/stl.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
def __init__(
    self,
    file_path: str,
    outer_diameter: float,
    length: float,
    inhibited_surfaces: grain_base.InhibitedSurfaces | None = None,
    grid_resolution: int = grain_fmm.DEFAULT_GRID_RESOLUTION,
) -> None:
    """
    Initialize an STL-backed FMM grain segment.

    Args:
        file_path: Path to a watertight STL mesh of the grain.
        outer_diameter: Outer diameter [m].
        length: Segment length [m].
        inhibited_surfaces: Surfaces inhibited from burning.
        grid_resolution: Grid points per axis of the cross-section.
    """
    self.file_path = file_path
    self.outer_diameter = outer_diameter
    self.length = length

    # "Cache" variables:
    self.face_area_interpolator = None

    super().__init__(
        length=length,
        outer_diameter=outer_diameter,
        inhibited_surfaces=inhibited_surfaces,
        grid_resolution=grid_resolution,
    )

generate_initial_face_map()

Generate the initial face map by voxelizing an STL mesh.

The trimesh voxel grid does not line up with the FMM grid, so it is resampled onto the canonical (normalized_length, grid_resolution, grid_resolution) shape the rest of the 3D machinery expects.

Source code in machwave/models/grain/fmm/stl.py 
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
def generate_initial_face_map(self) -> np.typing.NDArray[np.int_]:
    """
    Generate the initial face map by voxelizing an STL mesh.

    The trimesh voxel grid does not line up with the FMM grid, so it is
    resampled onto the canonical `(normalized_length, grid_resolution, grid_resolution)`
    shape the rest of the 3D machinery expects.
    """
    mesh = trimesh.load_mesh(self.file_path)
    assert isinstance(mesh, trimesh.Trimesh), "Expected a single Trimesh"
    assert mesh.is_watertight, "Mesh must be watertight"

    voxels = mesh.voxelized(pitch=self.get_grid_spacing())
    assert voxels is not None, "Voxelization failed"
    voxel_map = voxels.fill().matrix.view(np.ndarray).transpose().astype(np.int_)

    return _resample_nearest(voxel_map, self.get_coordinate_grids()[0].shape)

get_grid_spacing()

Return the grid spacing [m] used to voxelize the mesh.

Source code in machwave/models/grain/fmm/stl.py 
67
68
69
def get_grid_spacing(self) -> float:
    """Return the grid spacing [m] used to voxelize the mesh."""
    return self.outer_diameter / (self.grid_resolution - 1)

validate()

Validate STL grain segment geometry.

Source code in machwave/models/grain/fmm/stl.py 
58
59
60
61
62
63
64
65
def validate(self) -> None:
    """Validate STL grain segment geometry."""
    grain.GrainSegment.validate(self)
    if not self.grid_resolution >= 20:
        raise grain.GrainGeometryError(
            f"Grid resolution must be at least 20 for STL grains, got {self.grid_resolution}"
        )
    self._validate_axial_resolution()