alice-physics 1.1.0

//! Fluid Flow Visualization
//!
//! Generates visualization data from fluid particle simulations:
//!
//! - **Flow Arrows**: Grid-sampled velocity arrows for vector field display
//! - **Streamlines**: Particle traces through the velocity field
//!
//! All computations use deterministic 128-bit fixed-point arithmetic.

use crate::math::{Fix128, Vec3Fix};

#[cfg(not(feature = "std"))]
use alloc::vec;
#[cfg(not(feature = "std"))]
use alloc::vec::Vec;

// ============================================================================
// Flow Visualization Types
// ============================================================================

/// A velocity arrow for flow visualization.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct FlowArrow {
    /// Arrow base position
    pub position: Vec3Fix,
    /// Velocity direction (normalized)
    pub direction: Vec3Fix,
    /// Velocity magnitude
    pub magnitude: Fix128,
}

/// Collection of streamlines stored as a flat buffer with per-line lengths.
///
/// `points[offsets[i]..offsets[i+1]]` gives the polyline for streamline `i`.
#[derive(Clone, Debug)]
pub struct Streamlines {
    /// All polyline points packed contiguously.
    pub points: Vec<Vec3Fix>,
    /// Start index in `points` for each line; length = num_lines + 1.
    pub offsets: Vec<usize>,
}

impl Streamlines {
    /// Number of streamlines.
    #[inline]
    #[must_use]
    pub fn len(&self) -> usize {
        self.offsets.len().saturating_sub(1)
    }

    /// Returns `true` if there are no streamlines.
    #[inline]
    #[must_use]
    pub fn is_empty(&self) -> bool {
        self.len() == 0
    }

    /// Get the polyline for streamline `i`.
    #[inline]
    #[must_use]
    pub fn line(&self, i: usize) -> &[Vec3Fix] {
        let s = self.offsets[i];
        let e = self.offsets[i + 1];
        &self.points[s..e]
    }
}

/// Configuration for flow visualization.
#[derive(Clone, Debug)]
pub struct FlowVizConfig {
    /// Grid resolution (cells per axis)
    pub grid_resolution: usize,
    /// Minimum corner of the visualization volume
    pub bounds_min: Vec3Fix,
    /// Maximum corner of the visualization volume
    pub bounds_max: Vec3Fix,
    /// Scale factor for arrow lengths
    pub arrow_scale: Fix128,
}

// ============================================================================
// Flow Arrow Generation
// ============================================================================

/// Generate flow arrows from fluid particle data.
///
/// Samples velocity on a regular 3D grid by averaging nearby fluid particle
/// velocities using a radial kernel. Each grid cell with non-zero velocity
/// produces a flow arrow.
///
/// # Arguments
///
/// - `fluid_positions`: Fluid particle positions
/// - `fluid_velocities`: Fluid particle velocities
/// - `config`: Visualization configuration
#[must_use]
pub fn generate_flow_arrows(
    fluid_positions: &[Vec3Fix],
    fluid_velocities: &[Vec3Fix],
    config: &FlowVizConfig,
) -> Vec<FlowArrow> {
    if fluid_positions.is_empty() {
        return Vec::new();
    }

    let res = config.grid_resolution.max(1);
    let min = config.bounds_min;
    let max = config.bounds_max;

    let dx = (max.x - min.x) / Fix128::from_int(res as i64);
    let dy = (max.y - min.y) / Fix128::from_int(res as i64);
    let dz = (max.z - min.z) / Fix128::from_int(res as i64);

    // Cell size determines the search radius for averaging
    let cell_size = dx;
    let radius_sq = cell_size * cell_size;

    let mut arrows = Vec::new();

    for iz in 0..res {
        for iy in 0..res {
            for ix in 0..res {
                // Cell center
                let center = Vec3Fix::new(
                    min.x + dx * Fix128::from_int(ix as i64) + dx.half(),
                    min.y + dy * Fix128::from_int(iy as i64) + dy.half(),
                    min.z + dz * Fix128::from_int(iz as i64) + dz.half(),
                );

                // Average velocity of nearby particles
                let mut avg_vel = Vec3Fix::ZERO;
                let mut count = Fix128::ZERO;

                for (pi, pos) in fluid_positions.iter().enumerate() {
                    let diff = *pos - center;
                    let dist_sq = diff.dot(diff);

                    if dist_sq < radius_sq {
                        if pi < fluid_velocities.len() {
                            avg_vel = avg_vel + fluid_velocities[pi];
                        }
                        count = count + Fix128::ONE;
                    }
                }

                if count.is_zero() {
                    continue;
                }

                avg_vel = avg_vel / count;
                let magnitude = avg_vel.length();

                if magnitude.is_zero() {
                    continue;
                }

                let direction = avg_vel / magnitude;

                arrows.push(FlowArrow {
                    position: center,
                    direction,
                    magnitude: magnitude * config.arrow_scale,
                });
            }
        }
    }

    arrows
}

// ============================================================================
// Streamline Generation
// ============================================================================

/// Generate streamlines by tracing particles through a velocity field.
///
/// Starting from each seed point, advances forward through the velocity
/// field (interpolated from nearby fluid particles) for the given number
/// of steps. Returns a [`Streamlines`] struct with all polylines in a
/// single flat buffer.
///
/// # Arguments
///
/// - `fluid_positions`: Fluid particle positions
/// - `fluid_velocities`: Fluid particle velocities
/// - `seed_points`: Starting positions for streamlines
/// - `steps`: Number of integration steps per streamline
/// - `dt`: Time step for Euler integration
#[must_use]
pub fn generate_streamlines(
    fluid_positions: &[Vec3Fix],
    fluid_velocities: &[Vec3Fix],
    seed_points: &[Vec3Fix],
    steps: usize,
    dt: Fix128,
) -> Streamlines {
    let num_seeds = seed_points.len();
    let mut points: Vec<Vec3Fix> = Vec::with_capacity(num_seeds * (steps + 1));
    let mut offsets: Vec<usize> = Vec::with_capacity(num_seeds + 1);

    if fluid_positions.is_empty() || dt.is_zero() {
        for &s in seed_points {
            offsets.push(points.len());
            points.push(s);
        }
        offsets.push(points.len());
        return Streamlines { points, offsets };
    }

    let influence_radius = Fix128::ONE;
    let radius_sq = influence_radius * influence_radius;

    for &seed in seed_points {
        offsets.push(points.len());
        let mut pos = seed;
        points.push(pos);

        for _ in 0..steps {
            let vel = interpolate_velocity(
                pos,
                fluid_positions,
                fluid_velocities,
                radius_sq,
                influence_radius,
            );

            if vel.length_squared().is_zero() {
                break;
            }

            pos = pos + vel * dt;
            points.push(pos);
        }
    }
    offsets.push(points.len());

    Streamlines { points, offsets }
}

/// Interpolate velocity at a point using inverse-distance weighting.
fn interpolate_velocity(
    pos: Vec3Fix,
    fluid_positions: &[Vec3Fix],
    fluid_velocities: &[Vec3Fix],
    radius_sq: Fix128,
    radius: Fix128,
) -> Vec3Fix {
    let mut weighted_vel = Vec3Fix::ZERO;
    let mut total_weight = Fix128::ZERO;

    for (i, fp) in fluid_positions.iter().enumerate() {
        let diff = pos - *fp;
        let dist_sq = diff.dot(diff);

        if dist_sq >= radius_sq || i >= fluid_velocities.len() {
            continue;
        }

        // Linear weight: 1 - dist/radius
        let dist = dist_sq.sqrt();
        let weight = Fix128::ONE - dist / radius;

        if weight > Fix128::ZERO {
            weighted_vel = weighted_vel + fluid_velocities[i] * weight;
            total_weight = total_weight + weight;
        }
    }

    if total_weight.is_zero() {
        Vec3Fix::ZERO
    } else {
        weighted_vel / total_weight
    }
}

// ============================================================================
// Tests
// ============================================================================

#[cfg(test)]
mod tests {
    use super::*;

    fn basic_config() -> FlowVizConfig {
        FlowVizConfig {
            grid_resolution: 2,
            bounds_min: Vec3Fix::from_int(-2, -2, -2),
            bounds_max: Vec3Fix::from_int(2, 2, 2),
            arrow_scale: Fix128::ONE,
        }
    }

    #[test]
    fn test_flow_arrows_empty_fluid() {
        let config = basic_config();
        let arrows = generate_flow_arrows(&[], &[], &config);
        assert!(arrows.is_empty());
    }

    #[test]
    fn test_flow_arrows_uniform_field() {
        let config = FlowVizConfig {
            grid_resolution: 2,
            bounds_min: Vec3Fix::from_int(-1, -1, -1),
            bounds_max: Vec3Fix::from_int(1, 1, 1),
            arrow_scale: Fix128::ONE,
        };
        // Particles everywhere moving in +X
        let mut positions = Vec::new();
        let mut velocities = Vec::new();
        for ix in -1..=1 {
            for iy in -1..=1 {
                for iz in -1..=1 {
                    positions.push(Vec3Fix::from_int(ix, iy, iz));
                    velocities.push(Vec3Fix::from_int(1, 0, 0));
                }
            }
        }

        let arrows = generate_flow_arrows(&positions, &velocities, &config);
        // Should produce some arrows
        assert!(!arrows.is_empty());

        // All arrows should point in +X direction
        for arrow in &arrows {
            assert!(arrow.direction.x > Fix128::ZERO);
            assert!(arrow.magnitude > Fix128::ZERO);
        }
    }

    #[test]
    fn test_flow_arrows_magnitude_scaled() {
        let config = FlowVizConfig {
            grid_resolution: 1,
            bounds_min: Vec3Fix::from_int(-1, -1, -1),
            bounds_max: Vec3Fix::from_int(1, 1, 1),
            arrow_scale: Fix128::from_int(2),
        };
        let positions = vec![Vec3Fix::ZERO];
        let velocities = vec![Vec3Fix::from_int(3, 0, 0)];

        let arrows = generate_flow_arrows(&positions, &velocities, &config);
        if !arrows.is_empty() {
            // Magnitude should be velocity * scale
            assert!(arrows[0].magnitude > Fix128::from_int(5));
        }
    }

    #[test]
    fn test_streamlines_empty_fluid() {
        let seeds = vec![Vec3Fix::ZERO];
        let lines = generate_streamlines(&[], &[], &seeds, 10, Fix128::from_ratio(1, 60));
        assert_eq!(lines.len(), 1);
        assert_eq!(lines.line(0).len(), 1); // Only seed point
    }

    #[test]
    fn test_streamlines_single_particle() {
        let positions = vec![Vec3Fix::ZERO];
        let velocities = vec![Vec3Fix::from_int(1, 0, 0)];
        let seeds = vec![Vec3Fix::ZERO];

        let lines = generate_streamlines(
            &positions,
            &velocities,
            &seeds,
            5,
            Fix128::from_ratio(1, 10),
        );
        assert_eq!(lines.len(), 1);
        // Should have seed + some trace points
        assert!(lines.line(0).len() > 1);
        // Should have moved in +X direction
        let last = lines.line(0).last().unwrap();
        assert!(last.x > Fix128::ZERO);
    }

    #[test]
    fn test_streamlines_zero_dt() {
        let positions = vec![Vec3Fix::ZERO];
        let velocities = vec![Vec3Fix::from_int(1, 0, 0)];
        let seeds = vec![Vec3Fix::ZERO];

        let lines = generate_streamlines(&positions, &velocities, &seeds, 10, Fix128::ZERO);
        assert_eq!(lines.line(0).len(), 1); // Only seed
    }

    #[test]
    fn test_streamlines_multiple_seeds() {
        let positions = vec![Vec3Fix::ZERO, Vec3Fix::from_int(1, 0, 0)];
        let velocities = vec![Vec3Fix::UNIT_X, Vec3Fix::UNIT_X];
        let seeds = vec![
            Vec3Fix::from_int(-1, 0, 0),
            Vec3Fix::from_int(0, 1, 0),
            Vec3Fix::from_int(0, 0, -1),
        ];

        let lines = generate_streamlines(
            &positions,
            &velocities,
            &seeds,
            3,
            Fix128::from_ratio(1, 10),
        );
        assert_eq!(lines.len(), 3);
    }

    #[test]
    fn test_interpolate_velocity_no_neighbors() {
        let vel = interpolate_velocity(
            Vec3Fix::from_int(100, 100, 100),
            &[Vec3Fix::ZERO],
            &[Vec3Fix::from_int(1, 0, 0)],
            Fix128::ONE,
            Fix128::ONE,
        );
        assert!(vel.length_squared().is_zero());
    }

    #[test]
    fn test_flow_arrow_direction_normalized() {
        let config = FlowVizConfig {
            grid_resolution: 1,
            bounds_min: Vec3Fix::from_int(-1, -1, -1),
            bounds_max: Vec3Fix::from_int(1, 1, 1),
            arrow_scale: Fix128::ONE,
        };
        let positions = vec![Vec3Fix::ZERO];
        let velocities = vec![Vec3Fix::from_int(3, 4, 0)];

        let arrows = generate_flow_arrows(&positions, &velocities, &config);
        if !arrows.is_empty() {
            let len = arrows[0].direction.length();
            let eps = Fix128::from_ratio(1, 100);
            let diff = (len - Fix128::ONE).abs();
            assert!(diff < eps, "Direction should be normalized");
        }
    }
}