Skip to content

models.propulsion.grain.fmm

machwave.models.propulsion.grain.fmm

This module contains classes that model grain geometries using the Fast Marching Method (FMM). Both 2D and 3D geometries are supported, along with STL files.

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.

This class was inspired by the Andrew Reilley's software openMotor, in particular the fmm module. openMotor's repository can be accessed at: https://github.com/reilleya/openMotor

Source code in machwave/models/propulsion/grain/fmm/base.py 
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 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
class FMMGrainSegment(GrainSegment, ABC):
    """
    Fast Marching Method (FMM) implementation of a grain segment.

    This class was inspired by the Andrew Reilley's software openMotor, in
    particular the fmm module.
    openMotor's repository can be accessed at:
    https://github.com/reilleya/openMotor
    """

    def __init__(
        self,
        map_dim: int,
        length: float,
        outer_diameter: float,
        inhibited_ends: int = 0,
        density_ratio: float = 1.0,
    ) -> None:
        self.map_dim = map_dim

        # "Cache" variables:
        self.maps = None
        self.mask = None
        self.masked_face = None
        self.regression_map = None
        self.web_thickness = None

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

    @abstractmethod
    def get_initial_face_map(self) -> np.typing.NDArray[np.int_]:
        """
        Method needs to be implemented for each and every geometry.
        """
        pass

    @abstractmethod
    def get_maps(self) -> tuple:
        """
        Returns:
            - 2D: (map_x, map_y)
            - 3D: (map_x, map_y, map_z)
        """
        pass

    @abstractmethod
    def get_mask(self) -> np.ndarray:
        """
        Implementation varies depending if the geometry is 2D or 3D.
        """
        pass

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

        This method ensures the grain map dimension meets the minimum
        required size. If validation fails, a GrainGeometryError is raised.

        Raises:
            GrainGeometryError: If the grain map dimension is below the valid threshold.
        """
        super().validate()
        if not self.map_dim >= MINIMUM_MAP_DIMENSION:
            raise GrainGeometryError(
                f"Map dimension must be at least {MINIMUM_MAP_DIMENSION}, "
                f"got {self.map_dim}"
            )

    def normalize(self, value: int | float) -> float:
        """
        Converts a raw dimensional value into a normalized scale based on the
        object's outer diameter.

        Args:
            value: The dimensional value (e.g., length) to normalize.

        Returns:
            A float representing the dimension as a fraction of the object's
            half-diameter.
        """
        return value / (0.5 * self.outer_diameter)

    def denormalize(
        self, value: int | float | NDArray[np.float64]
    ) -> float | NDArray[np.float64]:
        """
        Converts a normalized input value into an actual dimension based on the
        object's outer diameter.

        Args:
            value: A numeric value representing a normalized quantity.

        Returns:
            The denormalized value as a float, calculated by scaling `value` with
            the object's outer diameter.
        """
        return (value / 2) * (self.outer_diameter)

    def map_to_area(
        self, value: float | NDArray[np.float64]
    ) -> float | NDArray[np.float64]:
        """
        Convert a pixel-area value into square meters.

        The conversion is based on the ratio of this object's outer diameter
        (squared) to the total pixel map dimension (squared).

        Args:
            value: The area in pixel units.

        Returns:
            The corresponding area in square meters.
        """
        return (self.outer_diameter**2) * (value / (self.map_dim**2))

    def map_to_length(
        self, value: float | NDArray[np.float64]
    ) -> float | NDArray[np.float64]:
        """
        Convert a pixel-distance value into meters.

        The conversion is based on the ratio of this object's outer diameter
        to its total map dimension.

        Args:
            value: The distance in pixel units.

        Returns:
            The corresponding distance in meters.
        """
        return self.outer_diameter * (value / self.map_dim)

    def get_empty_face_map(self) -> np.ndarray:
        """
        Return a new face map consisting entirely of ones.

        The shape of the array matches the first element in the object's stored maps.
        """
        return np.ones_like(self.get_maps()[0])

    def get_masked_face(self) -> np.ndarray:
        """
        Return a masked representation of the face map.

        The mask is circular and normalized to the map dimensions. If a mask
        has not been created yet, it is generated by combining the initial face
        map with the circular mask.
        """
        if self.masked_face is None:
            self.masked_face = np.ma.MaskedArray(
                self.get_initial_face_map(), self.get_mask()
            )
        return self.masked_face

    def get_cell_size(self) -> float:
        """
        Return the size of each grid cell in normalized coordinates.

        The value is derived by taking 1 divided by the map dimension.
        """
        return 1 / self.map_dim

    def get_regression_map(self):
        """
        Calculate and return the distance map for grain regression.

        This uses the fast marching method (scikit-fmm) on the masked face.
        Each value represents the distance from the initial face along the
        cross-section of the grain.
        """
        if self.regression_map is None:
            self.regression_map = (
                skfmm.distance(self.get_masked_face(), dx=self.get_cell_size()) * 2
            )
        return self.regression_map

    def get_web_thickness(self) -> float:
        """
        Return the maximum thickness of the grain web in real units.

        The web thickness is the largest distance from the center of the
        grain segment, derived from the distance map and converted to a
        real-world measurement.
        """
        if self.web_thickness is None:
            self.web_thickness = float(
                self.denormalize(float(np.amax(self.get_regression_map())))
            )
        return float(self.web_thickness)

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

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

        Args:
            web_distance: The depth of regression into the grain web.
            *args: Additional positional arguments.
            **kwargs: Additional keyword arguments.

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

    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()
        valid = np.logical_not(self.get_mask())

        # Create a masked array, where ~valid cells are masked out
        maskarr = np.ma.MaskedArray(
            (regression_map > web_distance_normalized).astype(np.int64), mask=~valid
        )

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

denormalize(value)

Converts a normalized input value into an actual dimension based on the object's outer diameter.

Parameters:

Name Type Description Default
value int | float | NDArray[float64]

A numeric value representing a normalized quantity.

 required

Returns:

Type Description
float | NDArray[float64]

The denormalized value as a float, calculated by scaling value with
float | NDArray[float64]

the object's outer diameter.
Source code in machwave/models/propulsion/grain/fmm/base.py 
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
def denormalize(
    self, value: int | float | NDArray[np.float64]
) -> float | NDArray[np.float64]:
    """
    Converts a normalized input value into an actual dimension based on the
    object's outer diameter.

    Args:
        value: A numeric value representing a normalized quantity.

    Returns:
        The denormalized value as a float, calculated by scaling `value` with
        the object's outer diameter.
    """
    return (value / 2) * (self.outer_diameter)

get_cell_size()

Return the size of each grid cell in normalized coordinates.

The value is derived by taking 1 divided by the map dimension.

Source code in machwave/models/propulsion/grain/fmm/base.py 
172
173
174
175
176
177
178
def get_cell_size(self) -> float:
    """
    Return the size of each grid cell in normalized coordinates.

    The value is derived by taking 1 divided by the map dimension.
    """
    return 1 / self.map_dim

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

Return the contours of the regression map after a specified web distance.

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

The depth of regression into the grain web.

 required
*args

Additional positional arguments.

 ()
**kwargs

Additional keyword arguments.

 {}

Returns:

Type Description
list[NDArray[float64]]

An array representing the computed contours of the grain regression.
Source code in machwave/models/propulsion/grain/fmm/base.py 
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
@abstractmethod
def get_contours(
    self, web_distance: float, *args, **kwargs
) -> list[NDArray[np.float64]]:
    """
    Return the contours of the regression map after a specified web distance.

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

    Args:
        web_distance: The depth of regression into the grain web.
        *args: Additional positional arguments.
        **kwargs: Additional keyword arguments.

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

get_empty_face_map()

Return a new face map consisting entirely of ones.

The shape of the array matches the first element in the object's stored maps.

Source code in machwave/models/propulsion/grain/fmm/base.py 
150
151
152
153
154
155
156
def get_empty_face_map(self) -> np.ndarray:
    """
    Return a new face map consisting entirely of ones.

    The shape of the array matches the first element in the object's stored maps.
    """
    return np.ones_like(self.get_maps()[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/propulsion/grain/fmm/base.py 
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
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()
    valid = np.logical_not(self.get_mask())

    # Create a masked array, where ~valid cells are masked out
    maskarr = np.ma.MaskedArray(
        (regression_map > web_distance_normalized).astype(np.int64), mask=~valid
    )

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

get_initial_face_map() abstractmethod

Method needs to be implemented for each and every geometry.

Source code in machwave/models/propulsion/grain/fmm/base.py 
46
47
48
49
50
51
@abstractmethod
def get_initial_face_map(self) -> np.typing.NDArray[np.int_]:
    """
    Method needs to be implemented for each and every geometry.
    """
    pass

get_maps() abstractmethod

Returns:

Type Description
tuple
  • 2D: (map_x, map_y)
tuple
  • 3D: (map_x, map_y, map_z)
Source code in machwave/models/propulsion/grain/fmm/base.py 
53
54
55
56
57
58
59
60
@abstractmethod
def get_maps(self) -> tuple:
    """
    Returns:
        - 2D: (map_x, map_y)
        - 3D: (map_x, map_y, map_z)
    """
    pass

get_mask() abstractmethod

Implementation varies depending if the geometry is 2D or 3D.

Source code in machwave/models/propulsion/grain/fmm/base.py 
62
63
64
65
66
67
@abstractmethod
def get_mask(self) -> np.ndarray:
    """
    Implementation varies depending if the geometry is 2D or 3D.
    """
    pass

get_masked_face()

Return a masked representation of the face map.

The mask is circular and normalized to the map dimensions. If a mask has not been created yet, it is generated by combining the initial face map with the circular mask.

Source code in machwave/models/propulsion/grain/fmm/base.py 
158
159
160
161
162
163
164
165
166
167
168
169
170
def get_masked_face(self) -> np.ndarray:
    """
    Return a masked representation of the face map.

    The mask is circular and normalized to the map dimensions. If a mask
    has not been created yet, it is generated by combining the initial face
    map with the circular mask.
    """
    if self.masked_face is None:
        self.masked_face = np.ma.MaskedArray(
            self.get_initial_face_map(), self.get_mask()
        )
    return self.masked_face

get_regression_map()

Calculate and return the distance map for grain regression.

This uses the fast marching method (scikit-fmm) on the masked face. Each value represents the distance from the initial face along the cross-section of the grain.

Source code in machwave/models/propulsion/grain/fmm/base.py 
180
181
182
183
184
185
186
187
188
189
190
191
192
def get_regression_map(self):
    """
    Calculate and return the distance map for grain regression.

    This uses the fast marching method (scikit-fmm) on the masked face.
    Each value represents the distance from the initial face along the
    cross-section of the grain.
    """
    if self.regression_map is None:
        self.regression_map = (
            skfmm.distance(self.get_masked_face(), dx=self.get_cell_size()) * 2
        )
    return self.regression_map

get_web_thickness()

Return the maximum thickness of the grain web in real units.

The web thickness is the largest distance from the center of the grain segment, derived from the distance map and converted to a real-world measurement.

Source code in machwave/models/propulsion/grain/fmm/base.py 
194
195
196
197
198
199
200
201
202
203
204
205
206
def get_web_thickness(self) -> float:
    """
    Return the maximum thickness of the grain web in real units.

    The web thickness is the largest distance from the center of the
    grain segment, derived from the distance map and converted to a
    real-world measurement.
    """
    if self.web_thickness is None:
        self.web_thickness = float(
            self.denormalize(float(np.amax(self.get_regression_map())))
        )
    return float(self.web_thickness)

map_to_area(value)

Convert a pixel-area value into square meters.

The conversion is based on the ratio of this object's outer diameter (squared) to the total pixel map dimension (squared).

Parameters:

Name Type Description Default
value float | NDArray[float64]

The area in pixel units.

 required

Returns:

Type Description
float | NDArray[float64]

The corresponding area in square meters.
Source code in machwave/models/propulsion/grain/fmm/base.py 
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def map_to_area(
    self, value: float | NDArray[np.float64]
) -> float | NDArray[np.float64]:
    """
    Convert a pixel-area value into square meters.

    The conversion is based on the ratio of this object's outer diameter
    (squared) to the total pixel map dimension (squared).

    Args:
        value: The area in pixel units.

    Returns:
        The corresponding area in square meters.
    """
    return (self.outer_diameter**2) * (value / (self.map_dim**2))

map_to_length(value)

Convert a pixel-distance value into meters.

The conversion is based on the ratio of this object's outer diameter to its total map dimension.

Parameters:

Name Type Description Default
value float | NDArray[float64]

The distance in pixel units.

 required

Returns:

Type Description
float | NDArray[float64]

The corresponding distance in meters.
Source code in machwave/models/propulsion/grain/fmm/base.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
def map_to_length(
    self, value: float | NDArray[np.float64]
) -> float | NDArray[np.float64]:
    """
    Convert a pixel-distance value into meters.

    The conversion is based on the ratio of this object's outer diameter
    to its total map dimension.

    Args:
        value: The distance in pixel units.

    Returns:
        The corresponding distance in meters.
    """
    return self.outer_diameter * (value / self.map_dim)

normalize(value)

Converts a raw dimensional value into a normalized scale based on the object's outer diameter.

Parameters:

Name Type Description Default
value int | float

The dimensional value (e.g., length) to normalize.

 required

Returns:

Type Description
float

A float representing the dimension as a fraction of the object's
float

half-diameter.
Source code in machwave/models/propulsion/grain/fmm/base.py 
86
87
88
89
90
91
92
93
94
95
96
97
98
def normalize(self, value: int | float) -> float:
    """
    Converts a raw dimensional value into a normalized scale based on the
    object's outer diameter.

    Args:
        value: The dimensional value (e.g., length) to normalize.

    Returns:
        A float representing the dimension as a fraction of the object's
        half-diameter.
    """
    return value / (0.5 * self.outer_diameter)

validate()

Validates the internal geometry of the grain.

This method ensures the grain map dimension meets the minimum required size. If validation fails, a GrainGeometryError is raised.

Raises:

Type Description
GrainGeometryError

If the grain map dimension is below the valid threshold.
Source code in machwave/models/propulsion/grain/fmm/base.py 
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
def validate(self) -> None:
    """
    Validates the internal geometry of the grain.

    This method ensures the grain map dimension meets the minimum
    required size. If validation fails, a GrainGeometryError is raised.

    Raises:
        GrainGeometryError: If the grain map dimension is below the valid threshold.
    """
    super().validate()
    if not self.map_dim >= MINIMUM_MAP_DIMENSION:
        raise GrainGeometryError(
            f"Map dimension must be at least {MINIMUM_MAP_DIMENSION}, "
            f"got {self.map_dim}"
        )

FMMGrainSegment2D

Bases: FMMGrainSegment, GrainSegment2D, ABC

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

This class was inspired by Andrew Reilley's openMotor software, in particular the fmm module. See: https://github.com/reilleya/openMotor

Source code in machwave/models/propulsion/grain/fmm/_2d.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
 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
class FMMGrainSegment2D(FMMGrainSegment, GrainSegment2D, ABC):
    """
    Fast Marching Method (FMM) implementation for 2D grain segment.

    This class was inspired by Andrew Reilley's openMotor software,
    in particular the fmm module. See:
    https://github.com/reilleya/openMotor
    """

    def __init__(
        self,
        length: float,
        outer_diameter: float,
        inhibited_ends: int = 0,
        map_dim: int = 1000,
        density_ratio: float = 1.0,
    ) -> None:
        self.face_area_interp_func: Callable[[float], float] | None = None
        self.burn_area_interp_func: Callable[[float], float] | None = None
        super().__init__(
            length=length,
            outer_diameter=outer_diameter,
            inhibited_ends=inhibited_ends,
            map_dim=map_dim,
            density_ratio=density_ratio,
        )

    def get_maps(self) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
        """
        Return a tuple of two 2D arrays (map_x, map_y).
        Each is of shape (map_dim, map_dim), ranging from -1 to 1.
        """
        if self.maps is None:
            map_x, map_y = np.meshgrid(
                np.linspace(-1, 1, self.map_dim, dtype=np.float64),
                np.linspace(-1, 1, self.map_dim, dtype=np.float64),
            )
            self.maps = (map_x, map_y)
        return self.maps

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

    def get_contours(self, web_distance: float) -> list[NDArray[np.float64]]:
        """
        Return a list of contour arrays for the given web distance.
        Each contour is typically an (N,2) array of (row, col) points.
        """
        # get_contours is imported from machwave.core.math.geometric
        map_dist = self.normalize(web_distance)
        return get_contours(self.get_regression_map(), map_dist)

    def get_port_area(self, web_distance: float) -> float:
        """
        Return the grain's port area (open cross-sectional area) at the given web distance.
        Could be a scalar or array, depending on how the computations are done.
        """
        face_area = self.get_face_area(web_distance)
        return get_circle_area(self.outer_diameter) - face_area

    def get_face_area_interp_func(self) -> Callable[[float], float]:
        """
        Build and return an interpolation function that, given a normalized
        web distance, returns the face area in square meters.
        """
        if self.face_area_interp_func is None:
            regression_map = self.get_regression_map()
            valid = np.logical_not(self.get_mask())

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

            step_count = int(max_dist * self.map_dim) + 2
            distances = np.arange(step_count, dtype=np.float64) / self.map_dim

            n_le = np.searchsorted(values, distances, side="right")
            counts = float(values.size) - n_le.astype(np.float64)
            face_area_values = np.asarray(self.map_to_area(counts), dtype=np.float64)

            # Smooth + interpolate (adapt for small arrays)
            smoothed = face_area_values
            if face_area_values.size >= 7:
                window_length = min(31, int(face_area_values.size))
                if window_length % 2 == 0:
                    window_length -= 1
                polyorder = min(5, window_length - 2)
                if window_length >= 3 and polyorder >= 1:
                    smoothed = savgol_filter(face_area_values, window_length, polyorder)

            smoothed_arr = np.asarray(smoothed, dtype=np.float64)
            self.face_area_interp_func = interp1d(
                distances,
                smoothed_arr,
                bounds_error=False,
                fill_value=(float(smoothed_arr[0]), float(smoothed_arr[-1])),  # type: ignore[arg-type]
                assume_sorted=True,
            )

        return self.face_area_interp_func

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

    def get_burn_area_interp_func(self) -> Callable[[float], float]:
        """Return a cached interpolator for burn area [m^2] vs normalized web."""

        if self.burn_area_interp_func is None:
            regression_map = self.get_regression_map()
            valid = np.logical_not(self.get_mask())
            values = np.asarray(regression_map[valid], dtype=np.float64).ravel()
            if values.size == 0:
                self.burn_area_interp_func = 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_interp_func

            values.sort()
            max_dist = float(values[-1])
            step_count = int(max_dist * self.map_dim) + 2
            distances = np.arange(step_count, dtype=np.float64) / self.map_dim

            n_le = np.searchsorted(values, distances, side="right")
            counts = float(values.size) - n_le.astype(np.float64)
            face_area_values = np.asarray(self.map_to_area(counts), dtype=np.float64)

            perimeter_values = np.empty_like(distances, dtype=np.float64)
            for i, dist in enumerate(distances):
                contours = get_contours(regression_map, float(dist))
                perimeter_values[i] = float(
                    sum(
                        self.map_to_length(get_length(contour, self.map_dim))
                        for contour in contours
                    )
                )

            web_distances = np.asarray(self.denormalize(distances), dtype=np.float64)
            length_values = np.asarray(
                [self.get_length(float(wd)) for wd in web_distances], dtype=np.float64
            )
            core_area_values = perimeter_values * length_values
            total_face_area_values = (2 - self.inhibited_ends) * face_area_values
            burn_area_values = core_area_values + total_face_area_values

            smoothed = burn_area_values
            if burn_area_values.size >= 7:
                window_length = min(31, int(burn_area_values.size))
                if window_length % 2 == 0:
                    window_length -= 1
                polyorder = min(5, window_length - 2)
                if window_length >= 3 and polyorder >= 1:
                    smoothed = savgol_filter(burn_area_values, window_length, polyorder)

            smoothed_arr = np.asarray(smoothed, dtype=np.float64)
            self.burn_area_interp_func = interp1d(
                distances,
                smoothed_arr,
                bounds_error=False,
                fill_value=(float(smoothed_arr[0]), float(smoothed_arr[-1])),  # type: ignore[arg-type]
                assume_sorted=True,
            )

        return self.burn_area_interp_func

    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
        map_distance = self.normalize(web_distance)
        value = float(self.get_burn_area_interp_func()(map_distance))
        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)
        # Sum the lengths of all contour segments
        return float(
            sum(
                self.map_to_length(get_length(contour, self.map_dim))
                for contour in contours
            )
        )

    def get_core_area(self, web_distance: float) -> float:
        """
        Calculate the core (internal) area at the given web distance by
        multiplying perimeter by grain segment length (a 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 GrainGeometryError(
                "The web distance traveled is greater than the grain segment's web thickness."
            )

    def _get_active_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.
        """
        face_map = self.get_face_map(web_distance)
        mask = face_map == 1  # active material only
        y_indices, x_indices = np.where(mask)

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

        return y_indices, x_indices

    def _indices_to_normalized_coords(
        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_coords, y_coords) in normalized units.
        """
        center_shift = self.map_dim / 2
        x_coords = (x_indices - center_shift).astype(np.float64)
        y_coords = (y_indices - center_shift).astype(np.float64)
        return x_coords, y_coords

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

        Args:
            web_distance: Web distance traveled [m].

        Returns:
            Center of gravity [x, y, z] in meters from the port of the segment.
        """
        self._validate_web_distance(web_distance)
        y_indices, x_indices = self._get_active_material_indices(web_distance)
        x_coords, y_coords = self._indices_to_normalized_coords(y_indices, x_indices)

        # Convert to physical coordinates
        x_phys = np.asarray(self.map_to_length(x_coords), dtype=np.float64)
        y_phys = np.asarray(self.map_to_length(y_coords), dtype=np.float64)
        # For 2D grain, axial position is constant at segment center
        z_phys = np.full_like(x_phys, self.length / 2)

        return get_center_of_gravity(x_phys, y_phys, z_phys)

    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³].

        Returns:
            A 3x3 inertia tensor [kg-m^2].
        """
        self._validate_web_distance(web_distance)
        y_indices, x_indices = self._get_active_material_indices(web_distance)
        x_coords, y_coords = self._indices_to_normalized_coords(y_indices, x_indices)

        # Get the center of gravity to use as reference point
        cog = self.get_center_of_gravity(web_distance)
        current_length = self.get_length(web_distance)

        # Convert to meters
        x_phys = self.map_to_length(x_coords)
        y_phys = self.map_to_length(y_coords)

        # Shift to CoG frame
        x_rel = x_phys - cog[1]
        y_rel = y_phys - cog[2]
        z_rel = np.zeros_like(x_rel)  # 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
        n_elements = len(x_indices)
        element_mass = total_mass / n_elements

        # Get base inertia from point masses (radial only)
        moi = get_moment_of_inertia_tensor(x_rel, y_rel, z_rel, element_mass)

        # For 2D grain (uniform along axial direction), add axial contribution
        # For a uniform rod of length L with mass M: I_axial = M*L²/12 (about CoG)
        I_axial = total_mass * current_length**2 / 12

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

        return moi

get_burn_area(web_distance)

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

Source code in machwave/models/propulsion/grain/fmm/_2d.py 
202
203
204
205
206
207
208
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
    map_distance = self.normalize(web_distance)
    value = float(self.get_burn_area_interp_func()(map_distance))
    return max(0.0, value)

get_burn_area_interp_func()

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

Source code in machwave/models/propulsion/grain/fmm/_2d.py 
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
def get_burn_area_interp_func(self) -> Callable[[float], float]:
    """Return a cached interpolator for burn area [m^2] vs normalized web."""

    if self.burn_area_interp_func is None:
        regression_map = self.get_regression_map()
        valid = np.logical_not(self.get_mask())
        values = np.asarray(regression_map[valid], dtype=np.float64).ravel()
        if values.size == 0:
            self.burn_area_interp_func = 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_interp_func

        values.sort()
        max_dist = float(values[-1])
        step_count = int(max_dist * self.map_dim) + 2
        distances = np.arange(step_count, dtype=np.float64) / self.map_dim

        n_le = np.searchsorted(values, distances, side="right")
        counts = float(values.size) - n_le.astype(np.float64)
        face_area_values = np.asarray(self.map_to_area(counts), dtype=np.float64)

        perimeter_values = np.empty_like(distances, dtype=np.float64)
        for i, dist in enumerate(distances):
            contours = get_contours(regression_map, float(dist))
            perimeter_values[i] = float(
                sum(
                    self.map_to_length(get_length(contour, self.map_dim))
                    for contour in contours
                )
            )

        web_distances = np.asarray(self.denormalize(distances), dtype=np.float64)
        length_values = np.asarray(
            [self.get_length(float(wd)) for wd in web_distances], dtype=np.float64
        )
        core_area_values = perimeter_values * length_values
        total_face_area_values = (2 - self.inhibited_ends) * face_area_values
        burn_area_values = core_area_values + total_face_area_values

        smoothed = burn_area_values
        if burn_area_values.size >= 7:
            window_length = min(31, int(burn_area_values.size))
            if window_length % 2 == 0:
                window_length -= 1
            polyorder = min(5, window_length - 2)
            if window_length >= 3 and polyorder >= 1:
                smoothed = savgol_filter(burn_area_values, window_length, polyorder)

        smoothed_arr = np.asarray(smoothed, dtype=np.float64)
        self.burn_area_interp_func = interp1d(
            distances,
            smoothed_arr,
            bounds_error=False,
            fill_value=(float(smoothed_arr[0]), float(smoothed_arr[-1])),  # type: ignore[arg-type]
            assume_sorted=True,
        )

    return self.burn_area_interp_func

get_center_of_gravity(web_distance)

Calculates the center of gravity of a 2D FMM grain segment at a web distance.

Parameters:

Name Type Description Default
web_distance float

Web distance traveled [m].

 required

Returns:

Type Description
NDArray[float64]

Center of gravity [x, y, z] in meters from the port of the segment.
Source code in machwave/models/propulsion/grain/fmm/_2d.py 
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
def get_center_of_gravity(self, web_distance: float) -> NDArray[np.float64]:
    """
    Calculates the center of gravity of a 2D FMM grain segment at a web
    distance.

    Args:
        web_distance: Web distance traveled [m].

    Returns:
        Center of gravity [x, y, z] in meters from the port of the segment.
    """
    self._validate_web_distance(web_distance)
    y_indices, x_indices = self._get_active_material_indices(web_distance)
    x_coords, y_coords = self._indices_to_normalized_coords(y_indices, x_indices)

    # Convert to physical coordinates
    x_phys = np.asarray(self.map_to_length(x_coords), dtype=np.float64)
    y_phys = np.asarray(self.map_to_length(y_coords), dtype=np.float64)
    # For 2D grain, axial position is constant at segment center
    z_phys = np.full_like(x_phys, self.length / 2)

    return get_center_of_gravity(x_phys, y_phys, z_phys)

get_contours(web_distance)

Return a list of contour arrays for the given web distance. Each contour is typically an (N,2) array of (row, col) points.

Source code in machwave/models/propulsion/grain/fmm/_2d.py 
72
73
74
75
76
77
78
79
def get_contours(self, web_distance: float) -> list[NDArray[np.float64]]:
    """
    Return a list of contour arrays for the given web distance.
    Each contour is typically an (N,2) array of (row, col) points.
    """
    # get_contours is imported from machwave.core.math.geometric
    map_dist = self.normalize(web_distance)
    return get_contours(self.get_regression_map(), map_dist)

get_core_area(web_distance)

Calculate the core (internal) area at the given web distance by multiplying perimeter by grain segment length (a 2D approximation).

Source code in machwave/models/propulsion/grain/fmm/_2d.py 
223
224
225
226
227
228
def get_core_area(self, web_distance: float) -> float:
    """
    Calculate the core (internal) area at the given web distance by
    multiplying perimeter by grain segment length (a 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/propulsion/grain/fmm/_2d.py 
210
211
212
213
214
215
216
217
218
219
220
221
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)
    # Sum the lengths of all contour segments
    return float(
        sum(
            self.map_to_length(get_length(contour, self.map_dim))
            for contour in contours
        )
    )

get_face_area(web_distance)

Return the face area at the given web distance.

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

get_face_area_interp_func()

Build and return an interpolation function that, given a normalized web distance, returns the face area in square meters.

Source code in machwave/models/propulsion/grain/fmm/_2d.py 
 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
def get_face_area_interp_func(self) -> Callable[[float], float]:
    """
    Build and return an interpolation function that, given a normalized
    web distance, returns the face area in square meters.
    """
    if self.face_area_interp_func is None:
        regression_map = self.get_regression_map()
        valid = np.logical_not(self.get_mask())

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

        step_count = int(max_dist * self.map_dim) + 2
        distances = np.arange(step_count, dtype=np.float64) / self.map_dim

        n_le = np.searchsorted(values, distances, side="right")
        counts = float(values.size) - n_le.astype(np.float64)
        face_area_values = np.asarray(self.map_to_area(counts), dtype=np.float64)

        # Smooth + interpolate (adapt for small arrays)
        smoothed = face_area_values
        if face_area_values.size >= 7:
            window_length = min(31, int(face_area_values.size))
            if window_length % 2 == 0:
                window_length -= 1
            polyorder = min(5, window_length - 2)
            if window_length >= 3 and polyorder >= 1:
                smoothed = savgol_filter(face_area_values, window_length, polyorder)

        smoothed_arr = np.asarray(smoothed, dtype=np.float64)
        self.face_area_interp_func = interp1d(
            distances,
            smoothed_arr,
            bounds_error=False,
            fill_value=(float(smoothed_arr[0]), float(smoothed_arr[-1])),  # type: ignore[arg-type]
            assume_sorted=True,
        )

    return self.face_area_interp_func

get_maps()

Return a tuple of two 2D arrays (map_x, map_y). Each is of shape (map_dim, map_dim), ranging from -1 to 1.

Source code in machwave/models/propulsion/grain/fmm/_2d.py 
50
51
52
53
54
55
56
57
58
59
60
61
def get_maps(self) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
    """
    Return a tuple of two 2D arrays (map_x, map_y).
    Each is of shape (map_dim, map_dim), ranging from -1 to 1.
    """
    if self.maps is None:
        map_x, map_y = np.meshgrid(
            np.linspace(-1, 1, self.map_dim, dtype=np.float64),
            np.linspace(-1, 1, self.map_dim, dtype=np.float64),
        )
        self.maps = (map_x, map_y)
    return self.maps

get_mask()

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

Source code in machwave/models/propulsion/grain/fmm/_2d.py 
63
64
65
66
67
68
69
70
def get_mask(self) -> NDArray[np.bool_]:
    """
    Return a boolean mask indicating which points lie outside the unit circle.
    """
    if self.mask is None:
        map_x, map_y = self.get_maps()
        self.mask = (map_x**2 + map_y**2) > 1
    return self.mask

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³].

 required

Returns:

Type Description
NDArray[float64]

A 3x3 inertia tensor [kg-m^2].
Source code in machwave/models/propulsion/grain/fmm/_2d.py 
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
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³].

    Returns:
        A 3x3 inertia tensor [kg-m^2].
    """
    self._validate_web_distance(web_distance)
    y_indices, x_indices = self._get_active_material_indices(web_distance)
    x_coords, y_coords = self._indices_to_normalized_coords(y_indices, x_indices)

    # Get the center of gravity to use as reference point
    cog = self.get_center_of_gravity(web_distance)
    current_length = self.get_length(web_distance)

    # Convert to meters
    x_phys = self.map_to_length(x_coords)
    y_phys = self.map_to_length(y_coords)

    # Shift to CoG frame
    x_rel = x_phys - cog[1]
    y_rel = y_phys - cog[2]
    z_rel = np.zeros_like(x_rel)  # 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
    n_elements = len(x_indices)
    element_mass = total_mass / n_elements

    # Get base inertia from point masses (radial only)
    moi = get_moment_of_inertia_tensor(x_rel, y_rel, z_rel, element_mass)

    # For 2D grain (uniform along axial direction), add axial contribution
    # For a uniform rod of length L with mass M: I_axial = M*L²/12 (about CoG)
    I_axial = total_mass * current_length**2 / 12

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

    return moi

get_port_area(web_distance)

Return the grain's port area (open cross-sectional area) at the given web distance. Could be a scalar or array, depending on how the computations are done.

Source code in machwave/models/propulsion/grain/fmm/_2d.py 
81
82
83
84
85
86
87
def get_port_area(self, web_distance: float) -> float:
    """
    Return the grain's port area (open cross-sectional area) at the given web distance.
    Could be a scalar or array, depending on how the computations are done.
    """
    face_area = self.get_face_area(web_distance)
    return get_circle_area(self.outer_diameter) - face_area

FMMGrainSegment3D

Bases: FMMGrainSegment, GrainSegment3D, ABC

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

This class was inspired by the Andrew Reilley's software openMotor, in particular the fmm module. openMotor's repository can be accessed at: https://github.com/reilleya/openMotor

Source code in machwave/models/propulsion/grain/fmm/_3d.py 
 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
class FMMGrainSegment3D(FMMGrainSegment, GrainSegment3D, ABC):
    """
    Fast Marching Method (FMM) implementation for 3D grain segment.

    This class was inspired by the Andrew Reilley's software openMotor, in
    particular the fmm module.
    openMotor's repository can be accessed at:
    https://github.com/reilleya/openMotor
    """

    def __init__(
        self,
        length: float,
        outer_diameter: float,
        inhibited_ends: int = 0,
        map_dim: int = 100,
        density_ratio: float = 1.0,
    ) -> None:
        self.burn_area_interp_func: Callable[[float], float] | None = None
        super().__init__(
            length=length,
            outer_diameter=outer_diameter,
            inhibited_ends=inhibited_ends,
            map_dim=map_dim,
            density_ratio=density_ratio,
        )

    def get_port_area(self, web_distance: float, z: float) -> 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 (in meters) along the grain, where z=0 is the top
                and z=self.length is the bottom (or vice versa, depending on
                geometry setup).

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

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

        face_area = float(self.map_to_area(float(np.count_nonzero(solid[z_index]))))
        return get_circle_area(self.outer_diameter) - face_area

    def get_normalized_length(self) -> int:
        return int(self.map_dim * self.length / self.outer_diameter)

    def get_maps(
        self,
    ) -> tuple[
        NDArray[np.float64],
        NDArray[np.float64],
        NDArray[np.float64],
    ]:
        if self.maps is None:
            map_y, map_z, map_x = np.meshgrid(
                np.linspace(-1, 1, self.map_dim),
                np.linspace(1, 0, self.get_normalized_length()),  # z axis
                np.linspace(-1, 1, self.map_dim),
            )

            self.maps = (map_x, map_y, map_z)

        return self.maps

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

        return self.mask

    def get_contours(
        self, web_distance: float, length_normalized: float
    ) -> list[NDArray[np.float64]]:
        map_dist = self.normalize(web_distance)
        valid = np.logical_not(self.get_mask())
        boolean_3d = np.logical_and(self.get_regression_map() > map_dist, valid)

        z_index = int(round(length_normalized))
        boolean_slice_2d = boolean_3d[z_index]

        return get_contours(
            boolean_slice_2d,
            map_dist,
        )

    def _get_burn_area_uncached(self, *, map_dist: float) -> float:
        valid = np.logical_not(self.get_mask())
        boolean_3d = np.logical_and(self.get_regression_map() > map_dist, valid)

        web_distance = float(self.denormalize(map_dist))
        length_factor = self.get_length(web_distance=web_distance) / self.map_dim

        total = 0.0
        for z_index in range(self.get_normalized_length()):
            boolean_slice_2d = boolean_3d[z_index]
            contours = get_contours(boolean_slice_2d, map_dist)
            perimeter = sum(
                self.map_to_length(get_length(contour, self.map_dim))
                for contour in contours
            )
            total += float(perimeter) * float(length_factor)

        return float(total)

    def get_burn_area_interp_func(self) -> Callable[[float], float]:
        """Return a cached interpolator for burn area [m^2] vs normalized web."""

        if self.burn_area_interp_func is None:
            regression_map = self.get_regression_map()
            valid = np.logical_not(self.get_mask())

            values = np.asarray(regression_map[valid], dtype=np.float64).ravel()
            if values.size == 0:
                self.burn_area_interp_func = 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_interp_func

            values.sort()
            max_dist = float(values[-1])
            oversample = 3
            denom = float(self.map_dim * oversample)
            step_count = int(max_dist * denom) + 2
            distances = np.arange(step_count, dtype=np.float64) / denom

            burn_area_values = np.empty_like(distances, dtype=np.float64)
            for i, dist in enumerate(distances):
                burn_area_values[i] = self._get_burn_area_uncached(map_dist=float(dist))

            burn_area_values = np.asarray(burn_area_values, dtype=np.float64)
            self.burn_area_interp_func = interp1d(
                distances,
                burn_area_values,
                bounds_error=False,
                fill_value=(float(burn_area_values[0]), float(burn_area_values[-1])),  # type: ignore[arg-type]
                assume_sorted=True,
            )

        return self.burn_area_interp_func

    def get_burn_area(self, web_distance: float) -> float:
        if web_distance > self.get_web_thickness():
            return 0.0
        map_distance = self.normalize(web_distance)
        value = float(self.get_burn_area_interp_func()(map_distance))
        return max(0.0, value)

    def get_volume_per_element(self) -> float:
        return (float(self.denormalize(self.get_cell_size())) * 2) ** 3

    def get_volume(self, web_distance: float) -> float:
        face_map = self.get_face_map(web_distance=web_distance)
        active_elements = np.count_nonzero(face_map == 1)
        volume_per_element = self.get_volume_per_element()
        return active_elements * volume_per_element

    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 GrainGeometryError(
                "The web distance traveled is greater than the grain "
                "segment's web thickness."
            )

    def _get_active_material_indices(
        self, web_distance: float
    ) -> tuple[NDArray[np.int_], 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 (z_indices, y_indices, x_indices) for active material.

        Raises:
            GrainGeometryError: If no active material is found.
        """
        face_map = self.get_face_map(web_distance)
        mask = face_map == 1  # active material only
        z_indices, y_indices, x_indices = np.where(mask)

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

        return z_indices, y_indices, x_indices

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

        Args:
            z_indices: Z-axis (axial) indices from face map.
            y_indices: Y-axis indices from face map.
            x_indices: X-axis indices from face map.

        Returns:
            Tuple of (x_coords, y_coords, z_coords) in normalized units.
        """
        center_shift = self.map_dim / 2
        x_coords = (x_indices - center_shift).astype(np.float64)
        y_coords = (y_indices - center_shift).astype(np.float64)
        z_coords = z_indices.astype(np.float64)
        return x_coords, y_coords, z_coords

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

        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)
        z_indices, y_indices, x_indices = self._get_active_material_indices(
            web_distance
        )
        x_coords, y_coords, z_coords = self._indices_to_normalized_coords(
            z_indices, y_indices, x_indices
        )

        # Convert normalized coordinates into physical meters
        x_phys = np.asarray(self.map_to_length(x_coords), dtype=np.float64)
        y_phys = np.asarray(self.map_to_length(y_coords), dtype=np.float64)
        z_phys = np.asarray(self.map_to_length(z_coords), dtype=np.float64)

        return get_center_of_gravity(x_phys, y_phys, z_phys)

    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)
        z_indices, y_indices, x_indices = self._get_active_material_indices(
            web_distance
        )
        x_coords, y_coords, z_coords = self._indices_to_normalized_coords(
            z_indices, y_indices, x_indices
        )

        # Get the center of gravity
        cog = self.get_center_of_gravity(web_distance)

        # Convert to meters
        x_phys = self.map_to_length(x_coords)
        y_phys = self.map_to_length(y_coords)
        z_phys = self.map_to_length(z_coords)

        # Shift to CoG frame
        x_rel = x_phys - cog[1]
        y_rel = y_phys - cog[2]
        z_rel = z_phys - cog[0]

        # Calculate element mass
        element_volume = self.get_volume_per_element()
        element_mass = element_volume * ideal_density * self.density_ratio

        # Use core function to compute inertia tensor
        return get_moment_of_inertia_tensor(x_rel, y_rel, z_rel, element_mass)

get_burn_area_interp_func()

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

Source code in machwave/models/propulsion/grain/fmm/_3d.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
def get_burn_area_interp_func(self) -> Callable[[float], float]:
    """Return a cached interpolator for burn area [m^2] vs normalized web."""

    if self.burn_area_interp_func is None:
        regression_map = self.get_regression_map()
        valid = np.logical_not(self.get_mask())

        values = np.asarray(regression_map[valid], dtype=np.float64).ravel()
        if values.size == 0:
            self.burn_area_interp_func = 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_interp_func

        values.sort()
        max_dist = float(values[-1])
        oversample = 3
        denom = float(self.map_dim * oversample)
        step_count = int(max_dist * denom) + 2
        distances = np.arange(step_count, dtype=np.float64) / denom

        burn_area_values = np.empty_like(distances, dtype=np.float64)
        for i, dist in enumerate(distances):
            burn_area_values[i] = self._get_burn_area_uncached(map_dist=float(dist))

        burn_area_values = np.asarray(burn_area_values, dtype=np.float64)
        self.burn_area_interp_func = interp1d(
            distances,
            burn_area_values,
            bounds_error=False,
            fill_value=(float(burn_area_values[0]), float(burn_area_values[-1])),  # type: ignore[arg-type]
            assume_sorted=True,
        )

    return self.burn_area_interp_func

get_center_of_gravity(web_distance)

Calculates the center of gravity of a 3D FMM grain segment at a web distance.

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/propulsion/grain/fmm/_3d.py 
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
def get_center_of_gravity(self, web_distance: float) -> NDArray[np.float64]:
    """
    Calculates the center of gravity of a 3D FMM grain segment at a web
    distance.

    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)
    z_indices, y_indices, x_indices = self._get_active_material_indices(
        web_distance
    )
    x_coords, y_coords, z_coords = self._indices_to_normalized_coords(
        z_indices, y_indices, x_indices
    )

    # Convert normalized coordinates into physical meters
    x_phys = np.asarray(self.map_to_length(x_coords), dtype=np.float64)
    y_phys = np.asarray(self.map_to_length(y_coords), dtype=np.float64)
    z_phys = np.asarray(self.map_to_length(z_coords), dtype=np.float64)

    return get_center_of_gravity(x_phys, y_phys, z_phys)

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/propulsion/grain/fmm/_3d.py 
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
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)
    z_indices, y_indices, x_indices = self._get_active_material_indices(
        web_distance
    )
    x_coords, y_coords, z_coords = self._indices_to_normalized_coords(
        z_indices, y_indices, x_indices
    )

    # Get the center of gravity
    cog = self.get_center_of_gravity(web_distance)

    # Convert to meters
    x_phys = self.map_to_length(x_coords)
    y_phys = self.map_to_length(y_coords)
    z_phys = self.map_to_length(z_coords)

    # Shift to CoG frame
    x_rel = x_phys - cog[1]
    y_rel = y_phys - cog[2]
    z_rel = z_phys - cog[0]

    # Calculate element mass
    element_volume = self.get_volume_per_element()
    element_mass = element_volume * ideal_density * self.density_ratio

    # Use core function to compute inertia tensor
    return get_moment_of_inertia_tensor(x_rel, y_rel, z_rel, element_mass)

get_port_area(web_distance, z)

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 (in meters) along the grain, where z=0 is the top and z=self.length is the bottom (or vice versa, depending on geometry setup).

 required

Returns:

Type Description
float

A float representing the port area at the specified z slice, in m².
Source code in machwave/models/propulsion/grain/fmm/_3d.py 
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
def get_port_area(self, web_distance: float, z: float) -> 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 (in meters) along the grain, where z=0 is the top
            and z=self.length is the bottom (or vice versa, depending on
            geometry setup).

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

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

    face_area = float(self.map_to_area(float(np.count_nonzero(solid[z_index]))))
    return get_circle_area(self.outer_diameter) - face_area

FMMSTLGrainSegment

Bases: FMMGrainSegment3D, ABC

Fast Marching Method (FMM) implementation for a grain segment obtained from an STL file.

Source code in machwave/models/propulsion/grain/fmm/stl.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
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
class FMMSTLGrainSegment(FMMGrainSegment3D, ABC):
    """
    Fast Marching Method (FMM) implementation for a grain segment obtained
    from an STL file.
    """

    def __init__(
        self,
        file_path: str,
        outer_diameter: float,
        length: float,
        inhibited_ends: int = 0,
        map_dim: int = 50,
    ) -> None:
        self.file_path = file_path
        self.outer_diameter = outer_diameter
        self.length = length

        # "Cache" variables:
        self.face_area_interp_func = None

        super().__init__(
            length=length,
            outer_diameter=outer_diameter,
            inhibited_ends=inhibited_ends,
            map_dim=map_dim,
        )

    def validate(self) -> None:
        if not self.map_dim >= 20:
            raise GrainGeometryError(
                f"Map dimension must be at least 20 for STL grains, got {self.map_dim}"
            )

    def get_voxel_size(self) -> float:
        """
        NOTE: Only returns correct voxel size if map_dim is an odd number.

        :return: the voxel edge size.
        :rtype: float
        """
        return self.outer_diameter / int(self.map_dim - 1)

    def get_initial_face_map(self) -> np.typing.NDArray[np.int_]:
        """
        Generate a map by voxelizing an STL file. Uses trimesh library.

        NOTE: Still needs to convert boolean matrix to masked array.
        """
        mesh: trimesh.Trimesh = trimesh.load_mesh(self.file_path)
        assert mesh.is_watertight, "Mesh must be watertight"

        volume = mesh.voxelized(pitch=self.get_voxel_size()).fill()
        voxel_map: np.typing.NDArray[np.int_] = (
            volume.matrix.view(np.ndarray).transpose().astype(np.int_)
        )

        assert voxel_map.shape == self.get_maps()[0].shape, (
            "Generated map shape mismatch"
        )

        return voxel_map

get_initial_face_map()

Generate a map by voxelizing an STL file. Uses trimesh library.

NOTE: Still needs to convert boolean matrix to masked array.

Source code in machwave/models/propulsion/grain/fmm/stl.py 
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def get_initial_face_map(self) -> np.typing.NDArray[np.int_]:
    """
    Generate a map by voxelizing an STL file. Uses trimesh library.

    NOTE: Still needs to convert boolean matrix to masked array.
    """
    mesh: trimesh.Trimesh = trimesh.load_mesh(self.file_path)
    assert mesh.is_watertight, "Mesh must be watertight"

    volume = mesh.voxelized(pitch=self.get_voxel_size()).fill()
    voxel_map: np.typing.NDArray[np.int_] = (
        volume.matrix.view(np.ndarray).transpose().astype(np.int_)
    )

    assert voxel_map.shape == self.get_maps()[0].shape, (
        "Generated map shape mismatch"
    )

    return voxel_map

get_voxel_size()

NOTE: Only returns correct voxel size if map_dim is an odd number.

:return: the voxel edge size. :rtype: float

Source code in machwave/models/propulsion/grain/fmm/stl.py 
44
45
46
47
48
49
50
51
def get_voxel_size(self) -> float:
    """
    NOTE: Only returns correct voxel size if map_dim is an odd number.

    :return: the voxel edge size.
    :rtype: float
    """
    return self.outer_diameter / int(self.map_dim - 1)