oxiphysics_geometry/level_set/
functions.rs

1//! Auto-generated module
2//!
3//! 🤖 Generated with [SplitRS](https://github.com/cool-japan/splitrs)
4
5use super::types::LevelSetField;
6
7/// Godunov upwind gradient magnitude for reinitialization.
8///
9/// Computes the Godunov upwind approximation of |∇φ| needed for the
10/// Sussman reinitialization scheme.
11pub(super) fn godunov_grad_mag(
12    field: &LevelSetField,
13    ix: usize,
14    iy: usize,
15    iz: usize,
16    sign: f64,
17) -> f64 {
18    let dx = field.dx;
19    let phi_c = field.get(ix, iy, iz);
20    let d_xm = if ix > 0 {
21        (phi_c - field.get(ix - 1, iy, iz)) / dx
22    } else {
23        0.0
24    };
25    let d_xp = if ix + 1 < field.nx {
26        (field.get(ix + 1, iy, iz) - phi_c) / dx
27    } else {
28        0.0
29    };
30    let d_ym = if iy > 0 {
31        (phi_c - field.get(ix, iy - 1, iz)) / dx
32    } else {
33        0.0
34    };
35    let d_yp = if iy + 1 < field.ny {
36        (field.get(ix, iy + 1, iz) - phi_c) / dx
37    } else {
38        0.0
39    };
40    let d_zm = if iz > 0 {
41        (phi_c - field.get(ix, iy, iz - 1)) / dx
42    } else {
43        0.0
44    };
45    let d_zp = if iz + 1 < field.nz {
46        (field.get(ix, iy, iz + 1) - phi_c) / dx
47    } else {
48        0.0
49    };
50    let gx = if sign > 0.0 {
51        d_xm.max(0.0).powi(2).max((-d_xp).max(0.0).powi(2))
52    } else {
53        (-d_xm).max(0.0).powi(2).max(d_xp.max(0.0).powi(2))
54    };
55    let gy = if sign > 0.0 {
56        d_ym.max(0.0).powi(2).max((-d_yp).max(0.0).powi(2))
57    } else {
58        (-d_ym).max(0.0).powi(2).max(d_yp.max(0.0).powi(2))
59    };
60    let gz = if sign > 0.0 {
61        d_zm.max(0.0).powi(2).max((-d_zp).max(0.0).powi(2))
62    } else {
63        (-d_zm).max(0.0).powi(2).max(d_zp.max(0.0).powi(2))
64    };
65    (gx + gy + gz).sqrt()
66}
67/// Linear interpolation of vertex position on a cube edge.
68pub(super) fn interp_vertex(pa: [f64; 3], pb: [f64; 3], va: f64, vb: f64, iso: f64) -> [f64; 3] {
69    let dv = vb - va;
70    let t = if dv.abs() < 1e-30 {
71        0.5
72    } else {
73        (iso - va) / dv
74    };
75    [
76        pa[0] + t * (pb[0] - pa[0]),
77        pa[1] + t * (pb[1] - pa[1]),
78        pa[2] + t * (pb[2] - pa[2]),
79    ]
80}
81/// Edges of a unit cube: (corner_a, corner_b) pairs for the 12 edges.
82pub(super) const MC_EDGE_CORNERS: [(usize, usize); 12] = [
83    (0, 1),
84    (1, 2),
85    (2, 3),
86    (3, 0),
87    (4, 5),
88    (5, 6),
89    (6, 7),
90    (7, 4),
91    (0, 4),
92    (1, 5),
93    (2, 6),
94    (3, 7),
95];
96pub(super) const MC_EDGE_TABLE: [u16; 256] = [
97    0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a,
98    0xd03, 0xe09, 0xf00, 0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895,
99    0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, 0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435,
100    0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 0x3a0, 0x2a9, 0x1a3, 0x0aa,
101    0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, 0x460,
102    0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963,
103    0xa69, 0xb60, 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff,
104    0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c,
105    0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6,
106    0x2cf, 0x1c5, 0x0cc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, 0x8c0, 0x9c9,
107    0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0x0cc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9,
108    0x7c0, 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x055, 0x35f, 0x256,
109    0x55a, 0x453, 0x759, 0x650, 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc,
110    0x3f5, 0x0ff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f,
111    0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x066, 0x76a, 0x663, 0x569, 0x460, 0xca0, 0xda9, 0xea3,
112    0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0x0aa, 0x1a3, 0x2a9, 0x3a0,
113    0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x835, 0xb3f, 0xa36, 0x53c, 0x435, 0x73f, 0x636, 0x13a,
114    0x033, 0x339, 0x230, 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795,
115    0x49f, 0x596, 0x29a, 0x393, 0x099, 0x190, 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905,
116    0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x000,
117];
118/// Marching cubes triangle table (Lorensen & Cline, 1987).
119///
120/// For each of the 256 cube configurations, lists the edge indices of the
121/// triangles to generate.  Each row ends with -1 as sentinel.
122/// (Abbreviated to key configurations; full 256-entry tables follow the same pattern.)
123pub(super) const MC_TRI_TABLE: [[i8; 16]; 256] = [
124    [
125        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
126    ],
127    [0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
128    [0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
129    [1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
130    [1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
131    [0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
132    [9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
133    [2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1],
134    [3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
135    [0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
136    [1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
137    [1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1],
138    [3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
139    [0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1],
140    [3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1],
141    [9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
142    [4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
143    [4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
144    [0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
145    [4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1],
146    [1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
147    [3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1],
148    [9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1],
149    [2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1],
150    [8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
151    [11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1],
152    [9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1],
153    [4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1],
154    [3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1],
155    [1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1],
156    [4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1],
157    [4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1],
158    [
159        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
160    ],
161    [
162        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
163    ],
164    [
165        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
166    ],
167    [
168        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
169    ],
170    [
171        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
172    ],
173    [
174        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
175    ],
176    [
177        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
178    ],
179    [
180        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
181    ],
182    [
183        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
184    ],
185    [
186        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
187    ],
188    [
189        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
190    ],
191    [
192        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
193    ],
194    [
195        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
196    ],
197    [
198        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
199    ],
200    [
201        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
202    ],
203    [
204        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
205    ],
206    [
207        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
208    ],
209    [
210        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
211    ],
212    [
213        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
214    ],
215    [
216        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
217    ],
218    [
219        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
220    ],
221    [
222        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
223    ],
224    [
225        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
226    ],
227    [
228        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
229    ],
230    [
231        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
232    ],
233    [
234        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
235    ],
236    [
237        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
238    ],
239    [
240        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
241    ],
242    [
243        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
244    ],
245    [
246        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
247    ],
248    [
249        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
250    ],
251    [
252        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
253    ],
254    [
255        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
256    ],
257    [
258        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
259    ],
260    [
261        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
262    ],
263    [
264        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
265    ],
266    [
267        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
268    ],
269    [
270        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
271    ],
272    [
273        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
274    ],
275    [
276        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
277    ],
278    [
279        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
280    ],
281    [
282        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
283    ],
284    [
285        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
286    ],
287    [
288        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
289    ],
290    [
291        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
292    ],
293    [
294        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
295    ],
296    [
297        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
298    ],
299    [
300        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
301    ],
302    [
303        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
304    ],
305    [
306        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
307    ],
308    [
309        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
310    ],
311    [
312        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
313    ],
314    [
315        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
316    ],
317    [
318        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
319    ],
320    [
321        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
322    ],
323    [
324        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
325    ],
326    [
327        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
328    ],
329    [
330        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
331    ],
332    [
333        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
334    ],
335    [
336        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
337    ],
338    [
339        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
340    ],
341    [
342        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
343    ],
344    [
345        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
346    ],
347    [
348        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
349    ],
350    [
351        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
352    ],
353    [
354        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
355    ],
356    [
357        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
358    ],
359    [
360        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
361    ],
362    [
363        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
364    ],
365    [
366        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
367    ],
368    [
369        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
370    ],
371    [
372        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
373    ],
374    [
375        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
376    ],
377    [
378        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
379    ],
380    [
381        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
382    ],
383    [
384        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
385    ],
386    [
387        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
388    ],
389    [
390        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
391    ],
392    [
393        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
394    ],
395    [
396        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
397    ],
398    [
399        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
400    ],
401    [
402        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
403    ],
404    [
405        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
406    ],
407    [
408        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
409    ],
410    [
411        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
412    ],
413    [
414        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
415    ],
416    [
417        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
418    ],
419    [
420        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
421    ],
422    [
423        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
424    ],
425    [
426        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
427    ],
428    [
429        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
430    ],
431    [
432        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
433    ],
434    [
435        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
436    ],
437    [
438        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
439    ],
440    [
441        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
442    ],
443    [
444        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
445    ],
446    [
447        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
448    ],
449    [
450        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
451    ],
452    [
453        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
454    ],
455    [
456        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
457    ],
458    [
459        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
460    ],
461    [
462        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
463    ],
464    [
465        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
466    ],
467    [
468        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
469    ],
470    [
471        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
472    ],
473    [
474        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
475    ],
476    [
477        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
478    ],
479    [
480        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
481    ],
482    [
483        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
484    ],
485    [
486        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
487    ],
488    [
489        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
490    ],
491    [
492        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
493    ],
494    [
495        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
496    ],
497    [
498        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
499    ],
500    [
501        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
502    ],
503    [
504        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
505    ],
506    [
507        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
508    ],
509    [
510        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
511    ],
512    [
513        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
514    ],
515    [
516        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
517    ],
518    [
519        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
520    ],
521    [
522        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
523    ],
524    [
525        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
526    ],
527    [
528        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
529    ],
530    [
531        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
532    ],
533    [
534        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
535    ],
536    [
537        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
538    ],
539    [
540        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
541    ],
542    [
543        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
544    ],
545    [
546        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
547    ],
548    [
549        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
550    ],
551    [
552        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
553    ],
554    [
555        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
556    ],
557    [
558        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
559    ],
560    [
561        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
562    ],
563    [
564        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
565    ],
566    [
567        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
568    ],
569    [
570        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
571    ],
572    [
573        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
574    ],
575    [
576        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
577    ],
578    [
579        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
580    ],
581    [
582        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
583    ],
584    [
585        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
586    ],
587    [
588        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
589    ],
590    [
591        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
592    ],
593    [
594        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
595    ],
596    [
597        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
598    ],
599    [
600        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
601    ],
602    [
603        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
604    ],
605    [
606        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
607    ],
608    [
609        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
610    ],
611    [
612        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
613    ],
614    [
615        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
616    ],
617    [
618        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
619    ],
620    [
621        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
622    ],
623    [
624        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
625    ],
626    [
627        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
628    ],
629    [
630        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
631    ],
632    [
633        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
634    ],
635    [
636        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
637    ],
638    [
639        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
640    ],
641    [
642        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
643    ],
644    [
645        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
646    ],
647    [
648        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
649    ],
650    [
651        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
652    ],
653    [
654        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
655    ],
656    [
657        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
658    ],
659    [
660        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
661    ],
662    [
663        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
664    ],
665    [
666        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
667    ],
668    [
669        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
670    ],
671    [
672        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
673    ],
674    [
675        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
676    ],
677    [
678        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
679    ],
680    [
681        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
682    ],
683    [
684        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
685    ],
686    [
687        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
688    ],
689    [
690        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
691    ],
692    [
693        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
694    ],
695    [
696        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
697    ],
698    [
699        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
700    ],
701    [
702        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
703    ],
704    [
705        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
706    ],
707    [
708        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
709    ],
710    [
711        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
712    ],
713    [
714        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
715    ],
716    [
717        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
718    ],
719    [
720        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
721    ],
722    [
723        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
724    ],
725    [
726        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
727    ],
728    [
729        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
730    ],
731    [
732        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
733    ],
734    [
735        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
736    ],
737    [
738        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
739    ],
740    [
741        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
742    ],
743    [
744        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
745    ],
746    [
747        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
748    ],
749    [
750        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
751    ],
752    [
753        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
754    ],
755    [
756        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
757    ],
758    [
759        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
760    ],
761    [
762        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
763    ],
764    [
765        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
766    ],
767    [
768        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
769    ],
770    [
771        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
772    ],
773    [
774        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
775    ],
776    [
777        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
778    ],
779    [
780        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
781    ],
782    [
783        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
784    ],
785    [
786        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
787    ],
788    [
789        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
790    ],
791    [
792        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
793    ],
794    [
795        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
796    ],
797    [
798        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
799    ],
800    [
801        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
802    ],
803    [
804        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
805    ],
806    [
807        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
808    ],
809    [
810        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
811    ],
812    [
813        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
814    ],
815    [
816        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
817    ],
818    [
819        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
820    ],
821    [
822        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
823    ],
824    [
825        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
826    ],
827    [
828        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
829    ],
830];
831#[cfg(test)]
832mod tests {
833    use super::*;
834    use crate::level_set::IsosurfaceExtraction;
835    use crate::level_set::LevelSetBoolean;
836    use crate::level_set::LevelSetEvolution;
837    use crate::level_set::LevelSetReinit;
838    use crate::level_set::SignedDistanceTransform;
839    fn make_sphere_field(radius: f64) -> LevelSetField {
840        let mut f = LevelSetField::new(8, 8, 8, 0.25, [-1.0, -1.0, -1.0]);
841        f.init_sphere([0.0, 0.0, 0.0], radius);
842        f
843    }
844    #[test]
845    fn test_field_new_dimensions() {
846        let f = LevelSetField::new(4, 5, 6, 0.1, [0.0, 0.0, 0.0]);
847        assert_eq!(f.nx, 4);
848        assert_eq!(f.ny, 5);
849        assert_eq!(f.nz, 6);
850        assert_eq!(f.data.len(), 4 * 5 * 6);
851    }
852    #[test]
853    fn test_field_set_get() {
854        let mut f = LevelSetField::new(4, 4, 4, 0.1, [0.0, 0.0, 0.0]);
855        f.set(1, 2, 3, -0.5);
856        assert!((f.get(1, 2, 3) + 0.5).abs() < 1e-12);
857    }
858    #[test]
859    fn test_field_idx_linear() {
860        let f = LevelSetField::new(4, 4, 4, 0.1, [0.0, 0.0, 0.0]);
861        let idx = f.idx(1, 2, 3);
862        assert_eq!(idx, 57);
863    }
864    #[test]
865    fn test_field_node_pos() {
866        let f = LevelSetField::new(5, 5, 5, 0.5, [0.0, 0.0, 0.0]);
867        let p = f.node_pos(2, 0, 0);
868        assert!((p[0] - 1.0).abs() < 1e-12);
869        assert!(p[1].abs() < 1e-12);
870    }
871    #[test]
872    fn test_field_out_of_bounds_get() {
873        let f = LevelSetField::new(4, 4, 4, 0.1, [0.0, 0.0, 0.0]);
874        assert_eq!(f.get(10, 0, 0), f64::MAX);
875    }
876    #[test]
877    fn test_field_out_of_bounds_set_noop() {
878        let mut f = LevelSetField::new(4, 4, 4, 0.1, [0.0, 0.0, 0.0]);
879        f.set(100, 0, 0, 42.0);
880        assert_eq!(f.data.iter().filter(|&&v| v == 42.0).count(), 0);
881    }
882    #[test]
883    fn test_field_init_sphere_center_negative() {
884        let f = make_sphere_field(0.5);
885        let center_val = f.interpolate([0.0, 0.0, 0.0]);
886        assert!(
887            center_val < 0.0,
888            "center should be inside sphere: {center_val}"
889        );
890    }
891    #[test]
892    fn test_field_init_sphere_far_point_positive() {
893        let f = make_sphere_field(0.5);
894        let far_val = f.get(0, 0, 0);
895        assert!(far_val > 0.0, "corner should be outside sphere: {far_val}");
896    }
897    #[test]
898    fn test_field_init_box_inside_negative() {
899        let mut f = LevelSetField::new(8, 8, 8, 0.25, [-1.0, -1.0, -1.0]);
900        f.init_box([-0.3, -0.3, -0.3], [0.3, 0.3, 0.3]);
901        let val = f.interpolate([0.0, 0.0, 0.0]);
902        assert!(val < 0.0, "inside box should be negative: {val}");
903    }
904    #[test]
905    fn test_field_count_inside_sphere() {
906        let f = make_sphere_field(0.5);
907        let n_inside = f.count_inside();
908        assert!(n_inside > 0, "some nodes should be inside sphere");
909        assert!(n_inside < f.data.len(), "not all nodes inside sphere");
910    }
911    #[test]
912    fn test_field_min_max_sphere() {
913        let f = make_sphere_field(0.5);
914        let vmin = f.min_value();
915        let vmax = f.max_value();
916        assert!(vmin < 0.0, "min should be negative: {vmin}");
917        assert!(vmax > 0.0, "max should be positive: {vmax}");
918    }
919    #[test]
920    fn test_field_gradient_sphere_on_surface() {
921        let f = make_sphere_field(0.5);
922        let mid = 4usize;
923        let g = f.gradient(mid, mid, mid);
924        let mag = (g[0] * g[0] + g[1] * g[1] + g[2] * g[2]).sqrt();
925        assert!(mag.is_finite(), "gradient magnitude should be finite");
926    }
927    #[test]
928    fn test_field_interpolate_on_node() {
929        let mut f = LevelSetField::new(5, 5, 5, 0.5, [0.0, 0.0, 0.0]);
930        f.set(1, 1, 1, 3.125);
931        let p = f.node_pos(1, 1, 1);
932        let val = f.interpolate(p);
933        assert!(val.is_finite());
934    }
935    #[test]
936    fn test_reinit_sphere_preserves_zero_set() {
937        let mut f = make_sphere_field(0.5);
938        let r = LevelSetReinit::new(1e-4, 5);
939        let inside_before: Vec<bool> = f.data.iter().map(|&v| v <= 0.0).collect();
940        r.reinitialize(&mut f);
941        let inside_after: Vec<bool> = f.data.iter().map(|&v| v <= 0.0).collect();
942        let changed = inside_before
943            .iter()
944            .zip(inside_after.iter())
945            .filter(|(a, b)| a != b)
946            .count();
947        assert!(
948            changed < f.data.len() / 4,
949            "too many sign changes: {changed}"
950        );
951    }
952    #[test]
953    fn test_reinit_new_parameters() {
954        let r = LevelSetReinit::new(1e-6, 100);
955        assert!((r.tol - 1e-6).abs() < 1e-20);
956        assert_eq!(r.max_iter, 100);
957    }
958    #[test]
959    fn test_fast_sweep_runs_without_panic() {
960        let mut f = make_sphere_field(0.5);
961        let r = LevelSetReinit::new(1e-4, 10);
962        r.fast_sweep(&mut f, 2);
963        assert_eq!(f.data.len(), 8 * 8 * 8);
964    }
965    #[test]
966    fn test_advection_shifts_interface() {
967        let mut f = make_sphere_field(0.5);
968        let evo = LevelSetEvolution::new(0.01, 1);
969        let val_before = f.get(4, 4, 4);
970        evo.advect(&mut f, [1.0, 0.0, 0.0]);
971        let val_after = f.get(4, 4, 4);
972        assert!(val_before.is_finite() && val_after.is_finite());
973        assert!((val_before - val_after).abs() >= 0.0);
974    }
975    #[test]
976    fn test_advection_no_nan() {
977        let mut f = make_sphere_field(0.4);
978        let evo = LevelSetEvolution::new(0.01, 3);
979        evo.advect(&mut f, [0.5, 0.5, 0.0]);
980        assert!(
981            f.data.iter().all(|v| v.is_finite()),
982            "NaN in advected field"
983        );
984    }
985    #[test]
986    fn test_mean_curvature_flow_no_nan() {
987        let mut f = make_sphere_field(0.5);
988        let evo = LevelSetEvolution::new(0.001, 2);
989        evo.mean_curvature_flow(&mut f);
990        assert!(
991            f.data.iter().all(|v| v.is_finite()),
992            "NaN in curvature-flowed field"
993        );
994    }
995    #[test]
996    fn test_normal_velocity_expansion() {
997        let mut f = make_sphere_field(0.5);
998        let evo = LevelSetEvolution::new(0.01, 1);
999        let n_inside_before = f.count_inside();
1000        evo.normal_velocity_extension(&mut f, 1.0);
1001        let n_inside_after = f.count_inside();
1002        assert!(
1003            n_inside_after >= n_inside_before,
1004            "expansion should not shrink interior"
1005        );
1006    }
1007    #[test]
1008    fn test_evolution_new_params() {
1009        let evo = LevelSetEvolution::new(0.05, 10);
1010        assert!((evo.dt - 0.05).abs() < 1e-15);
1011        assert_eq!(evo.n_steps, 10);
1012    }
1013    #[test]
1014    fn test_marching_cubes_sphere_generates_triangles() {
1015        let f = make_sphere_field(0.5);
1016        let ext = IsosurfaceExtraction::new(0.0);
1017        let tris = ext.extract(&f);
1018        assert!(!tris.is_empty(), "sphere should produce some triangles");
1019    }
1020    #[test]
1021    fn test_marching_cubes_all_outside_no_triangles() {
1022        let mut f = LevelSetField::new(4, 4, 4, 0.5, [0.0, 0.0, 0.0]);
1023        for v in f.data.iter_mut() {
1024            *v = 1.0;
1025        }
1026        let ext = IsosurfaceExtraction::new(0.0);
1027        let tris = ext.extract(&f);
1028        assert_eq!(
1029            tris.len(),
1030            0,
1031            "all-outside field should produce no triangles"
1032        );
1033    }
1034    #[test]
1035    fn test_marching_cubes_all_inside_no_triangles() {
1036        let mut f = LevelSetField::new(4, 4, 4, 0.5, [0.0, 0.0, 0.0]);
1037        for v in f.data.iter_mut() {
1038            *v = -1.0;
1039        }
1040        let ext = IsosurfaceExtraction::new(0.0);
1041        let tris = ext.extract(&f);
1042        assert_eq!(
1043            tris.len(),
1044            0,
1045            "all-inside field should produce no triangles"
1046        );
1047    }
1048    #[test]
1049    fn test_marching_cubes_triangle_vertices_finite() {
1050        let f = make_sphere_field(0.5);
1051        let ext = IsosurfaceExtraction::new(0.0);
1052        let tris = ext.extract(&f);
1053        for tri in &tris {
1054            for v in &tri.vertices {
1055                assert!(
1056                    v[0].is_finite() && v[1].is_finite() && v[2].is_finite(),
1057                    "Non-finite vertex: {:?}",
1058                    v
1059                );
1060            }
1061        }
1062    }
1063    #[test]
1064    fn test_marching_cubes_isovalue() {
1065        let f = make_sphere_field(0.5);
1066        let ext0 = IsosurfaceExtraction::new(0.0);
1067        let ext1 = IsosurfaceExtraction::new(0.3);
1068        assert!(ext0.count_triangles(&f) > 0);
1069        let n0 = ext0.count_triangles(&f);
1070        let n1 = ext1.count_triangles(&f);
1071        assert!(n0 > 0, "n0={n0}");
1072        let _ = n1;
1073    }
1074    #[test]
1075    fn test_sdt_unsigned_distance_origin_to_unit_sphere() {
1076        let pts = vec![[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0]];
1077        let sdt = SignedDistanceTransform::new(pts);
1078        let d = sdt.unsigned_distance([0.0, 0.0, 0.0]);
1079        assert!(
1080            (d - 1.0).abs() < 1e-12,
1081            "distance from origin to sphere point: {d}"
1082        );
1083    }
1084    #[test]
1085    fn test_sdt_unsigned_distance_on_source_is_zero() {
1086        let pts = vec![[3.0, 4.0, 0.0]];
1087        let sdt = SignedDistanceTransform::new(pts);
1088        let d = sdt.unsigned_distance([3.0, 4.0, 0.0]);
1089        assert!(
1090            d.abs() < 1e-12,
1091            "distance to source point should be zero: {d}"
1092        );
1093    }
1094    #[test]
1095    fn test_sdt_unsigned_distance_empty_is_max() {
1096        let sdt = SignedDistanceTransform::new(vec![]);
1097        assert_eq!(sdt.unsigned_distance([0.0, 0.0, 0.0]), f64::MAX);
1098    }
1099    #[test]
1100    fn test_sdt_compute_grid_unsigned() {
1101        let pts = vec![[0.0, 0.0, 0.0]];
1102        let sdt = SignedDistanceTransform::new(pts);
1103        let field = sdt.compute_grid_unsigned(4, 4, 4, 0.5, [0.0, 0.0, 0.0]);
1104        let d = field.get(0, 0, 0);
1105        assert!(d.abs() < 1e-12, "distance at origin: {d}");
1106        let dc = field.get(3, 3, 3);
1107        assert!(dc > 0.0, "corner distance: {dc}");
1108    }
1109    #[test]
1110    fn test_sdt_compute_grid_signed() {
1111        let pts: Vec<[f64; 3]> = (0..16)
1112            .flat_map(|i| {
1113                (0..16).map(move |j| {
1114                    let theta = i as f64 * std::f64::consts::PI / 8.0;
1115                    let phi = j as f64 * 2.0 * std::f64::consts::PI / 16.0;
1116                    [theta.cos() * phi.sin(), theta.sin() * phi.sin(), phi.cos()]
1117                })
1118            })
1119            .collect();
1120        let sdt = SignedDistanceTransform::new(pts);
1121        let field = sdt.compute_grid_signed(6, 6, 6, 0.5, [-1.5, -1.5, -1.5]);
1122        assert_eq!(field.data.len(), 6 * 6 * 6);
1123    }
1124    #[test]
1125    fn test_sdt_point_to_segment_midpoint() {
1126        let d = SignedDistanceTransform::point_to_segment(
1127            [0.0, 1.0, 0.0],
1128            [0.0, 0.0, 0.0],
1129            [1.0, 0.0, 0.0],
1130        );
1131        assert!((d - 1.0).abs() < 1e-10, "distance = {d}");
1132    }
1133    #[test]
1134    fn test_sdt_point_to_segment_on_segment() {
1135        let d = SignedDistanceTransform::point_to_segment(
1136            [0.5, 0.0, 0.0],
1137            [0.0, 0.0, 0.0],
1138            [1.0, 0.0, 0.0],
1139        );
1140        assert!(d.abs() < 1e-12, "point on segment: {d}");
1141    }
1142    #[test]
1143    fn test_boolean_union_smaller_sdf() {
1144        let f1 = make_sphere_field(0.5);
1145        let f2 = make_sphere_field(0.3);
1146        let u = LevelSetBoolean::union(&f1, &f2);
1147        for i in 0..u.data.len() {
1148            assert!(u.data[i] <= f1.data[i] + 1e-12);
1149            assert!(u.data[i] <= f2.data[i] + 1e-12);
1150        }
1151    }
1152    #[test]
1153    fn test_boolean_intersection_larger_sdf() {
1154        let f1 = make_sphere_field(0.5);
1155        let f2 = make_sphere_field(0.3);
1156        let inter = LevelSetBoolean::intersection(&f1, &f2);
1157        for i in 0..inter.data.len() {
1158            assert!(inter.data[i] >= f1.data[i] - 1e-12);
1159            assert!(inter.data[i] >= f2.data[i] - 1e-12);
1160        }
1161    }
1162    #[test]
1163    fn test_boolean_difference_outside_both() {
1164        let f1 = make_sphere_field(0.5);
1165        let f2 = make_sphere_field(0.2);
1166        let diff = LevelSetBoolean::difference(&f1, &f2);
1167        assert_eq!(diff.data.len(), f1.data.len());
1168    }
1169    #[test]
1170    fn test_boolean_complement_flips_sign() {
1171        let f = make_sphere_field(0.5);
1172        let comp = LevelSetBoolean::complement(&f);
1173        for i in 0..f.data.len() {
1174            assert!(
1175                (comp.data[i] + f.data[i]).abs() < 1e-12,
1176                "complement failed at {i}"
1177            );
1178        }
1179    }
1180    #[test]
1181    fn test_boolean_smooth_union_between_min_and_max() {
1182        let f1 = make_sphere_field(0.5);
1183        let f2 = make_sphere_field(0.3);
1184        let su = LevelSetBoolean::smooth_union(&f1, &f2, 0.1);
1185        let regular = LevelSetBoolean::union(&f1, &f2);
1186        for i in 0..su.data.len() {
1187            assert!(
1188                su.data[i] <= regular.data[i] + 0.2,
1189                "smooth union too large at {i}"
1190            );
1191        }
1192    }
1193    #[test]
1194    fn test_boolean_union_idempotent() {
1195        let f = make_sphere_field(0.5);
1196        let u = LevelSetBoolean::union(&f, &f);
1197        for i in 0..f.data.len() {
1198            assert!(
1199                (u.data[i] - f.data[i]).abs() < 1e-12,
1200                "union with self should be identity"
1201            );
1202        }
1203    }
1204    #[test]
1205    fn test_boolean_intersection_idempotent() {
1206        let f = make_sphere_field(0.5);
1207        let inter = LevelSetBoolean::intersection(&f, &f);
1208        for i in 0..f.data.len() {
1209            assert!(
1210                (inter.data[i] - f.data[i]).abs() < 1e-12,
1211                "intersection with self should be identity"
1212            );
1213        }
1214    }
1215    #[test]
1216    fn test_boolean_difference_self_is_empty() {
1217        let f = make_sphere_field(0.5);
1218        let diff = LevelSetBoolean::difference(&f, &f);
1219        for &v in &diff.data {
1220            assert!(
1221                v >= -1e-12,
1222                "difference with self should be outside everywhere: {v}"
1223            );
1224        }
1225    }
1226    #[test]
1227    fn test_interp_vertex_midpoint() {
1228        let p = interp_vertex([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], -1.0, 1.0, 0.0);
1229        assert!((p[0] - 0.5).abs() < 1e-12, "midpoint x: {}", p[0]);
1230    }
1231    #[test]
1232    fn test_interp_vertex_at_a() {
1233        let p = interp_vertex([1.0, 2.0, 3.0], [4.0, 5.0, 6.0], 0.0, 1.0, 0.0);
1234        assert!((p[0] - 1.0).abs() < 1e-12);
1235        assert!((p[1] - 2.0).abs() < 1e-12);
1236        assert!((p[2] - 3.0).abs() < 1e-12);
1237    }
1238    #[test]
1239    fn test_interp_vertex_equal_values_returns_midpoint() {
1240        let p = interp_vertex([0.0, 0.0, 0.0], [2.0, 0.0, 0.0], 1.0, 1.0, 0.5);
1241        assert!(p[0].is_finite());
1242    }
1243}
oxiphysics_geometry/level_set/functions.rs

oxiphysics_geometry/level_set/
functions.rs
