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), instead of re-meshing the geometry at every regression step.

  • 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
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 cells_to_meters(
        self, value: float | NDArray[np.float64]
    ) -> float | NDArray[np.float64]:
        """
        Convert a pixel-distance value to [m].

        Args:
            value: Distance in pixel units.

        Returns:
            Distance [m].
        """
        return self.outer_diameter * (value / self.grid_resolution)

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

        Args:
            value: Area in pixel units.

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

    def get_normalized_spacing(self) -> float:
        """Return the cell size in normalized coordinates (`1 / grid_resolution`)."""
        return 1 / self.grid_resolution

    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."""
        return skfmm.distance(masked_face, dx=self.get_normalized_spacing()) * 2

    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 pixel-distance value to [m].

Parameters:

Name Type Description Default
value float | NDArray[float64]

Distance in pixel units.

 required

Returns:

Type Description
float | NDArray[float64]

Distance [m].
Source code in machwave/models/grain/fmm/base.py 
112
113
114
115
116
117
118
119
120
121
122
123
124
def cells_to_meters(
    self, value: float | NDArray[np.float64]
) -> float | NDArray[np.float64]:
    """
    Convert a pixel-distance value to [m].

    Args:
        value: Distance in pixel units.

    Returns:
        Distance [m].
    """
    return self.outer_diameter * (value / self.grid_resolution)

cells_to_square_meters(value)

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

Parameters:

Name Type Description Default
value float | NDArray[float64]

Area in pixel units.

 required

Returns:

Type Description
float | NDArray[float64]

Area [m^2].
Source code in machwave/models/grain/fmm/base.py 
126
127
128
129
130
131
132
133
134
135
136
137
138
def cells_to_square_meters(
    self, value: float | NDArray[np.float64]
) -> float | NDArray[np.float64]:
    """
    Convert a pixel-area value to [m^2].

    Args:
        value: Area in pixel units.

    Returns:
        Area [m^2].
    """
    return (self.outer_diameter**2) * (value / (self.grid_resolution**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 
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
@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 
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
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 
171
172
173
174
175
176
177
178
179
180
181
182
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_normalized_spacing()

Return the cell size in normalized coordinates (1 / grid_resolution).

Source code in machwave/models/grain/fmm/base.py 
140
141
142
def get_normalized_spacing(self) -> float:
    """Return the cell size in normalized coordinates (`1 / grid_resolution`)."""
    return 1 / self.grid_resolution

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_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 
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
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 
227
228
229
230
231
232
233
234
235
236
237
238
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
390
391
392
393
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, self.grid_resolution)
                        )
                        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, self.grid_resolution)
                )
                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 / 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 
221
222
223
224
225
226
227
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
218
219
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, self.grid_resolution)
                    )
                    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 
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
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 
241
242
243
244
245
246
247
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 
229
230
231
232
233
234
235
236
237
238
239
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, self.grid_resolution)
            )
            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 
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
def 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 
 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
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
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

        # 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 int(self.grid_resolution * self.length / self.outer_diameter)

    def get_coordinate_grids(
        self,
    ) -> tuple[
        NDArray[np.float64],
        NDArray[np.float64],
        NDArray[np.float64],
    ]:
        if self.coordinate_grids is None:
            map_y, map_z, map_x = np.meshgrid(
                np.linspace(-1, 1, self.grid_resolution),
                np.linspace(1, 0, self.get_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_]]:
        if face_map.shape[0] <= 2:
            return face_map, excluded_mask

        inner_surface_inhibited_cells = np.zeros_like(face_map, dtype=bool)
        inner_surface_inhibited_cells[1:-1] = face_map[1:-1] == 0
        inner_surface_inhibited_cells[0] = inner_surface_inhibited_cells[1]
        inner_surface_inhibited_cells[-1] = inner_surface_inhibited_cells[-2]

        # super() runs binary_erosion on the full 3D volume, which includes the
        # end-face layers in the outer_surface_boundary. Apply end inhibition AFTER so
        # those cells are not overwritten back to 0 by the outer-surface logic.
        face_map, excluded_mask = super()._apply_surface_inhibition(
            face_map, excluded_mask
        )

        if self.inhibited_surfaces.upper_end:
            end_face_inhibited_cells = (
                face_map[-1] == 0
            ) & ~inner_surface_inhibited_cells[-1]
            face_map[-1][end_face_inhibited_cells] = 1

        if self.inhibited_surfaces.lower_end:
            end_face_inhibited_cells = (
                face_map[0] == 0
            ) & ~inner_surface_inhibited_cells[0]
            face_map[0][end_face_inhibited_cells] = 1

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

        return face_map, excluded_mask

    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.length / max(self.get_axial_resolution() - 1, 1)
        radial_grid_spacing = self.outer_diameter / (self.grid_resolution - 1)
        distance = skfmm.distance(
            masked_face,
            dx=[axial_grid_spacing, radial_grid_spacing, radial_grid_spacing],  # type: ignore[arg-type]
        )
        return distance * (2.0 / self.outer_diameter)

    def get_contours(
        self, web_distance: float, axial_position_normalized: float
    ) -> list[NDArray[np.float64]]:
        iso_level = self.normalize(web_distance)
        axial_index = int(round(axial_position_normalized))
        axial_slice = self.get_regression_map()[axial_index]
        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_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.length / max(self.get_axial_resolution() - 1, 1)
            radial_grid_spacing = self.outer_diameter / (self.grid_resolution - 1)
            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_level_count = 80
            iso_levels = np.linspace(
                max_iso_level / iso_level_count, max_iso_level, iso_level_count
            )
            iso_surface_areas = np.array(
                [
                    self._compute_iso_surface_area(
                        distance_field, float(level), grid_spacing
                    )
                    for level in iso_levels
                ],
                dtype=np.float64,
            )
            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

    @staticmethod
    def _compute_iso_surface_area(
        distance_field: NDArray[np.float64],
        level: float,
        grid_spacing: tuple[float, float, float],
    ) -> float:
        """Marching-cubes area [m^2] of one regression iso-level, 0 if empty."""
        try:
            vertices, faces, _, _ = measure.marching_cubes(
                distance_field, level=level, spacing=grid_spacing
            )
        except (ValueError, RuntimeError):
            return 0.0
        return float(measure.mesh_surface_area(vertices, faces))

    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)

        normalized_z = z / self.length
        max_index = self.get_axial_resolution() - 1
        axial_index = int(round(normalized_z * max_index))
        axial_index = (
            0
            if axial_index < 0
            else (max_index if axial_index > max_index else axial_index)
        )

        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_voxel_volume(self) -> float:
        return (float(self.denormalize(self.get_normalized_spacing())) * 2) ** 3

    def get_volume(self, web_distance: float) -> float:
        solid_voxel_count = int(np.count_nonzero(self._get_solid_mask(web_distance)))
        return solid_voxel_count * self.get_voxel_volume()

    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 / 2
        x = np.asarray(
            self.cells_to_meters(np.arange(x_count) - grid_center), dtype=np.float64
        )
        y = np.asarray(
            self.cells_to_meters(np.arange(y_count) - grid_center), dtype=np.float64
        )
        z = np.asarray(self.cells_to_meters(np.arange(axial_count)), dtype=np.float64)

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

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 
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
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_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.length / max(self.get_axial_resolution() - 1, 1)
        radial_grid_spacing = self.outer_diameter / (self.grid_resolution - 1)
        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_level_count = 80
        iso_levels = np.linspace(
            max_iso_level / iso_level_count, max_iso_level, iso_level_count
        )
        iso_surface_areas = np.array(
            [
                self._compute_iso_surface_area(
                    distance_field, float(level), grid_spacing
                )
                for level in iso_levels
            ],
            dtype=np.float64,
        )
        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 
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
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_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 
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
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_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 
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
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)

    normalized_z = z / self.length
    max_index = self.get_axial_resolution() - 1
    axial_index = int(round(normalized_z * max_index))
    axial_index = (
        0
        if axial_index < 0
        else (max_index if axial_index > max_index else axial_index)
    )

    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

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