feat(nexus-mpm): add Material Point Method solver with one-way coupling

Cargo.toml

@@ -10,6 +10,10 @@ members = [
     "crates/nexus_rbd3d",
     "crates/nexus_rbd_shaders2d",
     "crates/nexus_rbd_shaders3d",
+    "crates/nexus_mpm_shaders2d",
+    "crates/nexus_mpm_shaders3d",
+    "crates/nexus_mpm2d",
+    "crates/nexus_mpm3d",
     "crates/nexus_python3d",
 ]
 resolver = "2"
@@ -68,13 +72,17 @@ crunchy = "0.2.4"
 # Shader crates
 nexus_rbd_shaders2d = { version = "0.4.0", path = "crates/nexus_rbd_shaders2d" }
 nexus_rbd_shaders3d = { version = "0.4.0", path = "crates/nexus_rbd_shaders3d" }
+nexus_mpm_shaders2d = { version = "0.4.0", path = "crates/nexus_mpm_shaders2d" }
+nexus_mpm_shaders3d = { version = "0.4.0", path = "crates/nexus_mpm_shaders3d" }
 
-# Internal crates. rbd is pulled with default-features off so dependents can
-# pick their own feature set, mirroring the rapier/parry pattern.
+# Internal crates. rbd is pulled both with defaults and with default-features
+# off (mpm wants just the dim feature), so it mirrors the rapier/parry pattern.
 nexus2d = { version = "0.4.0", path = "crates/nexus2d" }
 nexus3d = { version = "0.4.0", path = "crates/nexus3d" }
 nexus_rbd2d = { version = "0.4.0", path = "crates/nexus_rbd2d", default-features = false }
 nexus_rbd3d = { version = "0.4.0", path = "crates/nexus_rbd3d", default-features = false }
+nexus_mpm2d = { version = "0.4.0", path = "crates/nexus_mpm2d" }
+nexus_mpm3d = { version = "0.4.0", path = "crates/nexus_mpm3d" }
 nexus_viewer2d = { version = "0.4.0", path = "crates/nexus_viewer2d" }
 nexus_viewer3d = { version = "0.4.0", path = "crates/nexus_viewer3d" }
 

crates/examples2d/Cargo.toml

@@ -14,7 +14,7 @@ metal = ["nexus_viewer2d/metal"]
 
 [dependencies]
 glamx = { workspace = true }
-nexus2d = { workspace = true, features = [ "rbd" ] }
+nexus2d = { workspace = true, features = [ "rbd", "mpm" ] }
 nexus_viewer2d = { workspace = true }
 nexus_rbd2d = { workspace = true, features = ["default"] }
 rapier2d = { workspace = true, features = ["default"] }
@@ -24,6 +24,7 @@ kiss3d = { workspace = true }
 oorandom = { workspace = true }
 rand = "0.10"
 anyhow = { workspace = true }
+pollster = { workspace = true }
 
 [target.'cfg(target_arch = "wasm32")'.dependencies]
 getrandom = { workspace = true }

crates/examples2d/all_examples2.rs

@@ -2,17 +2,28 @@ use inflector::Inflector;
 use nexus_viewer2d::{BackendType, DemoKind, NexusViewer};
 use nexus2d::prelude::{NexusPipeline, NexusPipelineMask};
 
-mod balls2;
-mod boxes2;
-mod boxes_and_balls2;
-mod compound2;
-mod dynamic_rbd2;
-mod joint_ball2;
-mod joint_fixed2;
-mod joint_prismatic2;
-mod polyline2;
-mod primitives2;
-mod pyramid2;
+mod rbd_balls2;
+mod rbd_boxes2;
+mod rbd_boxes_and_balls2;
+mod rbd_compound2;
+mod rbd_dynamic2;
+mod rbd_joint_ball2;
+mod rbd_joint_fixed2;
+mod rbd_joint_prismatic2;
+mod rbd_polyline2;
+mod rbd_primitives2;
+mod rbd_pyramid2;
+
+// MPM examples.
+mod mpm_centilever_beam2;
+mod mpm_cohesion_sweep2;
+mod mpm_dam_break2;
+mod mpm_elastic_cut2;
+mod mpm_elasticity2;
+mod mpm_emitter2;
+mod mpm_hourglass2;
+mod mpm_sand2;
+mod mpm_snowball2;
 
 /// Declares the demo registry: a `(name, kind)` list for the picker UI and a
 /// name -> `run()` dispatcher. Keeping both in one macro keeps them in sync.
@@ -42,17 +53,27 @@ macro_rules! demos {
 }
 
 demos! {
-    "Balls" => Rbd : balls2,
-    "Boxes" => Rbd : boxes2,
-    "Boxes & balls" => Rbd : boxes_and_balls2,
-    "Compound" => Rbd : compound2,
-    "Dynamic insertion" => Rbd : dynamic_rbd2,
-    "Pyramid" => Rbd : pyramid2,
-    "Primitives" => Rbd : primitives2,
-    "Polyline" => Rbd : polyline2,
-    "Joints (spherical)" => Rbd : joint_ball2,
-    "Joints (prismatic)" => Rbd : joint_prismatic2,
-    "Joints (fixed)" => Rbd : joint_fixed2,
+    "Balls" => Rbd : rbd_balls2,
+    "Boxes" => Rbd : rbd_boxes2,
+    "Boxes & balls" => Rbd : rbd_boxes_and_balls2,
+    "Compound" => Rbd : rbd_compound2,
+    "Dynamic insertion" => Rbd : rbd_dynamic2,
+    "Pyramid" => Rbd : rbd_pyramid2,
+    "Primitives" => Rbd : rbd_primitives2,
+    "Polyline" => Rbd : rbd_polyline2,
+    "Joints (spherical)" => Rbd : rbd_joint_ball2,
+    "Joints (prismatic)" => Rbd : rbd_joint_prismatic2,
+    "Joints (fixed)" => Rbd : rbd_joint_fixed2,
+    // MPM demos.
+    "Cantilever beam" => Mpm : mpm_centilever_beam2,
+    "Sand" => Mpm : mpm_sand2,
+    "Sand emitter" => Mpm : mpm_emitter2,
+    "Elasticity" => Mpm : mpm_elasticity2,
+    "Elastic cut" => Mpm : mpm_elastic_cut2,
+    "Dam break" => Mpm : mpm_dam_break2,
+    "Cohesion sweep" => Mpm : mpm_cohesion_sweep2,
+    "Snowballs" => Mpm : mpm_snowball2,
+    "Hourglass" => Mpm : mpm_hourglass2,
 }
 
 struct CliOptions {

crates/examples2d/mpm_centilever_beam2.rs

@@ -0,0 +1,70 @@
+use khal::backend::GpuTimestamps;
+use nexus_viewer2d::NexusViewer;
+use nexus2d::mpm::solver::{BoundaryCondition, Particle, ParticleModel, SimulationParams};
+use nexus2d::prelude::{NexusPipeline, NexusState, RbdCoupling};
+
+use rapier2d::prelude::{ColliderBuilder, Pose, RigidBodyBuilder};
+
+pub async fn run(
+    viewer: &mut NexusViewer,
+    pipeline: &mut NexusPipeline,
+) -> anyhow::Result<NexusState> {
+    let mut state = NexusState::default();
+
+    let width = 10.0;
+    let height = 2.0;
+    let fixed_part = 1.0;
+    let cell_width = 0.2;
+    let particle_per_cell_dim = 2;
+    let young_modulus = 1.0e8;
+    let poisson_ratio = 0.3;
+
+    let diameter = cell_width / particle_per_cell_dim as f32;
+    let ni = ((width + fixed_part) / diameter).ceil() as usize;
+    let nj = (height / diameter).ceil() as usize;
+
+    let mut particles = vec![];
+    for i in 0..ni {
+        for j in 0..nj {
+            let position = glamx::vec2(i as f32, j as f32) * diameter;
+            let density = 1000.0;
+            let radius = diameter / 2.0;
+            let model = ParticleModel::elastic_neo_hookean(young_modulus, poisson_ratio);
+            particles.push(Particle::new(position, radius, density, model));
+        }
+    }
+
+    let params = SimulationParams {
+        gravity: glamx::vec2(0.0, -9.81),
+        padding: 0.0,
+        dt: 1.0 / 60.0,
+    };
+    state.set_mpm_params(viewer.backend(), params, cell_width)?;
+    state.set_mpm_substeps(150);
+    state.add_particles(viewer.backend(), particles)?;
+
+    // Fixed anchor the beam is cantilevered from (boundary coupled to the MPM
+    // continuum).
+    let body = RigidBodyBuilder::fixed()
+        .translation(glamx::vec2(0.0, height / 2.0))
+        .build();
+    let collider = ColliderBuilder::cuboid(fixed_part, height).build();
+    let shape = collider.shared_shape().clone();
+    let stick = BoundaryCondition::stick();
+    let handle = state.insert_rigid_body(body, collider, RbdCoupling::MpmOneWay(stick));
+    viewer.insert_shape(handle, &shape, Pose::IDENTITY);
+
+    let mut timestamps = GpuTimestamps::new(viewer.backend(), 2048);
+    state.finalize(viewer.backend()).await?;
+
+    while viewer.render_frame().await {
+        if viewer.simulating() {
+            pipeline
+                .simulate(viewer.backend(), &mut state, Some(&mut timestamps))
+                .await?;
+        }
+        viewer.sync(&mut state, Some(&mut timestamps)).await?;
+    }
+
+    Ok(state)
+}

crates/examples2d/mpm_cohesion_sweep2.rs

@@ -0,0 +1,99 @@
+//! Five identical granular columns released at once, differing only in cohesion.
+
+use khal::backend::GpuTimestamps;
+use nexus_viewer2d::NexusViewer;
+use nexus2d::mpm::solver::{BoundaryCondition, Particle, ParticleModel, SimulationParams};
+use nexus2d::prelude::{NexusPipeline, NexusState, RbdCoupling};
+
+use glamx::{Vec4, vec2};
+use rapier2d::prelude::{ColliderBuilder, Pose, RigidBodyBuilder};
+
+const DENSITY: f32 = 1600.0;
+const YOUNG_MODULUS: f32 = 1.0e7;
+const POISSON_RATIO: f32 = 0.2;
+
+/// Cohesion of each column, left to right.
+const COHESIONS: [f32; 5] = [0.0, 0.001, 0.003, 0.01, 0.03];
+
+/// Half-width of each column.
+const COLUMN_HALF_WIDTH: f32 = 3.0;
+/// Height of each column.
+const COLUMN_HEIGHT: f32 = 20.0;
+/// Distance between column centers.
+const COLUMN_PITCH: f32 = 22.0;
+
+pub async fn run(
+    viewer: &mut NexusViewer,
+    pipeline: &mut NexusPipeline,
+) -> anyhow::Result<NexusState> {
+    let mut state = NexusState::default();
+
+    let cell_width = 0.2;
+    let radius = cell_width / 4.0;
+    let spacing = radius * 2.0;
+
+    let mut particles = vec![];
+    let nx = (COLUMN_HALF_WIDTH * 2.0 / spacing) as i32;
+    let ny = (COLUMN_HEIGHT / spacing) as i32;
+    let count = COHESIONS.len();
+    // Dry sand is pale, the most cohesive column is dark.
+    let shades: Vec<_> = (0..count)
+        .map(|c| {
+            let t = c as f32 / (count - 1) as f32;
+            Vec4::new(0.92 - 0.5 * t, 0.80 - 0.48 * t, 0.62 - 0.42 * t, 1.0)
+        })
+        .collect();
+    viewer.set_particle_group_colors(&shades);
+    for (c, cohesion) in COHESIONS.iter().enumerate() {
+        let center_x = (c as f32 - (count - 1) as f32 / 2.0) * COLUMN_PITCH;
+        let model = ParticleModel::cohesive_sand(YOUNG_MODULUS, POISSON_RATIO, *cohesion);
+
+        for i in 0..nx {
+            for j in 0..ny {
+                let position = vec2(
+                    center_x - COLUMN_HALF_WIDTH + (i as f32 + 0.5) * spacing,
+                    0.2 + (j as f32 + 0.5) * spacing,
+                );
+                particles.push(Particle::with_group(
+                    position, radius, DENSITY, model, c as u32,
+                ));
+            }
+        }
+    }
+
+    let params = SimulationParams {
+        gravity: vec2(0.0, -9.81),
+        padding: 0.0,
+        dt: 1.0 / 60.0,
+    };
+    state.set_mpm_params(viewer.backend(), params, cell_width)?;
+    state.set_mpm_substeps(15);
+    state.add_particles(viewer.backend(), particles)?;
+
+    /*
+     * Setup the floor.
+     */
+    let half_span = COLUMN_PITCH * count as f32 / 2.0 + 10.0;
+    let collider = ColliderBuilder::cuboid(half_span, 1.0).build();
+    let shape = collider.shared_shape().clone();
+    let body = RigidBodyBuilder::fixed()
+        .translation(vec2(0.0, -1.0))
+        .build();
+    let coupling = RbdCoupling::MpmOneWay(BoundaryCondition::separate(1.0));
+    let handle = state.insert_rigid_body(body, collider, coupling);
+    viewer.insert_shape(handle, &shape, Pose::IDENTITY);
+
+    let mut timestamps = GpuTimestamps::new(viewer.backend(), 2048);
+    state.finalize(viewer.backend()).await?;
+
+    while viewer.render_frame().await {
+        if viewer.simulating() {
+            pipeline
+                .simulate(viewer.backend(), &mut state, Some(&mut timestamps))
+                .await?;
+        }
+        viewer.sync(&mut state, Some(&mut timestamps)).await?;
+    }
+
+    Ok(state)
+}