vox_geometry_rust/

particle_system_solver3.rs

/*
 * // Copyright (c) 2021 Feng Yang
 * //
 * // I am making my contributions/submissions to this project solely in my
 * // personal capacity and am not conveying any rights to any intellectual
 * // property of any third parties.
 */

use crate::vector3::Vector3D;
use crate::particle_system_data3::*;
use crate::collider3::Collider3Ptr;
use crate::particle_emitter3::ParticleEmitter3Ptr;
use crate::vector_field3::VectorField3Ptr;
use crate::constant_vector_field3::ConstantVectorField3;
use crate::animation::{Animation, Frame};
use crate::physics_animation::{PhysicsAnimation, PhysicsAnimationData};
use log::info;
use std::time::SystemTime;
use std::sync::{RwLock, Arc};
use rayon::prelude::*;

///
/// # Basic 3-D particle system solver.
///
/// This class implements basic particle system solver. It includes gravity,
/// air drag, and collision. But it does not compute particle-to-particle
/// interaction. Thus, this solver is suitable for performing simple spray-like
/// simulations with low computational cost. This class can be further extend
/// to add more sophisticated simulations, such as SPH, to handle
/// particle-to-particle intersection.
///
/// \see        SphSolver3
///
pub struct ParticleSystemSolver3 {
    pub _drag_coefficient: f64,
    pub _restitution_coefficient: f64,
    pub _gravity: Vector3D,

    pub _particle_system_data: ParticleSystemData3Ptr,
    pub _new_positions: Vec<Vector3D>,
    pub _new_velocities: Vec<Vector3D>,
    pub _collider: Option<Collider3Ptr>,
    pub _emitter: Option<ParticleEmitter3Ptr>,
    pub _wind: VectorField3Ptr,

    pub _animation_data: PhysicsAnimationData,
}

impl ParticleSystemSolver3 {
    /// Constructs an empty solver.
    pub fn new_default() -> ParticleSystemSolver3 {
        return ParticleSystemSolver3::new(1.0e-3, 1.0e-3);
    }

    /// Constructs a solver with particle parameters.
    pub fn new(radius: f64,
               mass: f64) -> ParticleSystemSolver3 {
        let data = Arc::new(RwLock::new(ParticleSystemData3::new_default()));
        data.write().unwrap().set_radius(radius);
        data.write().unwrap().set_mass(mass);
        let wind = Arc::new(RwLock::new(ConstantVectorField3::new(Some(Vector3D::new_default()))));
        return ParticleSystemSolver3 {
            _drag_coefficient: 1.0e-4,
            _restitution_coefficient: 0.0,
            _gravity: Vector3D::new(0.0, crate::constants::K_GRAVITY, 0.0),
            _particle_system_data: data,
            _new_positions: vec![],
            _new_velocities: vec![],
            _collider: None,
            _emitter: None,
            _wind: wind,
            _animation_data: PhysicsAnimationData::new(),
        };
    }

    /// Returns builder fox ParticleSystemSolver3.
    pub fn builder() -> Builder {
        return Builder::new();
    }
}

impl ParticleSystemSolver3 {
    /// Returns the drag coefficient.
    pub fn drag_coefficient(&self) -> f64 {
        return self._drag_coefficient;
    }

    ///
    /// \brief      Sets the drag coefficient.
    ///
    /// The drag coefficient controls the amount of air-drag. The coefficient
    /// should be a positive number and 0 means no drag force.
    ///
    /// - parameter:  new_drag_coefficient The new drag coefficient.
    ///
    pub fn set_drag_coefficient(&mut self, new_drag_coefficient: f64) {
        self._drag_coefficient = f64::max(new_drag_coefficient, 0.0);
    }

    /// Sets the restitution coefficient.
    pub fn restitution_coefficient(&self) -> f64 {
        return self._restitution_coefficient;
    }

    ///
    /// \brief      Sets the restitution coefficient.
    ///
    /// The restitution coefficient controls the bouncy-ness of a particle when
    /// it hits a collider surface. The range of the coefficient should be 0 to
    /// 1 -- 0 means no bounce back and 1 means perfect reflection.
    ///
    /// - parameter:  new_restitution_coefficient The new restitution coefficient.
    ///
    pub fn set_restitution_coefficient(&mut self, new_restitution_coefficient: f64) {
        self._restitution_coefficient = crate::math_utils::clamp(new_restitution_coefficient, 0.0, 1.0);
    }

    /// Returns the gravity.
    pub fn gravity(&self) -> &Vector3D {
        return &self._gravity;
    }

    /// Sets the gravity.
    pub fn set_gravity(&mut self, new_gravity: &Vector3D) {
        self._gravity = *new_gravity;
    }

    ///
    /// \brief      Returns the particle system data.
    ///
    /// This function returns the particle system data. The data is created when
    /// this solver is constructed and also owned by the solver.
    ///
    /// \return     The particle system data.
    ///
    pub fn particle_system_data(&self) -> &ParticleSystemData3Ptr {
        return &self._particle_system_data;
    }

    /// Returns the collider.
    pub fn collider(&self) -> &Option<Collider3Ptr> {
        return &self._collider;
    }

    /// Sets the collider.
    pub fn set_collider(&mut self, new_collider: &Collider3Ptr) {
        self._collider = Some(new_collider.clone());
    }

    /// Returns the emitter.
    pub fn emitter(&self) -> &Option<ParticleEmitter3Ptr> {
        return &self._emitter;
    }

    /// Sets the emitter.
    pub fn set_emitter(&mut self, new_emitter: &ParticleEmitter3Ptr) {
        self._emitter = Some(new_emitter.clone());
        new_emitter.write().unwrap().set_target(self._particle_system_data.clone());
    }

    /// Returns the wind field.
    pub fn wind(&self) -> &VectorField3Ptr {
        return &self._wind;
    }

    ///
    /// \brief      Sets the wind.
    ///
    /// Wind can be applied to the particle system by setting a vector field to
    /// the solver.
    ///
    /// - parameter:  new_wind The new wind.
    ///
    pub fn set_wind(&mut self, new_wind: &VectorField3Ptr) {
        self._wind = new_wind.clone();
    }
}

impl Animation for ParticleSystemSolver3 {
    fn on_update(&mut self, frame: &Frame) {
        PhysicsAnimation::on_update(self, frame);
    }
}

impl PhysicsAnimation for ParticleSystemSolver3 {
    fn on_advance_time_step(&mut self, time_step_in_seconds: f64) {
        self.begin_advance_time_step(time_step_in_seconds);

        let mut timer = SystemTime::now();
        self.accumulate_forces(time_step_in_seconds);
        info!("Accumulating forces took {} seconds", timer.elapsed().unwrap().as_secs_f64());

        timer = SystemTime::now();
        self.time_integration(time_step_in_seconds);
        info!("Time integration took {} seconds", timer.elapsed().unwrap().as_secs_f64());

        timer = SystemTime::now();
        self.resolve_collision();
        info!("Resolving collision took {} seconds", timer.elapsed().unwrap().as_secs_f64());

        self.end_advance_time_step(time_step_in_seconds);
    }

    fn on_initialize(&mut self) {
        // When initializing the solver, update the collider and emitter state as
        // well since they also affects the initial condition of the simulation.
        let mut timer = SystemTime::now();
        self.update_collider(0.0);
        info!("Update collider took {} seconds", timer.elapsed().unwrap().as_secs_f64());

        timer = SystemTime::now();
        self.update_emitter(0.0);
        info!("Update emitter took {} seconds", timer.elapsed().unwrap().as_secs_f64());
    }

    fn view(&self) -> &PhysicsAnimationData {
        return &self._animation_data;
    }

    fn view_mut(&mut self) -> &mut PhysicsAnimationData {
        return &mut self._animation_data;
    }
}

impl ParticleSystemSolver3 {
    /// Accumulates forces applied to the particles.
    fn accumulate_forces(&mut self, _: f64) {
        // Add external forces
        self.accumulate_external_forces();
    }

    /// Called when a time-step is about to begin.
    fn on_begin_advance_time_step(&mut self, _: f64) {}

    /// Called after a time-step is completed.
    fn on_end_advance_time_step(&mut self, _: f64) {}

    /// Resolves any collisions occurred by the particles.
    fn resolve_collision(&mut self) {
        if let Some(collider) = &self._collider {
            let radius = self._particle_system_data.read().unwrap().radius();
            let restitution = self._restitution_coefficient;
            (&mut self._new_positions, &mut self._new_velocities).into_par_iter().for_each(|(pos, vel)| {
                collider.read().unwrap().resolve_collision(
                    radius,
                    restitution,
                    pos,
                    vel);
            });
        }
    }

    fn begin_advance_time_step(&mut self, time_step_in_seconds: f64) {
        // Clear forces
        self._particle_system_data.write().unwrap().forces_mut().fill(Vector3D::new_default());

        // Update collider and emitter
        let mut timer = SystemTime::now();
        self.update_collider(time_step_in_seconds);
        info!("Update collider took {} seconds", timer.elapsed().unwrap().as_secs_f64());

        timer = SystemTime::now();
        self.update_emitter(time_step_in_seconds);
        info!("Update emitter took {} seconds", timer.elapsed().unwrap().as_secs_f64());

        // Allocate buffers
        let n = self._particle_system_data.read().unwrap().number_of_particles();
        self._new_positions.resize(n, Vector3D::new_default());
        self._new_velocities.resize(n, Vector3D::new_default());

        self.on_begin_advance_time_step(time_step_in_seconds);
    }

    fn end_advance_time_step(&mut self, time_step_in_seconds: f64) {
        // Update data
        (self._particle_system_data.write().unwrap().positions_mut(),
         self._particle_system_data.write().unwrap().velocities_mut(),
         &self._new_positions, &self._new_velocities).into_par_iter().for_each(|(pos, vel, pos_new, vel_new)| {
            *pos = *pos_new;
            *vel = *vel_new;
        });

        self.on_end_advance_time_step(time_step_in_seconds);
    }

    fn accumulate_external_forces(&mut self) {
        let mass = self._particle_system_data.read().unwrap().mass();

        (self._particle_system_data.write().unwrap().forces_mut(),
         self._particle_system_data.read().unwrap().positions(),
         self._particle_system_data.read().unwrap().velocities()).into_par_iter().for_each(|(force, pos, vel)| {
            // Gravity
            let mut f = self._gravity * mass;

            // Wind forces
            let relative_vel = *vel - self._wind.read().unwrap().sample(pos);
            f += relative_vel * -self._drag_coefficient;

            *force += f;
        });
    }

    fn time_integration(&mut self, time_step_in_seconds: f64) {
        let mass = self._particle_system_data.read().unwrap().mass();

        (&mut self._new_positions, &mut self._new_velocities,
         self._particle_system_data.read().unwrap().forces(),
         self._particle_system_data.read().unwrap().positions(),
         self._particle_system_data.read().unwrap().velocities()).into_par_iter().for_each(|(new_pos, new_vel, force, pos, vel)| {
            // Integrate velocity first
            *new_vel = *vel + *force * time_step_in_seconds / mass;

            // Integrate position.
            *new_pos = *pos + *new_vel * time_step_in_seconds;
        });
    }

    fn update_collider(&mut self, time_step_in_seconds: f64) {
        if let Some(collider) = &self._collider {
            collider.write().unwrap().update(self.current_time_in_seconds(),
                                             time_step_in_seconds);
        }
    }

    fn update_emitter(&mut self, time_step_in_seconds: f64) {
        if let Some(emitter) = &self._emitter {
            emitter.write().unwrap().update(self.current_time_in_seconds(),
                                            time_step_in_seconds);
        }
    }
}

/// Shared pointer type for the ParticleSystemSolver3.
pub type ParticleSystemSolver3Ptr = Arc<RwLock<ParticleSystemSolver3>>;

//--------------------------------------------------------------------------------------------------
///
/// # Base class for particle-based solver builder.
///
pub trait ParticleSystemSolverBuilderBase3 {
    /// Returns builder with particle radius.
    fn with_radius(&mut self, radius: f64) -> &mut Self;

    /// Returns builder with mass per particle.
    fn with_mass(&mut self, mass: f64) -> &mut Self;
}

///
/// # Front-end to create ParticleSystemSolver3 objects step by step.
///
pub struct Builder {
    _radius: f64,
    _mass: f64,
}

impl Builder {
    /// Builds ParticleSystemSolver3.
    pub fn build(&mut self) -> ParticleSystemSolver3 {
        return ParticleSystemSolver3::new(self._radius, self._mass);
    }

    /// Builds shared pointer of ParticleSystemSolver3 instance.
    pub fn make_shared(&mut self) -> ParticleSystemSolver3Ptr {
        return ParticleSystemSolver3Ptr::new(RwLock::new(self.build()));
    }

    /// constructor
    pub fn new() -> Builder {
        return Builder {
            _radius: 1.0e-3,
            _mass: 1.0e-3,
        };
    }
}

impl ParticleSystemSolverBuilderBase3 for Builder {
    fn with_radius(&mut self, radius: f64) -> &mut Self {
        self._radius = radius;
        return self;
    }

    fn with_mass(&mut self, mass: f64) -> &mut Self {
        self._mass = mass;
        return self;
    }
}