//! The cut geometry of an embedded body on the grid: φ at the cell
//! corners; inside every cell the interface is the LINEAR interpolant on a
//! fixed Kuhn split into six tetrahedra (each face into two triangles along
//! the same diagonal from both sides), so face apertures and cell volumes
//! are exact for the interpolant, continuous in the corner values, and
//! consistent across shared faces — the wall polygon's vector area by
//! closure (`A_w n_w = −Σ_f A_f n_f`) telescopes exactly over a closed body.

use super::Grid;
use super::body::Body;

/// The cut data of one instant.
#[derive(Debug, Clone)]
pub struct CutGeometry {
    pub grid: Grid,
    /// φ at the corners, `(nx + 1) × (ny + 1) × (nz + 1)`.
    pub phi: Vec<f64>,
    /// Fluid area fraction of every u / v / w face (the staggered layouts).
    pub a_u: Vec<f64>,
    pub a_v: Vec<f64>,
    pub a_w: Vec<f64>,
    /// Fluid volume fraction of every cell.
    pub vol: Vec<f64>,
    /// The wall polygon's vector area per cell (outward from the fluid),
    /// `−Σ_f A_f n_f · face area`, in area units.
    pub wall: Vec<[f64; 3]>,
    /// φ at the face centres (the wall distance of a near-wall face).
    pub d_u: Vec<f64>,
    pub d_v: Vec<f64>,
    pub d_w: Vec<f64>,
    /// A lower bound on |φ| per corner (the narrow band of a moving body:
    /// a corner is re-evaluated only when its bound, decayed by the body's
    /// motion, comes within the band; far corners keep a stale value with
    /// the right sign, which is all their cells use).
    pub bound: Vec<f64>,
    /// P1-4: the corners re-evaluated by this build (all of them without a
    /// narrow band) — a cell whose corners are all untouched has unchanged
    /// apertures, volume and activity.
    pub touched: Vec<bool>,
    /// Per cell: any of its eight corners touched (computed once per build).
    pub touched_cell: Vec<bool>,
}

impl CutGeometry {
    #[inline]
    fn node(g: Grid, k: usize, j: usize, i: usize) -> usize {
        (k * (g.ny + 1) + j) * (g.nx + 1) + i
    }

    /// Build the cut data of `body` at time `t` (every corner evaluated).
    pub fn build(body: &Body, grid: Grid, t: f64) -> Self {
        Self::build_from(body, grid, t, None)
    }

    /// Build the cut data of `body` at `t`, re-evaluating only the corners
    /// of `prev` whose |φ| bound, decayed by `motion` (the body's largest
    /// displacement since `prev`), falls within `band` of the surface.
    /// Identical to [`Self::build`] in every cut cell.
    pub fn build_from(
        body: &Body,
        grid: Grid,
        t: f64,
        prev: Option<(&CutGeometry, f64, f64)>,
    ) -> Self {
        let g = grid;
        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
        let n_nodes = (nx + 1) * (ny + 1) * (nz + 1);
        // Every loop below is per-index with read-only inputs: rayon (S2-2b),
        // the same arithmetic per entry.
        use rayon::prelude::*;
        let mut phi = vec![0.0; n_nodes];
        let mut bound = vec![0.0; n_nodes];
        let mut touched = vec![true; n_nodes];
        phi.par_iter_mut()
            .zip(bound.par_iter_mut())
            .zip(touched.par_iter_mut())
            .enumerate()
            .for_each(|(n, ((phi_n, bound_n), touched_n))| {
                let (k, j, i) = (
                    n / ((ny + 1) * (nx + 1)),
                    (n / (nx + 1)) % (ny + 1),
                    n % (nx + 1),
                );
                if let Some((p, band, motion)) = prev {
                    let b = p.bound[n] - motion;
                    if b > band {
                        *phi_n = p.phi[n];
                        *bound_n = b;
                        *touched_n = false;
                        return;
                    }
                }
                let v = body.phi(i as f64 * dx, j as f64 * dy, k as f64 * dz, t);
                *phi_n = v;
                *bound_n = v.abs();
            });
        let corner = |k: usize, j: usize, i: usize| phi[Self::node(g, k, j, i)];
        // Face apertures: a face's two triangles along the diagonal from its
        // (0, 0) to its (1, 1) corner in the face's own (a, b) order — the
        // Kuhn split's diagonals: for an x-face (y, z), a y-face (x, z), a
        // z-face (x, y); the same triangles seen from either cell.
        let quad_fraction = quad_fraction;
        let mut a_u = vec![0.0; (nx + 1) * ny * nz];
        let mut a_v = vec![0.0; nx * (ny + 1) * nz];
        let mut a_w = vec![0.0; nx * ny * (nz + 1)];
        let mut d_u = vec![0.0; (nx + 1) * ny * nz];
        let mut d_v = vec![0.0; nx * (ny + 1) * nz];
        let mut d_w = vec![0.0; nx * ny * (nz + 1)];
        a_u.par_iter_mut()
            .zip(d_u.par_iter_mut())
            .enumerate()
            .for_each(|(f, (a, d))| {
                // x-face at i: corners (j, k), (j+1, k), (j, k+1), (j+1, k+1)
                let (k, j, i) = (f / (ny * (nx + 1)), (f / (nx + 1)) % ny, f % (nx + 1));
                let (q00, q10, q01, q11) = (
                    corner(k, j, i),
                    corner(k, j + 1, i),
                    corner(k + 1, j, i),
                    corner(k + 1, j + 1, i),
                );
                *a = quad_fraction(q00, q10, q01, q11);
                *d = 0.25 * (q00 + q10 + q01 + q11);
            });
        a_v.par_iter_mut()
            .zip(d_v.par_iter_mut())
            .enumerate()
            .for_each(|(f, (a, d))| {
                // y-face at j: corners (i, k), (i+1, k), (i, k+1), (i+1, k+1)
                let (k, j, i) = (f / ((ny + 1) * nx), (f / nx) % (ny + 1), f % nx);
                let (q00, q10, q01, q11) = (
                    corner(k, j, i),
                    corner(k, j, i + 1),
                    corner(k + 1, j, i),
                    corner(k + 1, j, i + 1),
                );
                *a = quad_fraction(q00, q10, q01, q11);
                *d = 0.25 * (q00 + q10 + q01 + q11);
            });
        a_w.par_iter_mut()
            .zip(d_w.par_iter_mut())
            .enumerate()
            .for_each(|(f, (a, d))| {
                // z-face at k: corners (i, j), (i+1, j), (i, j+1), (i+1, j+1)
                let (k, j, i) = (f / (ny * nx), (f / nx) % ny, f % nx);
                let (q00, q10, q01, q11) = (
                    corner(k, j, i),
                    corner(k, j, i + 1),
                    corner(k, j + 1, i),
                    corner(k, j + 1, i + 1),
                );
                *a = quad_fraction(q00, q10, q01, q11);
                *d = 0.25 * (q00 + q10 + q01 + q11);
            });
        // Cell volumes by the Kuhn split: the six tetrahedra around the
        // diagonal (0,0,0)–(1,1,1) in unit-cube coordinates.
        let mut vol = vec![0.0; g.cells()];
        let mut wall = vec![[0.0; 3]; g.cells()];
        const KUHN: [[[usize; 3]; 4]; 6] = [
            [[0, 0, 0], [1, 0, 0], [1, 1, 0], [1, 1, 1]],
            [[0, 0, 0], [1, 0, 0], [1, 0, 1], [1, 1, 1]],
            [[0, 0, 0], [0, 1, 0], [1, 1, 0], [1, 1, 1]],
            [[0, 0, 0], [0, 1, 0], [0, 1, 1], [1, 1, 1]],
            [[0, 0, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1]],
            [[0, 0, 0], [0, 0, 1], [0, 1, 1], [1, 1, 1]],
        ];
        let cell_fluid = |k: usize, j: usize, i: usize| -> f64 {
            let mut fluid = 0.0;
            for tet in &KUHN {
                let pts: Vec<[f64; 3]> = tet
                    .iter()
                    .map(|c| [c[0] as f64, c[1] as f64, c[2] as f64])
                    .collect();
                let vals: Vec<f64> = tet
                    .iter()
                    .map(|c| corner(k + c[2], j + c[1], i + c[0]))
                    .collect();
                fluid += tet_fluid_volume(&pts, &vals);
            }
            fluid
        };
        let (ax, ay, az) = (dy * dz, dx * dz, dx * dy);
        vol.par_iter_mut()
            .zip(wall.par_iter_mut())
            .enumerate()
            .for_each(|(idx, (v, w))| {
                let (k, j, i) = (idx / (ny * nx), (idx / nx) % ny, idx % nx);
                // The six tets fill the unit cube (volume 1).
                *v = cell_fluid(k, j, i);
                // Outward normals of the cell's faces times their fluid area,
                // summed; the wall closes the fluid part of the cell.
                let sx = (a_u[g.uface(k, j, i + 1)] - a_u[g.uface(k, j, i)]) * ax;
                let sy = (a_v[g.vface(k, j + 1, i)] - a_v[g.vface(k, j, i)]) * ay;
                let sz = (a_w[g.wface(k + 1, j, i)] - a_w[g.wface(k, j, i)]) * az;
                *w = [-sx, -sy, -sz];
            });
        let touched_cell: Vec<bool> = (0..g.cells())
            .into_par_iter()
            .map(|idx| {
                let (k, j, i) = (idx / (ny * nx), (idx / nx) % ny, idx % nx);
                let mut t = false;
                for dk in 0..2 {
                    for dj in 0..2 {
                        for di in 0..2 {
                            t |= touched[Self::node(g, k + dk, j + dj, i + di)];
                        }
                    }
                }
                t
            })
            .collect();
        Self {
            grid,
            phi,
            a_u,
            a_v,
            a_w,
            vol,
            wall,
            d_u,
            d_v,
            d_w,
            bound,
            touched,
            touched_cell,
        }
    }

    /// The cells with a touched corner (P1-4; the stored map).
    #[must_use]
    pub fn touched_cells(&self) -> Vec<bool> {
        if !self.touched_cell.is_empty() {
            return self.touched_cell.clone();
        }
        use rayon::prelude::*;
        let g = self.grid;
        let (nx, ny) = (g.nx, g.ny);
        (0..g.cells())
            .into_par_iter()
            .map(|idx| {
                let (k, j, i) = g.kji(idx);
                let mut t = false;
                for dk in 0..2 {
                    for dj in 0..2 {
                        for di in 0..2 {
                            t |= self.touched[Self::node(g, k + dk, j + dj, i + di)];
                        }
                    }
                }
                let _ = (nx, ny);
                t
            })
            .collect()
    }

    /// φ at the corner `(k, j, i)` of the corner lattice.
    #[inline]
    #[must_use]
    pub fn corner_phi(&self, k: usize, j: usize, i: usize) -> f64 {
        self.phi[Self::node(self.grid, k, j, i)]
    }

    /// Total fluid volume.
    #[must_use]
    pub fn fluid_volume(&self) -> f64 {
        let g = self.grid;
        self.vol.iter().sum::<f64>() * g.dx * g.dy * g.dz
    }

    /// Σ over cells of |wall vector area| (the wall's area up to the
    /// non-planarity of the per-cell polygon) and the closure vector Σ wall.
    #[must_use]
    pub fn wall_area_and_closure(&self) -> (f64, [f64; 3]) {
        let mut area = 0.0;
        let mut sum = [0.0; 3];
        for w in &self.wall {
            area += (w[0] * w[0] + w[1] * w[1] + w[2] * w[2]).sqrt();
            sum[0] += w[0];
            sum[1] += w[1];
            sum[2] += w[2];
        }
        (area, sum)
    }
}

/// The fluid fraction of a triangle from its three corner values of φ
/// (the linear interpolant; fluid where φ ≥ 0).
pub(super) fn tri_area_fraction(p0: f64, p1: f64, p2: f64) -> f64 {
    let v = [p0, p1, p2];
    let pos = v.iter().filter(|&&q| q >= 0.0).count();
    match pos {
        0 => 0.0,
        3 => 1.0,
        1 => {
            let a = v.iter().position(|&q| q >= 0.0).unwrap();
            let (b, c) = ((a + 1) % 3, (a + 2) % 3);
            (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c]))
        }
        _ => {
            let a = v.iter().position(|&q| q < 0.0).unwrap();
            let (b, c) = ((a + 1) % 3, (a + 2) % 3);
            1.0 - (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c]))
        }
    }
}

/// The fluid fraction of a quad from its corner values in (a, b) order
/// (q00, q10, q01, q11): the two triangles along the (0, 0)–(1, 1)
/// diagonal (the Kuhn split's).
pub(super) fn quad_fraction(q00: f64, q10: f64, q01: f64, q11: f64) -> f64 {
    0.5 * (tri_area_fraction(q00, q10, q11) + tri_area_fraction(q00, q11, q01))
}

fn det3(a: [f64; 3], b: [f64; 3], c: [f64; 3]) -> f64 {
    a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0])
        + a[2] * (b[0] * c[1] - b[1] * c[0])
}

fn tet_volume(p: [[f64; 3]; 4]) -> f64 {
    let e = |a: [f64; 3], b: [f64; 3]| [b[0] - a[0], b[1] - a[1], b[2] - a[2]];
    det3(e(p[0], p[1]), e(p[0], p[2]), e(p[0], p[3])).abs() / 6.0
}

fn lerp(a: [f64; 3], b: [f64; 3], t: f64) -> [f64; 3] {
    [
        a[0] + t * (b[0] - a[0]),
        a[1] + t * (b[1] - a[1]),
        a[2] + t * (b[2] - a[2]),
    ]
}

/// The volume of `{φ ≥ 0}` in a tetrahedron with the linear interpolant of
/// the corner values `v` (φ = 0 at a corner counts as fluid). Continuous in
/// `v`: every case's cut points move continuously and the cases agree on
/// their boundaries.
fn tet_fluid_volume(pts: &[[f64; 3]], v: &[f64]) -> f64 {
    let total = tet_volume([pts[0], pts[1], pts[2], pts[3]]);
    let pos: Vec<usize> = (0..4).filter(|&q| v[q] >= 0.0).collect();
    let neg: Vec<usize> = (0..4).filter(|&q| v[q] < 0.0).collect();
    let cut = |a: usize, b: usize| lerp(pts[a], pts[b], v[a] / (v[a] - v[b]));
    match pos.len() {
        0 => 0.0,
        4 => total,
        1 => {
            let a = pos[0];
            let (b, c, d) = (neg[0], neg[1], neg[2]);
            tet_volume([pts[a], cut(a, b), cut(a, c), cut(a, d)])
        }
        3 => {
            let a = neg[0];
            let (b, c, d) = (pos[0], pos[1], pos[2]);
            total - tet_volume([pts[a], cut(a, b), cut(a, c), cut(a, d)])
        }
        _ => {
            // Two fluid corners A, B; the fluid wedge {A, B, P_AC, P_AD, P_BC, P_BD}
            // as three tetrahedra.
            let (a, b) = (pos[0], pos[1]);
            let (c, d) = (neg[0], neg[1]);
            let (pac, pad, pbc, pbd) = (cut(a, c), cut(a, d), cut(b, c), cut(b, d));
            tet_volume([pts[a], pts[b], pbc, pbd])
                + tet_volume([pts[a], pac, pbc, pbd])
                + tet_volume([pts[a], pac, pad, pbd])
        }
    }
}