feat: add support for fixed particles

Author: sebcrozet

shaders/slosh/models/neo_hookean_elasticity.slang

@@ -2,6 +2,8 @@
 
 module neo_hookean_elasticity;
 
+import stensor.geometry.svd2;
+import stensor.geometry.svd3;
 import nexus.aliases;
 import slosh.models.utils;
 import slosh.models.interfaces;
@@ -41,4 +43,24 @@ public struct NeoHookeanModel {
         );
         return bound.timestep_bound(particle_density0, particle_def_grad_det, particle_velocity, cell_width);
     }
+
+    public func pos_energy_density(in def_grad: Mat<float>, elastic_hardening: float) -> float {
+        let j = determinant(def_grad);
+
+        #if DIM == 2
+        let svd = Svd2(def_grad);
+        #else
+        let svd = Svd3(def_grad);
+        #endif
+
+        let sig = max(svd.S - Vect<float>(1.0), Vect<float>(0.0));
+        let pos_dev_part = this.mu * elastic_hardening * dot(sig, sig);
+        let spherical_part = this.lambda * elastic_hardening * 0.5 * (j - 1.0) * (j - 1.0);
+
+        if (j < 1.0) {
+            return pos_dev_part;
+        } else {
+            return pos_dev_part + spherical_part;
+        }
+    }
 }

shaders/slosh/solver/particle.slang

@@ -30,7 +30,9 @@ public struct Dynamics {
     /// Rayleigh mass-proportional damping coefficient (1/s).
     public float damping;
     public float phase;
+    // TODO: merge these two into a `flags` field?
     public uint enabled;
+    public uint fixed;
 
     public func init_density() -> float {
         return this.mass / this.init_volume;

shaders/slosh/solver/particle_update.slang

@@ -53,6 +53,15 @@ func particle_update(
     var dynamics = particles_dyn[particle_id];
     let particle_pos = particles_pos[particle_id].pt;
 
+    // If the particle is fixed, clear its velocity and its affine matrix.
+    // This isn’t ideal (this should typically be handled on the grid) but we sometimes
+    // need sub-grid-sized fixed particles.
+    if (dynamics.fixed != 0) {
+        dynamics.velocity = Vect<float>(0.0);
+        dynamics.affine = Mat<float>(0.0);
+        return;
+    }
+
     /*
      * Advection.
      */

shaders_testbed/slosh_testbed/prep_readback2.slang

@@ -44,7 +44,7 @@ func prep_readback(
     if (particle_id < instances.getCount()) {
         let def_grad = particles_dyn[particle_id].def_grad;
         let init_deformation = diag(Vect<float>(particles_dyn[particle_id].init_radius * 2.0));
-        instances[particle_id].deformation = mul(init_deformation, def_grad);
+        instances[particle_id].deformation = mul(init_deformation, clamp(def_grad, Mat<float>(-4.0), Mat<float>(4.0)));
         instances[particle_id].position = float2(particles_pos[particle_id].pt);
 
         let color = base_colors[particle_id];

src/solver/particle.rs

@@ -53,6 +53,8 @@ pub struct ParticleDynamics {
     pub phase: f32,
     /// Whether this particle is active (1) or disabled (0).
     pub enabled: u32,
+    /// Whether this particle is fixed (1) or dynamic (0).
+    pub fixed: u32,
 }
 
 impl ParticleDynamics {
@@ -80,9 +82,15 @@ impl ParticleDynamics {
             cdf: Cdf::default(),
             phase: 1.0,
             enabled: 1,
+            fixed: 0,
         }
     }
 
+    /// Sets whether this particle is fixed (true) or dynamic (false).
+    pub fn set_fixed(&mut self, fixed: bool) {
+        self.fixed = fixed as u32;
+    }
+
     /// Sets the damping coefficient for this particle.
     pub fn set_damping(&mut self, damping: f32) {
         self.damping = damping;
@@ -481,6 +489,11 @@ impl<B: Backend, GpuModel: GpuParticleModelData> GpuParticles<B, GpuModel> {
         &self.models
     }
 
+    /// Returns mutable reference to material model buffer.
+    pub fn models_mut(&mut self) -> &mut GpuTensor<GpuModel, B> {
+        &mut self.models
+    }
+
     /// Returns reference to position buffer.
     pub fn positions(&self) -> &GpuTensor<ParticlePosition, B> {
         &self.positions
@@ -491,6 +504,11 @@ impl<B: Backend, GpuModel: GpuParticleModelData> GpuParticles<B, GpuModel> {
         &self.dynamics
     }
 
+    /// Returns mutable reference to dynamics buffer (velocity, deformation, mass).
+    pub fn dynamics_mut(&mut self) -> &mut GpuTensor<ParticleDynamics, B> {
+        &mut self.dynamics
+    }
+
     /// Returns reference to sorted particle ID buffer.
     pub fn sorted_ids(&self) -> &GpuTensor<u32, B> {
         &self.sorted_ids

src_testbed/lib.rs

@@ -82,6 +82,7 @@ impl<GpuModel: GpuParticleModelData> Stage<GpuModel> {
         let limits = Limits {
             max_storage_buffers_per_shader_stage: 11,
             max_compute_workgroup_storage_size: 32768, // Why do we need this if wgsparkl didn’t?
+            max_buffer_size: 1_000_000_000,
             max_storage_buffer_binding_size: 1_000_000_000,
             ..Limits::default()
         };