oxiphysics-geometry 0.1.2

//! Auto-generated module
//!
//! 🤖 Generated with [SplitRS](https://github.com/cool-japan/splitrs)

use super::types::LevelSetField;

/// Godunov upwind gradient magnitude for reinitialization.
///
/// Computes the Godunov upwind approximation of |∇φ| needed for the
/// Sussman reinitialization scheme.
pub(super) fn godunov_grad_mag(
    field: &LevelSetField,
    ix: usize,
    iy: usize,
    iz: usize,
    sign: f64,
) -> f64 {
    let dx = field.dx;
    let phi_c = field.get(ix, iy, iz);
    let d_xm = if ix > 0 {
        (phi_c - field.get(ix - 1, iy, iz)) / dx
    } else {
        0.0
    };
    let d_xp = if ix + 1 < field.nx {
        (field.get(ix + 1, iy, iz) - phi_c) / dx
    } else {
        0.0
    };
    let d_ym = if iy > 0 {
        (phi_c - field.get(ix, iy - 1, iz)) / dx
    } else {
        0.0
    };
    let d_yp = if iy + 1 < field.ny {
        (field.get(ix, iy + 1, iz) - phi_c) / dx
    } else {
        0.0
    };
    let d_zm = if iz > 0 {
        (phi_c - field.get(ix, iy, iz - 1)) / dx
    } else {
        0.0
    };
    let d_zp = if iz + 1 < field.nz {
        (field.get(ix, iy, iz + 1) - phi_c) / dx
    } else {
        0.0
    };
    let gx = if sign > 0.0 {
        d_xm.max(0.0).powi(2).max((-d_xp).max(0.0).powi(2))
    } else {
        (-d_xm).max(0.0).powi(2).max(d_xp.max(0.0).powi(2))
    };
    let gy = if sign > 0.0 {
        d_ym.max(0.0).powi(2).max((-d_yp).max(0.0).powi(2))
    } else {
        (-d_ym).max(0.0).powi(2).max(d_yp.max(0.0).powi(2))
    };
    let gz = if sign > 0.0 {
        d_zm.max(0.0).powi(2).max((-d_zp).max(0.0).powi(2))
    } else {
        (-d_zm).max(0.0).powi(2).max(d_zp.max(0.0).powi(2))
    };
    (gx + gy + gz).sqrt()
}
/// Linear interpolation of vertex position on a cube edge.
pub(super) fn interp_vertex(pa: [f64; 3], pb: [f64; 3], va: f64, vb: f64, iso: f64) -> [f64; 3] {
    let dv = vb - va;
    let t = if dv.abs() < 1e-30 {
        0.5
    } else {
        (iso - va) / dv
    };
    [
        pa[0] + t * (pb[0] - pa[0]),
        pa[1] + t * (pb[1] - pa[1]),
        pa[2] + t * (pb[2] - pa[2]),
    ]
}
/// Edges of a unit cube: (corner_a, corner_b) pairs for the 12 edges.
pub(super) const MC_EDGE_CORNERS: [(usize, usize); 12] = [
    (0, 1),
    (1, 2),
    (2, 3),
    (3, 0),
    (4, 5),
    (5, 6),
    (6, 7),
    (7, 4),
    (0, 4),
    (1, 5),
    (2, 6),
    (3, 7),
];
pub(super) const MC_EDGE_TABLE: [u16; 256] = [
    0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a,
    0xd03, 0xe09, 0xf00, 0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895,
    0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, 0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435,
    0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 0x3a0, 0x2a9, 0x1a3, 0x0aa,
    0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, 0x460,
    0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963,
    0xa69, 0xb60, 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff,
    0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c,
    0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6,
    0x2cf, 0x1c5, 0x0cc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, 0x8c0, 0x9c9,
    0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0x0cc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9,
    0x7c0, 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x055, 0x35f, 0x256,
    0x55a, 0x453, 0x759, 0x650, 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc,
    0x3f5, 0x0ff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f,
    0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x066, 0x76a, 0x663, 0x569, 0x460, 0xca0, 0xda9, 0xea3,
    0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0x0aa, 0x1a3, 0x2a9, 0x3a0,
    0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x835, 0xb3f, 0xa36, 0x53c, 0x435, 0x73f, 0x636, 0x13a,
    0x033, 0x339, 0x230, 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795,
    0x49f, 0x596, 0x29a, 0x393, 0x099, 0x190, 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905,
    0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x000,
];
/// Marching cubes triangle table (Lorensen & Cline, 1987).
///
/// For each of the 256 cube configurations, lists the edge indices of the
/// triangles to generate.  Each row ends with -1 as sentinel.
/// (Abbreviated to key configurations; full 256-entry tables follow the same pattern.)
pub(super) const MC_TRI_TABLE: [[i8; 16]; 256] = [
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1],
    [3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1],
    [3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1],
    [3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1],
    [9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1],
    [1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1],
    [9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1],
    [2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1],
    [8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
    [11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1],
    [9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1],
    [4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1],
    [3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1],
    [1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1],
    [4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1],
    [4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
    [
        -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
    ],
];
#[cfg(test)]
mod tests {
    use super::*;
    use crate::level_set::IsosurfaceExtraction;
    use crate::level_set::LevelSetBoolean;
    use crate::level_set::LevelSetEvolution;
    use crate::level_set::LevelSetReinit;
    use crate::level_set::SignedDistanceTransform;
    fn make_sphere_field(radius: f64) -> LevelSetField {
        let mut f = LevelSetField::new(8, 8, 8, 0.25, [-1.0, -1.0, -1.0]);
        f.init_sphere([0.0, 0.0, 0.0], radius);
        f
    }
    #[test]
    fn test_field_new_dimensions() {
        let f = LevelSetField::new(4, 5, 6, 0.1, [0.0, 0.0, 0.0]);
        assert_eq!(f.nx, 4);
        assert_eq!(f.ny, 5);
        assert_eq!(f.nz, 6);
        assert_eq!(f.data.len(), 4 * 5 * 6);
    }
    #[test]
    fn test_field_set_get() {
        let mut f = LevelSetField::new(4, 4, 4, 0.1, [0.0, 0.0, 0.0]);
        f.set(1, 2, 3, -0.5);
        assert!((f.get(1, 2, 3) + 0.5).abs() < 1e-12);
    }
    #[test]
    fn test_field_idx_linear() {
        let f = LevelSetField::new(4, 4, 4, 0.1, [0.0, 0.0, 0.0]);
        let idx = f.idx(1, 2, 3);
        assert_eq!(idx, 57);
    }
    #[test]
    fn test_field_node_pos() {
        let f = LevelSetField::new(5, 5, 5, 0.5, [0.0, 0.0, 0.0]);
        let p = f.node_pos(2, 0, 0);
        assert!((p[0] - 1.0).abs() < 1e-12);
        assert!(p[1].abs() < 1e-12);
    }
    #[test]
    fn test_field_out_of_bounds_get() {
        let f = LevelSetField::new(4, 4, 4, 0.1, [0.0, 0.0, 0.0]);
        assert_eq!(f.get(10, 0, 0), f64::MAX);
    }
    #[test]
    fn test_field_out_of_bounds_set_noop() {
        let mut f = LevelSetField::new(4, 4, 4, 0.1, [0.0, 0.0, 0.0]);
        f.set(100, 0, 0, 42.0);
        assert_eq!(f.data.iter().filter(|&&v| v == 42.0).count(), 0);
    }
    #[test]
    fn test_field_init_sphere_center_negative() {
        let f = make_sphere_field(0.5);
        let center_val = f.interpolate([0.0, 0.0, 0.0]);
        assert!(
            center_val < 0.0,
            "center should be inside sphere: {center_val}"
        );
    }
    #[test]
    fn test_field_init_sphere_far_point_positive() {
        let f = make_sphere_field(0.5);
        let far_val = f.get(0, 0, 0);
        assert!(far_val > 0.0, "corner should be outside sphere: {far_val}");
    }
    #[test]
    fn test_field_init_box_inside_negative() {
        let mut f = LevelSetField::new(8, 8, 8, 0.25, [-1.0, -1.0, -1.0]);
        f.init_box([-0.3, -0.3, -0.3], [0.3, 0.3, 0.3]);
        let val = f.interpolate([0.0, 0.0, 0.0]);
        assert!(val < 0.0, "inside box should be negative: {val}");
    }
    #[test]
    fn test_field_count_inside_sphere() {
        let f = make_sphere_field(0.5);
        let n_inside = f.count_inside();
        assert!(n_inside > 0, "some nodes should be inside sphere");
        assert!(n_inside < f.data.len(), "not all nodes inside sphere");
    }
    #[test]
    fn test_field_min_max_sphere() {
        let f = make_sphere_field(0.5);
        let vmin = f.min_value();
        let vmax = f.max_value();
        assert!(vmin < 0.0, "min should be negative: {vmin}");
        assert!(vmax > 0.0, "max should be positive: {vmax}");
    }
    #[test]
    fn test_field_gradient_sphere_on_surface() {
        let f = make_sphere_field(0.5);
        let mid = 4usize;
        let g = f.gradient(mid, mid, mid);
        let mag = (g[0] * g[0] + g[1] * g[1] + g[2] * g[2]).sqrt();
        assert!(mag.is_finite(), "gradient magnitude should be finite");
    }
    #[test]
    fn test_field_interpolate_on_node() {
        let mut f = LevelSetField::new(5, 5, 5, 0.5, [0.0, 0.0, 0.0]);
        f.set(1, 1, 1, 3.125);
        let p = f.node_pos(1, 1, 1);
        let val = f.interpolate(p);
        assert!(val.is_finite());
    }
    #[test]
    fn test_reinit_sphere_preserves_zero_set() {
        let mut f = make_sphere_field(0.5);
        let r = LevelSetReinit::new(1e-4, 5);
        let inside_before: Vec<bool> = f.data.iter().map(|&v| v <= 0.0).collect();
        r.reinitialize(&mut f);
        let inside_after: Vec<bool> = f.data.iter().map(|&v| v <= 0.0).collect();
        let changed = inside_before
            .iter()
            .zip(inside_after.iter())
            .filter(|(a, b)| a != b)
            .count();
        assert!(
            changed < f.data.len() / 4,
            "too many sign changes: {changed}"
        );
    }
    #[test]
    fn test_reinit_new_parameters() {
        let r = LevelSetReinit::new(1e-6, 100);
        assert!((r.tol - 1e-6).abs() < 1e-20);
        assert_eq!(r.max_iter, 100);
    }
    #[test]
    fn test_fast_sweep_runs_without_panic() {
        let mut f = make_sphere_field(0.5);
        let r = LevelSetReinit::new(1e-4, 10);
        r.fast_sweep(&mut f, 2);
        assert_eq!(f.data.len(), 8 * 8 * 8);
    }
    #[test]
    fn test_advection_shifts_interface() {
        let mut f = make_sphere_field(0.5);
        let evo = LevelSetEvolution::new(0.01, 1);
        let val_before = f.get(4, 4, 4);
        evo.advect(&mut f, [1.0, 0.0, 0.0]);
        let val_after = f.get(4, 4, 4);
        assert!(val_before.is_finite() && val_after.is_finite());
        assert!((val_before - val_after).abs() >= 0.0);
    }
    #[test]
    fn test_advection_no_nan() {
        let mut f = make_sphere_field(0.4);
        let evo = LevelSetEvolution::new(0.01, 3);
        evo.advect(&mut f, [0.5, 0.5, 0.0]);
        assert!(
            f.data.iter().all(|v| v.is_finite()),
            "NaN in advected field"
        );
    }
    #[test]
    fn test_mean_curvature_flow_no_nan() {
        let mut f = make_sphere_field(0.5);
        let evo = LevelSetEvolution::new(0.001, 2);
        evo.mean_curvature_flow(&mut f);
        assert!(
            f.data.iter().all(|v| v.is_finite()),
            "NaN in curvature-flowed field"
        );
    }
    #[test]
    fn test_normal_velocity_expansion() {
        let mut f = make_sphere_field(0.5);
        let evo = LevelSetEvolution::new(0.01, 1);
        let n_inside_before = f.count_inside();
        evo.normal_velocity_extension(&mut f, 1.0);
        let n_inside_after = f.count_inside();
        assert!(
            n_inside_after >= n_inside_before,
            "expansion should not shrink interior"
        );
    }
    #[test]
    fn test_evolution_new_params() {
        let evo = LevelSetEvolution::new(0.05, 10);
        assert!((evo.dt - 0.05).abs() < 1e-15);
        assert_eq!(evo.n_steps, 10);
    }
    #[test]
    fn test_marching_cubes_sphere_generates_triangles() {
        let f = make_sphere_field(0.5);
        let ext = IsosurfaceExtraction::new(0.0);
        let tris = ext.extract(&f);
        assert!(!tris.is_empty(), "sphere should produce some triangles");
    }
    #[test]
    fn test_marching_cubes_all_outside_no_triangles() {
        let mut f = LevelSetField::new(4, 4, 4, 0.5, [0.0, 0.0, 0.0]);
        for v in f.data.iter_mut() {
            *v = 1.0;
        }
        let ext = IsosurfaceExtraction::new(0.0);
        let tris = ext.extract(&f);
        assert_eq!(
            tris.len(),
            0,
            "all-outside field should produce no triangles"
        );
    }
    #[test]
    fn test_marching_cubes_all_inside_no_triangles() {
        let mut f = LevelSetField::new(4, 4, 4, 0.5, [0.0, 0.0, 0.0]);
        for v in f.data.iter_mut() {
            *v = -1.0;
        }
        let ext = IsosurfaceExtraction::new(0.0);
        let tris = ext.extract(&f);
        assert_eq!(
            tris.len(),
            0,
            "all-inside field should produce no triangles"
        );
    }
    #[test]
    fn test_marching_cubes_triangle_vertices_finite() {
        let f = make_sphere_field(0.5);
        let ext = IsosurfaceExtraction::new(0.0);
        let tris = ext.extract(&f);
        for tri in &tris {
            for v in &tri.vertices {
                assert!(
                    v[0].is_finite() && v[1].is_finite() && v[2].is_finite(),
                    "Non-finite vertex: {:?}",
                    v
                );
            }
        }
    }
    #[test]
    fn test_marching_cubes_isovalue() {
        let f = make_sphere_field(0.5);
        let ext0 = IsosurfaceExtraction::new(0.0);
        let ext1 = IsosurfaceExtraction::new(0.3);
        assert!(ext0.count_triangles(&f) > 0);
        let n0 = ext0.count_triangles(&f);
        let n1 = ext1.count_triangles(&f);
        assert!(n0 > 0, "n0={n0}");
        let _ = n1;
    }
    #[test]
    fn test_sdt_unsigned_distance_origin_to_unit_sphere() {
        let pts = vec![[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0]];
        let sdt = SignedDistanceTransform::new(pts);
        let d = sdt.unsigned_distance([0.0, 0.0, 0.0]);
        assert!(
            (d - 1.0).abs() < 1e-12,
            "distance from origin to sphere point: {d}"
        );
    }
    #[test]
    fn test_sdt_unsigned_distance_on_source_is_zero() {
        let pts = vec![[3.0, 4.0, 0.0]];
        let sdt = SignedDistanceTransform::new(pts);
        let d = sdt.unsigned_distance([3.0, 4.0, 0.0]);
        assert!(
            d.abs() < 1e-12,
            "distance to source point should be zero: {d}"
        );
    }
    #[test]
    fn test_sdt_unsigned_distance_empty_is_max() {
        let sdt = SignedDistanceTransform::new(vec![]);
        assert_eq!(sdt.unsigned_distance([0.0, 0.0, 0.0]), f64::MAX);
    }
    #[test]
    fn test_sdt_compute_grid_unsigned() {
        let pts = vec![[0.0, 0.0, 0.0]];
        let sdt = SignedDistanceTransform::new(pts);
        let field = sdt.compute_grid_unsigned(4, 4, 4, 0.5, [0.0, 0.0, 0.0]);
        let d = field.get(0, 0, 0);
        assert!(d.abs() < 1e-12, "distance at origin: {d}");
        let dc = field.get(3, 3, 3);
        assert!(dc > 0.0, "corner distance: {dc}");
    }
    #[test]
    fn test_sdt_compute_grid_signed() {
        let pts: Vec<[f64; 3]> = (0..16)
            .flat_map(|i| {
                (0..16).map(move |j| {
                    let theta = i as f64 * std::f64::consts::PI / 8.0;
                    let phi = j as f64 * 2.0 * std::f64::consts::PI / 16.0;
                    [theta.cos() * phi.sin(), theta.sin() * phi.sin(), phi.cos()]
                })
            })
            .collect();
        let sdt = SignedDistanceTransform::new(pts);
        let field = sdt.compute_grid_signed(6, 6, 6, 0.5, [-1.5, -1.5, -1.5]);
        assert_eq!(field.data.len(), 6 * 6 * 6);
    }
    #[test]
    fn test_sdt_point_to_segment_midpoint() {
        let d = SignedDistanceTransform::point_to_segment(
            [0.0, 1.0, 0.0],
            [0.0, 0.0, 0.0],
            [1.0, 0.0, 0.0],
        );
        assert!((d - 1.0).abs() < 1e-10, "distance = {d}");
    }
    #[test]
    fn test_sdt_point_to_segment_on_segment() {
        let d = SignedDistanceTransform::point_to_segment(
            [0.5, 0.0, 0.0],
            [0.0, 0.0, 0.0],
            [1.0, 0.0, 0.0],
        );
        assert!(d.abs() < 1e-12, "point on segment: {d}");
    }
    #[test]
    fn test_boolean_union_smaller_sdf() {
        let f1 = make_sphere_field(0.5);
        let f2 = make_sphere_field(0.3);
        let u = LevelSetBoolean::union(&f1, &f2);
        for i in 0..u.data.len() {
            assert!(u.data[i] <= f1.data[i] + 1e-12);
            assert!(u.data[i] <= f2.data[i] + 1e-12);
        }
    }
    #[test]
    fn test_boolean_intersection_larger_sdf() {
        let f1 = make_sphere_field(0.5);
        let f2 = make_sphere_field(0.3);
        let inter = LevelSetBoolean::intersection(&f1, &f2);
        for i in 0..inter.data.len() {
            assert!(inter.data[i] >= f1.data[i] - 1e-12);
            assert!(inter.data[i] >= f2.data[i] - 1e-12);
        }
    }
    #[test]
    fn test_boolean_difference_outside_both() {
        let f1 = make_sphere_field(0.5);
        let f2 = make_sphere_field(0.2);
        let diff = LevelSetBoolean::difference(&f1, &f2);
        assert_eq!(diff.data.len(), f1.data.len());
    }
    #[test]
    fn test_boolean_complement_flips_sign() {
        let f = make_sphere_field(0.5);
        let comp = LevelSetBoolean::complement(&f);
        for i in 0..f.data.len() {
            assert!(
                (comp.data[i] + f.data[i]).abs() < 1e-12,
                "complement failed at {i}"
            );
        }
    }
    #[test]
    fn test_boolean_smooth_union_between_min_and_max() {
        let f1 = make_sphere_field(0.5);
        let f2 = make_sphere_field(0.3);
        let su = LevelSetBoolean::smooth_union(&f1, &f2, 0.1);
        let regular = LevelSetBoolean::union(&f1, &f2);
        for i in 0..su.data.len() {
            assert!(
                su.data[i] <= regular.data[i] + 0.2,
                "smooth union too large at {i}"
            );
        }
    }
    #[test]
    fn test_boolean_union_idempotent() {
        let f = make_sphere_field(0.5);
        let u = LevelSetBoolean::union(&f, &f);
        for i in 0..f.data.len() {
            assert!(
                (u.data[i] - f.data[i]).abs() < 1e-12,
                "union with self should be identity"
            );
        }
    }
    #[test]
    fn test_boolean_intersection_idempotent() {
        let f = make_sphere_field(0.5);
        let inter = LevelSetBoolean::intersection(&f, &f);
        for i in 0..f.data.len() {
            assert!(
                (inter.data[i] - f.data[i]).abs() < 1e-12,
                "intersection with self should be identity"
            );
        }
    }
    #[test]
    fn test_boolean_difference_self_is_empty() {
        let f = make_sphere_field(0.5);
        let diff = LevelSetBoolean::difference(&f, &f);
        for &v in &diff.data {
            assert!(
                v >= -1e-12,
                "difference with self should be outside everywhere: {v}"
            );
        }
    }
    #[test]
    fn test_interp_vertex_midpoint() {
        let p = interp_vertex([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], -1.0, 1.0, 0.0);
        assert!((p[0] - 0.5).abs() < 1e-12, "midpoint x: {}", p[0]);
    }
    #[test]
    fn test_interp_vertex_at_a() {
        let p = interp_vertex([1.0, 2.0, 3.0], [4.0, 5.0, 6.0], 0.0, 1.0, 0.0);
        assert!((p[0] - 1.0).abs() < 1e-12);
        assert!((p[1] - 2.0).abs() < 1e-12);
        assert!((p[2] - 3.0).abs() < 1e-12);
    }
    #[test]
    fn test_interp_vertex_equal_values_returns_midpoint() {
        let p = interp_vertex([0.0, 0.0, 0.0], [2.0, 0.0, 0.0], 1.0, 1.0, 0.5);
        assert!(p[0].is_finite());
    }
}