//! The binary ghost wall — the 2D `EmbeddedMask` on the 3D grid. Cells are
//! fluid where φ(centre) > 0; an interior face between two cells that are
//! not both fluid is a ghost (within 1.5 h of the surface) or solid; a
//! ghost face is prescribed the least-squares linear reconstruction
//! through the fluid nodes of the trilinear stencil around the probe (one
//! cell beyond the face's mirror image along the normal) plus the foot of
//! the normal with its surface velocity — exact for linear fields — with
//! the profile along the normal as the fallback; the ghost faces bounding
//! fluid cells share a flux compatibility correction.

use super::Grid;
use super::body::Body;
use super::cut::CutGeometry;
use super::step::{Boundaries, Side};

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum FaceKind {
    Fluid,
    Ghost,
    Solid,
}

/// The wall treatment of an embedded body.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum WallScheme {
    /// The 2D binary mask ported (this module).
    #[default]
    GhostBinary,
    /// The apertured cut-cell wall (item 10).
    CutCell,
}

#[derive(Debug, Clone, Copy)]
pub(crate) struct StencilNode {
    pub(crate) idx: usize,
    pub(crate) x: f64,
    pub(crate) y: f64,
    pub(crate) z: f64,
    /// The trilinear weight (the fallback profile's interpolation).
    pub(crate) weight: f64,
    /// The z-direction weight alone: the least-squares weight of the node,
    /// so that a z-invariant field fits exactly as the 2D four-node fit
    /// (the two planes share the in-plane nodes' unit weight).
    pub(crate) zw: f64,
    pub(crate) fallback: Option<f64>,
}

#[derive(Debug, Clone)]
pub(super) struct Ghost {
    pub(super) idx: usize,
    x: f64,
    y: f64,
    z: f64,
    foot: (f64, f64, f64),
    u_surface: f64,
    s_face: f64,
    s_probe: f64,
    nodes: Vec<StencilNode>,
    /// Outward-from-fluid sign for the compatibility correction (0 when no
    /// fluid cell is adjacent).
    pub(super) flux_sign: f64,
}

#[derive(Clone)]
pub struct Mask {
    pub(super) grid: Grid,
    pub(super) periodic_z: bool,
    pub(super) cell_fluid: Vec<bool>,
    pub(super) u_kind: Vec<FaceKind>,
    pub(super) v_kind: Vec<FaceKind>,
    pub(super) w_kind: Vec<FaceKind>,
    pub(super) u_ghosts: Vec<Ghost>,
    pub(super) v_ghosts: Vec<Ghost>,
    pub(super) w_ghosts: Vec<Ghost>,
    pub(super) anchor: usize,
    pub(super) fluid_cells: usize,
    /// The cut geometry of the apertured wall (`WallScheme::CutCell`,
    /// `cutwall.rs`); `None` on the binary ghost wall.
    pub(super) cut: Option<CutGeometry>,
    /// The step-averaged apertures `½(αⁿ + αⁿ⁺¹)` of a moving cut wall
    /// (the space-time continuity: a cell's volume change over the step
    /// equals the flux through the apertures it had during it); `None` =
    /// the instantaneous ones.
    pub(super) step_apertures: Option<(Vec<f64>, Vec<f64>, Vec<f64>)>,
    /// The projection's space-time classification on a moving cut wall:
    /// a face is an unknown where its step-averaged aperture is positive,
    /// a cell has an equation where it holds fluid at either end of the
    /// step (a dying cell empties through the apertures it had); `None` =
    /// the instantaneous kinds.
    pub(super) step_open: Option<(Vec<bool>, Vec<bool>, Vec<bool>, Vec<bool>)>,
    /// Virtual merging (item 10b): the master cell of every small cell
    /// (`usize::MAX` = its own row) — a small cell shares its pressure
    /// unknown with its largest active face neighbour in the projection.
    pub(super) merge_master: Vec<usize>,
    /// The predictor's convection scheme (the exchange route replicates
    /// its limited fluxes on the faces next to prescribed ones).
    pub(super) scheme: crate::solvers::incompressible::ConvectionScheme,
    /// The fluid's density (the exchange route's convective flux).
    pub(super) density: f64,
    /// The cut wall's shear closure order (S2-4).
    pub(super) wall_order: u8,
    /// The oblique wall distance of the cut faces (S2-5).
    pub(super) wall_distance_oblique: bool,
    /// The transverse centroid correction and the fine distance floor (S2-6).
    pub(super) diffusion_transverse: bool,
    pub(super) distance_floor_fine: bool,
    pub(super) wall_advancing: bool,
    pub(super) exchange_convection_off: bool,
    /// The axis-distance implicit wall exchange (S2-5).
    pub(super) wall_exchange_axis: bool,
    /// S2-7: the momentum control volumes' side apertures from the
    /// interpolant on the sides' own corners (host prototype).
    pub(super) cv_sides_exact: bool,
    /// S2-7: the quadratic wall gradient's second point at the neighbour's
    /// own centroid distance (host prototype).
    pub(super) wall_order2_centroid: bool,
    /// S2-7b host prototypes: the axis foot's wall velocity in the solid
    /// exchange; the convective sides from the sides' own apertures.
    pub(super) wall_exchange_foot: bool,
    pub(super) conv_sides_exact: bool,
    /// The centroid prototype's pressure-gradient weights per u / v / w face.
    /// The centroid-distance cross diffusion (S2-5).
    pub(super) diffusion_centroid: bool,
    /// The open-part centroid shifts per u / v / w face, three components
    /// interleaved (built with `diffusion_centroid`).
    pub(super) face_shifts: Option<[Vec<f64>; 3]>,
    pub(super) grad_weights: Option<(Vec<f64>, Vec<f64>, Vec<f64>)>,
}

/// The z lattice position of a query: the lower plane index, the upper
/// plane index and the weight of the upper plane. Periodic z wraps; a wall
/// clamps (and a query outside the lattice takes the nearest plane).
pub(crate) fn z_planes(gz: f64, planes: usize, periodic: bool) -> (usize, usize, f64) {
    if planes <= 1 {
        return (0, 0, 0.0);
    }
    if periodic {
        let n = planes as f64;
        let w = gz.rem_euclid(n);
        let k0 = w.floor() as usize % planes;
        ((k0) % planes, (k0 + 1) % planes, w - w.floor())
    } else {
        let k0 = gz.floor().clamp(0.0, (planes - 2) as f64) as usize;
        (k0, k0 + 1, (gz - k0 as f64).clamp(0.0, 1.0))
    }
}

/// The nodes of velocity component `c` (0 u, 1 v, 2 w) around `point`
/// with trilinear weights (a single plane of nodes in a direction
/// collapses to one node), clamped in x and y, wrapped or clamped in z;
/// `fallback(idx)` supplies the value of a node that is not a fluid face.
pub(crate) fn stencil_nodes(
    point: (f64, f64, f64),
    c: usize,
    g: Grid,
    periodic_z: bool,
    fallback: impl Fn(usize) -> Option<f64>,
) -> Vec<StencilNode> {
    let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
    // w faces on a periodic grid: the plane k = nz is the plane 0.
    let (gx, gy, gz, max_i, max_j, planes_k) = match c {
        0 => (
            point.0 / dx,
            point.1 / dy - 0.5,
            point.2 / dz - 0.5,
            nx,
            ny - 1,
            nz,
        ),
        1 => (
            point.0 / dx - 0.5,
            point.1 / dy,
            point.2 / dz - 0.5,
            nx - 1,
            ny,
            nz,
        ),
        _ => (
            point.0 / dx - 0.5,
            point.1 / dy - 0.5,
            point.2 / dz,
            nx - 1,
            ny - 1,
            if periodic_z { nz } else { nz + 1 },
        ),
    };
    let i0 = gx.floor().clamp(0.0, (max_i.max(1) - 1) as f64) as usize;
    let j0 = gy.floor().clamp(0.0, (max_j.max(1) - 1) as f64) as usize;
    let fx = (gx - i0 as f64).clamp(0.0, 1.0);
    let fy = (gy - j0 as f64).clamp(0.0, 1.0);
    let (k0, k1, fz) = z_planes(gz, planes_k, periodic_z);
    let pos = |k: usize, j: usize, i: usize| match c {
        0 => (i as f64 * dx, (j as f64 + 0.5) * dy, (k as f64 + 0.5) * dz),
        1 => ((i as f64 + 0.5) * dx, j as f64 * dy, (k as f64 + 0.5) * dz),
        _ => ((i as f64 + 0.5) * dx, (j as f64 + 0.5) * dy, k as f64 * dz),
    };
    let index = |k: usize, j: usize, i: usize| match c {
        0 => g.uface(k, j, i),
        1 => g.vface(k, j, i),
        _ => g.wface(k, j, i),
    };
    let dirs = |m: usize, f: f64| {
        if m <= 1 {
            vec![(0usize, 1.0)]
        } else {
            vec![(0, 1.0 - f), (1, f)]
        }
    };
    let mut out = Vec::with_capacity(8);
    let kz = |dk: usize| if dk == 0 { k0 } else { k1 };
    for (dk, wk) in dirs(planes_k, fz) {
        for (dj, wj) in dirs(max_j, fy) {
            for (di, wi) in dirs(max_i, fx) {
                let (k, j, i) = (kz(dk), j0 + dj, i0 + di);
                let (x, y, z) = pos(k, j, i);
                let idx = index(k, j, i);
                out.push(StencilNode {
                    idx,
                    x,
                    y,
                    z,
                    weight: wi * wj * wk,
                    zw: wk,
                    fallback: fallback(idx),
                });
            }
        }
    }
    out
}

/// Weighted least-squares fit of `a + b (x−x0) + c (y−y0) + d (z−z0)`
/// through `pts` (each `(x, y, z, value, weight)`; zero-weight points drop
/// out) evaluated at `at`; directions in which every weighted point has
/// the same coordinate are dropped. `None` if the normal matrix is
/// singular relative to its scale.
pub(crate) fn linear_fit(pts_w: &[(f64, f64, f64, f64, f64)], at: (f64, f64, f64)) -> Option<f64> {
    let pts: Vec<(f64, f64, f64, f64, f64)> = pts_w.iter().copied().filter(|p| p.4 > 0.0).collect();
    let spread = |f: &dyn Fn(&(f64, f64, f64, f64, f64)) -> f64| {
        let (lo, hi) = pts
            .iter()
            .fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), p| {
                (lo.min(f(p)), hi.max(f(p)))
            });
        hi - lo > 1e-12 * (hi.abs() + lo.abs() + 1e-300)
    };
    let on = [true, spread(&|p| p.0), spread(&|p| p.1), spread(&|p| p.2)];
    let cols: Vec<usize> = (0..4).filter(|&c| on[c]).collect();
    let m = cols.len();
    if pts.len() < m {
        return None;
    }
    let mut a = vec![vec![0.0; m + 1]; m];
    for p in &pts {
        let full = [1.0, p.0 - at.0, p.1 - at.1, p.2 - at.2];
        let r: Vec<f64> = cols.iter().map(|&c| full[c]).collect();
        for i in 0..m {
            for j in 0..m {
                a[i][j] += p.4 * r[i] * r[j];
            }
            a[i][m] += p.4 * r[i] * p.3;
        }
    }
    let scale: f64 = (0..m).map(|i| a[i][i]).product();
    if scale <= 0.0 {
        return None;
    }
    let mut det = 1.0;
    for col in 0..m {
        let piv = (col..m)
            .max_by(|&p, &q| a[p][col].abs().partial_cmp(&a[q][col].abs()).unwrap())
            .unwrap();
        a.swap(col, piv);
        let d = a[col][col];
        det *= d;
        if d.abs() <= 1e-14 * scale.powf(1.0 / m as f64) {
            return None;
        }
        for r in col + 1..m {
            let f = a[r][col] / d;
            for c in col..=m {
                a[r][c] -= f * a[col][c];
            }
        }
    }
    if det.abs() <= 1e-10 * scale {
        return None;
    }
    let mut x = vec![0.0; m];
    for i in (0..m).rev() {
        let mut s = a[i][m];
        for j in i + 1..m {
            s -= a[i][j] * x[j];
        }
        x[i] = s / a[i][i];
    }
    Some(x[0])
}

impl Ghost {
    pub(super) fn reconstruct(&self, values: &[f64]) -> f64 {
        let mut pts: Vec<(f64, f64, f64, f64, f64)> = self
            .nodes
            .iter()
            .filter(|n| n.fallback.is_none())
            .map(|n| (n.x, n.y, n.z, values[n.idx], n.zw))
            .collect();
        pts.push((self.foot.0, self.foot.1, self.foot.2, self.u_surface, 1.0));
        if let Some(val) = linear_fit(&pts, (self.x, self.y, self.z)) {
            return val;
        }
        let mut probe = 0.0;
        for n in &self.nodes {
            probe += n.weight * n.fallback.unwrap_or_else(|| values[n.idx]);
        }
        self.u_surface + (probe - self.u_surface) * (self.s_face / self.s_probe)
    }
}

impl Mask {
    /// Classify the grid against `body` at `t`. A solid cell on a domain
    /// side is refused unless that side is a Velocity, SlipWall or Periodic
    /// side (an outlet may not be blocked).
    pub fn build(body: &Body, g: Grid, t: f64, b: Boundaries) -> Result<Self, String> {
        let (nx, ny, nz, dx, dy, dz) = (g.nx, g.ny, g.nz, g.dx, g.dy, g.dz);
        let periodic = b.z0 == Side::Periodic;
        let xc = |i: usize| (i as f64 + 0.5) * dx;
        let yc = |j: usize| (j as f64 + 0.5) * dy;
        let zc = |k: usize| (k as f64 + 0.5) * dz;
        let mut cell_fluid = vec![true; g.cells()];
        let mut fluid_cells = 0;
        let mut anchor = None;
        let allowed = |side: Side| matches!(side, Side::Velocity | Side::Periodic | Side::SlipWall);
        for k in 0..nz {
            for j in 0..ny {
                for i in 0..nx {
                    let fluid = body.phi(xc(i), yc(j), zc(k), t) > 0.0;
                    let idx = g.cell(k, j, i);
                    cell_fluid[idx] = fluid;
                    if fluid {
                        fluid_cells += 1;
                        if anchor.is_none() {
                            anchor = Some(idx);
                        }
                    } else {
                        let touches = (i == 0 && !allowed(b.x0))
                            || (i + 1 == nx && !allowed(b.x1))
                            || (j == 0 && !allowed(b.y0))
                            || (j + 1 == ny && !allowed(b.y1))
                            || (k == 0 && !allowed(b.z0))
                            || (k + 1 == nz && !allowed(b.z1));
                        if touches {
                            return Err(format!(
                                "embedded body reaches a domain side that is not a Velocity/Periodic side at cell ({k}, {j}, {i})"
                            ));
                        }
                    }
                }
            }
        }
        let Some(anchor) = anchor else {
            return Err("embedded body covers the whole domain".into());
        };
        let h_min = dx.min(dy).min(dz);
        let reach = 1.5 * h_min;
        let eps = 1e-6 * h_min;
        let is_fluid = |k: usize, j: usize, i: usize| cell_fluid[g.cell(k, j, i)];
        let mut u_kind = vec![FaceKind::Fluid; g.n_ufaces()];
        let mut v_kind = vec![FaceKind::Fluid; g.n_vfaces()];
        let mut w_kind = vec![FaceKind::Fluid; g.n_wfaces()];
        let kind_of = |phi: f64| {
            if phi > -reach {
                FaceKind::Ghost
            } else {
                FaceKind::Solid
            }
        };
        for k in 0..nz {
            for j in 0..ny {
                for i in 1..nx {
                    if !(is_fluid(k, j, i - 1) && is_fluid(k, j, i)) {
                        u_kind[g.uface(k, j, i)] =
                            kind_of(body.phi(i as f64 * dx, yc(j), zc(k), t));
                    }
                }
            }
            for j in 1..ny {
                for i in 0..nx {
                    if !(is_fluid(k, j - 1, i) && is_fluid(k, j, i)) {
                        v_kind[g.vface(k, j, i)] =
                            kind_of(body.phi(xc(i), j as f64 * dy, zc(k), t));
                    }
                }
            }
        }
        let w_range = if periodic { 0..nz } else { 1..nz };
        for k in w_range.clone() {
            let below = if k > 0 { k - 1 } else { nz - 1 };
            for j in 0..ny {
                for i in 0..nx {
                    if !(is_fluid(below, j, i) && is_fluid(k, j, i)) {
                        w_kind[g.wface(k, j, i)] =
                            kind_of(body.phi(xc(i), yc(j), k as f64 * dz, t));
                    }
                }
            }
        }
        if periodic {
            for j in 0..ny {
                for i in 0..nx {
                    w_kind[g.wface(nz, j, i)] = w_kind[g.wface(0, j, i)];
                }
            }
        }
        let face_pos = |c: usize, k: usize, j: usize, i: usize| match c {
            0 => (i as f64 * dx, yc(j), zc(k)),
            1 => (xc(i), j as f64 * dy, zc(k)),
            _ => (xc(i), yc(j), k as f64 * dz),
        };
        let kji_of = |c: usize, idx: usize| -> (usize, usize, usize) {
            match c {
                0 => (idx / ((nx + 1) * ny), (idx / (nx + 1)) % ny, idx % (nx + 1)),
                1 => (idx / (nx * (ny + 1)), (idx / nx) % (ny + 1), idx % nx),
                _ => (idx / (nx * ny), (idx / nx) % ny, idx % nx),
            }
        };
        let build_ghost = |c: usize,
                           kinds: &[FaceKind],
                           idx: usize,
                           (x, y, z): (f64, f64, f64),
                           flux_sign: f64|
         -> Ghost {
            let s_face = body.phi(x, y, z, t);
            let (n1, n2, n3) = body.normal(x, y, z, t, eps);
            let foot = (x - s_face * n1, y - s_face * n2, z - s_face * n3);
            let s_probe = s_face.abs() + h_min;
            let probe = (
                foot.0 + s_probe * n1,
                foot.1 + s_probe * n2,
                foot.2 + s_probe * n3,
            );
            let vel = body.surface_velocity(foot.0, foot.1, foot.2, t);
            let u_surface = [vel.0, vel.1, vel.2][c];
            let nodes = stencil_nodes(probe, c, g, periodic, |nidx| {
                if kinds[nidx] == FaceKind::Fluid {
                    None
                } else {
                    let (kk, jj, ii) = kji_of(c, nidx);
                    let (px, py, pz) = face_pos(c, kk, jj, ii);
                    let s = body.phi(px, py, pz, t);
                    let (m1, m2, m3) = body.normal(px, py, pz, t, eps);
                    let f = body.surface_velocity(px - s * m1, py - s * m2, pz - s * m3, t);
                    Some([f.0, f.1, f.2][c])
                }
            });
            Ghost {
                idx,
                x,
                y,
                z,
                foot,
                u_surface,
                s_face,
                s_probe,
                nodes,
                flux_sign,
            }
        };
        let mut u_ghosts = Vec::new();
        let mut v_ghosts = Vec::new();
        let mut w_ghosts = Vec::new();
        for k in 0..nz {
            for j in 0..ny {
                for i in 1..nx {
                    let idx = g.uface(k, j, i);
                    if u_kind[idx] == FaceKind::Ghost {
                        let sign = if is_fluid(k, j, i - 1) {
                            1.0
                        } else if is_fluid(k, j, i) {
                            -1.0
                        } else {
                            0.0
                        };
                        u_ghosts.push(build_ghost(0, &u_kind, idx, face_pos(0, k, j, i), sign));
                    }
                }
            }
            for j in 1..ny {
                for i in 0..nx {
                    let idx = g.vface(k, j, i);
                    if v_kind[idx] == FaceKind::Ghost {
                        let sign = if is_fluid(k, j - 1, i) {
                            1.0
                        } else if is_fluid(k, j, i) {
                            -1.0
                        } else {
                            0.0
                        };
                        v_ghosts.push(build_ghost(1, &v_kind, idx, face_pos(1, k, j, i), sign));
                    }
                }
            }
        }
        for k in w_range {
            let below = if k > 0 { k - 1 } else { nz - 1 };
            for j in 0..ny {
                for i in 0..nx {
                    let idx = g.wface(k, j, i);
                    if w_kind[idx] == FaceKind::Ghost {
                        let sign = if is_fluid(below, j, i) {
                            1.0
                        } else if is_fluid(k, j, i) {
                            -1.0
                        } else {
                            0.0
                        };
                        w_ghosts.push(build_ghost(2, &w_kind, idx, face_pos(2, k, j, i), sign));
                    }
                }
            }
        }
        Ok(Self {
            grid: g,
            periodic_z: periodic,
            cell_fluid,
            u_kind,
            v_kind,
            w_kind,
            u_ghosts,
            v_ghosts,
            w_ghosts,
            anchor,
            fluid_cells,
            cut: None,
            step_apertures: None,
            step_open: None,
            merge_master: Vec::new(),
            scheme: crate::solvers::incompressible::ConvectionScheme::Upwind,
            density: 1.0,
            wall_order: 1,
            wall_distance_oblique: false,
            diffusion_transverse: false,
            distance_floor_fine: false,
            wall_advancing: false,
            exchange_convection_off: false,
            wall_exchange_axis: false,
            cv_sides_exact: false,
            wall_order2_centroid: false,
            wall_exchange_foot: false,
            conv_sides_exact: false,
            grad_weights: None,
            diffusion_centroid: false,
            face_shifts: None,
        })
    }

    /// The step-averaged apertures of a moving cut wall (`None` at rest).
    #[must_use]
    pub fn step_apertures(&self) -> Option<&(Vec<f64>, Vec<f64>, Vec<f64>)> {
        self.step_apertures.as_ref()
    }

    /// The master of a virtually merged small cell.
    #[inline]
    #[must_use]
    pub fn master(&self, idx: usize) -> Option<usize> {
        match self.merge_master.get(idx) {
            Some(&m) if m != usize::MAX => Some(m),
            _ => None,
        }
    }
    #[must_use]
    pub fn merged_cells(&self) -> usize {
        self.merge_master
            .iter()
            .filter(|&&m| m != usize::MAX)
            .count()
    }

    // The projection's unknowns (the instantaneous kinds at rest).
    #[inline]
    #[must_use]
    pub fn u_open(&self, idx: usize) -> bool {
        self.step_open
            .as_ref()
            .map_or(self.u_kind[idx] == FaceKind::Fluid, |o| o.0[idx])
    }
    #[inline]
    #[must_use]
    pub fn v_open(&self, idx: usize) -> bool {
        self.step_open
            .as_ref()
            .map_or(self.v_kind[idx] == FaceKind::Fluid, |o| o.1[idx])
    }
    #[inline]
    #[must_use]
    pub fn w_open(&self, idx: usize) -> bool {
        self.step_open
            .as_ref()
            .map_or(self.w_kind[idx] == FaceKind::Fluid, |o| o.2[idx])
    }
    #[inline]
    #[must_use]
    pub fn cell_active(&self, idx: usize) -> bool {
        self.step_open
            .as_ref()
            .map_or(self.cell_fluid[idx], |o| o.3[idx])
    }

    /// The step-averaged aperture of a u / v / w face (the instantaneous
    /// one for a wall at rest).
    #[inline]
    #[must_use]
    pub fn au_step(&self, idx: usize) -> f64 {
        self.step_apertures
            .as_ref()
            .map_or_else(|| self.a_u(idx), |a| a.0[idx])
    }
    #[inline]
    #[must_use]
    pub fn av_step(&self, idx: usize) -> f64 {
        self.step_apertures
            .as_ref()
            .map_or_else(|| self.a_v(idx), |a| a.1[idx])
    }
    #[inline]
    #[must_use]
    pub fn aw_step(&self, idx: usize) -> f64 {
        self.step_apertures
            .as_ref()
            .map_or_else(|| self.a_w(idx), |a| a.2[idx])
    }

    /// The cut geometry (apertured wall only).
    #[must_use]
    pub fn cut(&self) -> Option<&CutGeometry> {
        self.cut.as_ref()
    }
    /// Fluid area fraction of a u / v / w face (1 on the binary wall).
    #[inline]
    #[must_use]
    pub fn a_u(&self, idx: usize) -> f64 {
        self.cut.as_ref().map_or(1.0, |c| c.a_u[idx])
    }
    #[inline]
    #[must_use]
    pub fn a_v(&self, idx: usize) -> f64 {
        self.cut.as_ref().map_or(1.0, |c| c.a_v[idx])
    }
    #[inline]
    #[must_use]
    pub fn a_w(&self, idx: usize) -> f64 {
        self.cut.as_ref().map_or(1.0, |c| c.a_w[idx])
    }
    /// Fluid volume fraction of a cell (1 on the binary wall).
    #[inline]
    #[must_use]
    pub fn grad_weight(&self, c: usize, idx: usize) -> f64 {
        self.grad_weights.as_ref().map_or(1.0, |w| match c {
            0 => w.0[idx],
            1 => w.1[idx],
            _ => w.2[idx],
        })
    }
    #[inline]
    #[must_use]
    pub fn vol(&self, idx: usize) -> f64 {
        self.cut.as_ref().map_or(1.0, |c| c.vol[idx])
    }

    #[inline]
    #[must_use]
    pub fn is_fluid_cell(&self, idx: usize) -> bool {
        self.cell_fluid[idx]
    }
    #[inline]
    #[must_use]
    pub fn u_kind(&self, idx: usize) -> FaceKind {
        self.u_kind[idx]
    }
    #[inline]
    #[must_use]
    pub fn v_kind(&self, idx: usize) -> FaceKind {
        self.v_kind[idx]
    }
    #[inline]
    #[must_use]
    pub fn w_kind(&self, idx: usize) -> FaceKind {
        self.w_kind[idx]
    }
    #[must_use]
    pub fn anchor(&self) -> usize {
        self.anchor
    }
    #[must_use]
    pub fn fluid_cells(&self) -> usize {
        self.fluid_cells
    }
    #[must_use]
    pub fn ghost_faces(&self) -> usize {
        self.u_ghosts.len() + self.v_ghosts.len() + self.w_ghosts.len()
    }
    #[must_use]
    pub fn grid(&self) -> Grid {
        self.grid
    }
    #[must_use]
    pub fn periodic_z(&self) -> bool {
        self.periodic_z
    }
}