diff --git a/Cargo.toml b/Cargo.toml
index bd7f8f9..02e56b8 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -4,6 +4,8 @@ members = [
     "crates/wgsparkl-testbed3d",
     "crates/wgsparkl2d",
     "crates/wgsparkl3d",
+    "examples3d",
+    "examples2d",
 ]
 resolver = "2"
 
@@ -23,15 +25,32 @@ rust.unexpected_cfgs = { level = "warn", check-cfg = [
 ] }
 
 [patch.crates-io]
-parry2d = { path = "../parry/crates/parry2d" }
-parry3d = { path = "../parry/crates/parry3d" }
-encase = { path = "../encase" }
-wgcore = { path = "../wgmath/crates/wgcore" }
-wgebra = { path = "../wgmath/crates/wgebra" }
-wgparry2d = { path = "../wgmath/crates/wgparry/crates/wgparry2d" }
-wgparry3d = { path = "../wgmath/crates/wgparry/crates/wgparry3d" }
-wgrapier2d = { path = "../wgmath/crates/wgrapier/crates/wgrapier2d" }
-wgrapier3d = { path = "../wgmath/crates/wgrapier/crates/wgrapier3d" }
+# parry2d = { path = "../parry/crates/parry2d" }
+# parry3d = { path = "../parry/crates/parry3d" }
+# encase = { path = "../encase" }
+# wgcore = { path = "../wgmath/crates/wgcore" }
+# wgebra = { path = "../wgmath/crates/wgebra" }
+# wgparry2d = { path = "../wgmath/crates/wgparry/crates/wgparry2d" }
+# wgparry3d = { path = "../wgmath/crates/wgparry/crates/wgparry3d" }
+# wgrapier2d = { path = "../wgmath/crates/wgrapier/crates/wgrapier2d" }
+# wgrapier3d = { path = "../wgmath/crates/wgrapier/crates/wgrapier3d" }
+
+parry2d = { git = "https://github.com/dimforge/parry.git", rev = "407fb449504dde7ba0a3d80a8de7e37ac991a77a" }
+parry3d = { git = "https://github.com/dimforge/parry.git", rev = "407fb449504dde7ba0a3d80a8de7e37ac991a77a" }
+encase = { git = "https://github.com/sebcrozet/encase.git", rev = "6755c1414c225cbeac9fde1ff9394b70d3b05404" }
+# wgcore = { git = "https://github.com/dimforge/wgmath.git", rev = "32b54c802be6040c5f003e7d3def4a33e4730d92" }
+wgcore = { git = "https://github.com/dimforge/wgmath.git", branch = "wgml-chat" }
+wgcore-derive = { git = "https://github.com/dimforge/wgmath.git", branch = "wgml-chat" }
+wgebra = { git = "https://github.com/dimforge/wgmath.git", branch = "wgml-chat" }
+wgparry3d = { git = "https://github.com/dimforge/wgmath.git", branch = "wgml-chat" }
+wgparry2d = { git = "https://github.com/dimforge/wgmath.git", branch = "wgml-chat" }
+wgrapier3d = { git = "https://github.com/dimforge/wgmath.git", branch = "wgml-chat" }
+wgrapier2d = { git = "https://github.com/dimforge/wgmath.git", branch = "wgml-chat" }
+rapier3d = { git = "https://github.com/vrixyz/rapier.git", branch = "uber_physics_context" }
+rapier2d = { git = "https://github.com/vrixyz/rapier.git", branch = "uber_physics_context" }
 
 [profile.release]
 opt-level = 'z'
+
+[profile.dev]
+opt-level = 1
diff --git a/build_wasm.sh b/build_wasm.sh
index 8b2f27d..42cdf7e 100755
--- a/build_wasm.sh
+++ b/build_wasm.sh
@@ -1,11 +1,11 @@
 #!/bin/sh
 
-cargo build --example elasticity2 --release --target wasm32-unknown-unknown --features dim2
-cargo build --example elasticity3 --release --target wasm32-unknown-unknown --features dim3
-wasm-bindgen --no-typescript --target web --out-dir dist2d --out-name elasticity2 ./target/wasm32-unknown-unknown/release/examples/elasticity2.wasm
-wasm-bindgen --no-typescript --target web --out-dir dist3d --out-name elasticity3 ./target/wasm32-unknown-unknown/release/examples/elasticity3.wasm
-wasm-opt -Oz -o ./dist2d/opt.wasm ./dist2d/elasticity2_bg.wasm && mv ./dist2d/opt.wasm ./dist2d/elasticity2_bg.wasm
-wasm-opt -Oz -o ./dist3d/opt.wasm ./dist3d/elasticity3_bg.wasm && mv ./dist3d/opt.wasm ./dist3d/elasticity3_bg.wasm
+# cargo build -p wgsparkl-examples-2d --release --target wasm32-unknown-unknown
+cargo build -p wgsparkl-examples-3d --release --target wasm32-unknown-unknown
+# wasm-bindgen --no-typescript --target web --out-dir dist2d --out-name wgsparkl-examples-2d ./target/wasm32-unknown-unknown/release/wgsparkl-examples-2d.wasm
+wasm-bindgen --no-typescript --target web --out-dir dist3d --out-name wgsparkl-examples-3d ./target/wasm32-unknown-unknown/release/wgsparkl-examples-3d.wasm
+# wasm-opt -Oz -o ./dist2d/opt.wasm ./dist2d/wgsparkl-examples-2d_bg.wasm && mv ./dist2d/opt.wasm ./dist2d/wgsparkl-examples-2d_bg.wasm
+wasm-opt -Oz -o ./dist3d/opt.wasm ./dist3d/wgsparkl-examples-3d_bg.wasm && mv ./dist3d/opt.wasm ./dist3d/wgsparkl-examples-3d_bg.wasm
 
-brotli ./dist2d/elasticity2_bg.wasm && mv ./dist2d/elasticity2_bg.wasm.br ./dist2d/elasticity2_bg.wasm
-brotli ./dist3d/elasticity3_bg.wasm && mv ./dist3d/elasticity3_bg.wasm.br ./dist3d/elasticity3_bg.wasm
\ No newline at end of file
+# brotli ./dist2d/wgsparkl-examples-2d_bg.wasm && mv ./dist2d/wgsparkl-examples-2d_bg.wasm.br ./dist2d/wgsparkl-examples-2d_bg.wasm
+brotli ./dist3d/wgsparkl-examples-3d_bg.wasm && mv ./dist3d/wgsparkl-examples-3d_bg.wasm.br ./dist3d/wgsparkl-examples-3d_bg.wasm
\ No newline at end of file
diff --git a/crates/wgsparkl-testbed2d/Cargo.toml b/crates/wgsparkl-testbed2d/Cargo.toml
index a402a4c..483c7e8 100644
--- a/crates/wgsparkl-testbed2d/Cargo.toml
+++ b/crates/wgsparkl-testbed2d/Cargo.toml
@@ -13,8 +13,8 @@ path = "../../src_testbed/lib.rs"
 required-features = ["dim2"]
 
 [features]
+default = ["dim2"]
 dim2 = []
-default_features = ["dim2"]
 
 [dependencies]
 nalgebra = { workspace = true, features = ["rand"] }
@@ -31,10 +31,17 @@ futures-test = "0.3"
 serial_test = "3"
 approx = "0.5"
 async-std = { version = "1", features = ["attributes"] }
-bevy = { version = "0.15.0", features = ["shader_format_glsl", "shader_format_spirv", "webgpu"] }
+bevy = { version = "0.15.0", features = [
+    "shader_format_glsl",
+    "shader_format_spirv",
+    "webgpu",
+] }
 #bevy_wasm_window_resize = "0.4"
 bevy_editor_cam = "0.5"
-bevy_egui = { version = "0.31", default-features = false, features = ["default_fonts", "render"] }
+bevy_egui = { version = "0.31", default-features = false, features = [
+    "default_fonts",
+    "render",
+] }
 bevy_sprite = "0.15"
 wgsparkl2d = { path = "../wgsparkl2d" }
-futures = "0.3"
\ No newline at end of file
+futures = "0.3"
diff --git a/crates/wgsparkl-testbed3d/Cargo.toml b/crates/wgsparkl-testbed3d/Cargo.toml
index 3eb114a..d652796 100644
--- a/crates/wgsparkl-testbed3d/Cargo.toml
+++ b/crates/wgsparkl-testbed3d/Cargo.toml
@@ -13,8 +13,8 @@ path = "../../src_testbed/lib.rs"
 required-features = ["dim3"]
 
 [features]
+default = ["dim3"]
 dim3 = []
-default_features = ["dim3"]
 
 [dependencies]
 nalgebra = { workspace = true, features = ["rand"] }
@@ -31,10 +31,17 @@ futures-test = "0.3"
 serial_test = "3"
 approx = "0.5"
 async-std = { version = "1", features = ["attributes"] }
-bevy = { version = "0.15.0", features = ["shader_format_glsl", "shader_format_spirv", "webgpu"] }
+bevy = { version = "0.15.0", features = [
+    "shader_format_glsl",
+    "shader_format_spirv",
+    "webgpu",
+] }
 bevy_editor_cam = "0.5"
 #bevy_mod_picking = { version = "0.20", default-features = false }
 #bevy_wasm_window_resize = "0.4"
-bevy_egui = { version = "0.31", default-features = false, features = ["default_fonts", "render"] }
+bevy_egui = { version = "0.31", default-features = false, features = [
+    "default_fonts",
+    "render",
+] }
 wgsparkl3d = { path = "../wgsparkl3d" }
 futures = "0.3"
diff --git a/crates/wgsparkl2d/Cargo.toml b/crates/wgsparkl2d/Cargo.toml
index 8645eba..8602b88 100644
--- a/crates/wgsparkl2d/Cargo.toml
+++ b/crates/wgsparkl2d/Cargo.toml
@@ -34,12 +34,3 @@ rapier2d = "0.23"
 # For wasm?
 getrandom = { version = "0.3.1", features = ["wasm_js"] }
 uuid = { version = "1", features = ["js"] }
-
-[dev-dependencies]
-nalgebra = { version = "0.33", features = ["rand"] }
-futures-test = "0.3"
-serial_test = "3"
-approx = "0.5"
-async-std = { version = "1", features = ["attributes"] }
-bevy = { version = "0.15.0", features = ["shader_format_glsl", "shader_format_spirv", "webgpu"] }
-wgsparkl_testbed2d = { path = "../wgsparkl-testbed2d", features = ["dim2"] }
diff --git a/crates/wgsparkl2d/examples/testbed2.rs b/crates/wgsparkl2d/examples/testbed2.rs
deleted file mode 100644
index e5373f5..0000000
--- a/crates/wgsparkl2d/examples/testbed2.rs
+++ /dev/null
@@ -1,53 +0,0 @@
-use wgsparkl_testbed2d::{wgsparkl, RapierData};
-
-use bevy::prelude::*;
-use bevy::render::renderer::RenderDevice;
-use nalgebra::{point, vector, Similarity2, Vector2};
-use rapier2d::prelude::{ColliderBuilder, RigidBodyBuilder};
-use wgebra::GpuSim2;
-use wgparry2d::parry::shape::Cuboid;
-use wgrapier2d::dynamics::{BodyDesc, GpuVelocity};
-use wgsparkl::models::DruckerPrager;
-use wgsparkl::solver::ParticlePhase;
-use wgsparkl::{
-    models::ElasticCoefficients,
-    pipeline::MpmData,
-    solver::{Particle, SimulationParams},
-};
-use wgsparkl_testbed2d::{init_testbed, AppState, PhysicsContext, SceneInits};
-
-mod elastic_cut2;
-mod elasticity2;
-mod sand2;
-
-pub fn main() {
-    let mut app = App::new();
-    init_testbed(&mut app);
-    app.add_systems(
-        Startup,
-        (register_scenes, start_default_scene)
-            .chain()
-            .after(wgsparkl_testbed2d::startup::setup_app),
-    );
-    app.run();
-}
-
-fn register_scenes(world: &mut World) {
-    let scenes = vec![
-        ("sand".to_string(), world.register_system(sand2::sand_demo)),
-        (
-            "elastic".to_string(),
-            world.register_system(elasticity2::elastic_demo),
-        ),
-        (
-            "elastic cut".to_string(),
-            world.register_system(elastic_cut2::elastic_cut_demo),
-        ),
-    ];
-    let mut inits = world.resource_mut::<SceneInits>();
-    inits.scenes = scenes;
-}
-
-fn start_default_scene(mut commands: Commands, scenes: Res<SceneInits>) {
-    scenes.init_scene(&mut commands, 0);
-}
diff --git a/crates/wgsparkl3d/Cargo.toml b/crates/wgsparkl3d/Cargo.toml
index 6654d2e..e77cfac 100644
--- a/crates/wgsparkl3d/Cargo.toml
+++ b/crates/wgsparkl3d/Cargo.toml
@@ -15,6 +15,9 @@ required-features = ["dim3"]
 [features]
 default = ["dim3"]
 dim3 = []
+load_gltf = ["obj-rs", "image"]
+load_obj = ["gltf", "image"]
+load_bevy = ["bevy", "image"]
 
 [dependencies]
 nalgebra = { workspace = true }
@@ -33,11 +36,19 @@ wgrapier3d = "0.2"
 # TODO: make rapier optional?
 rapier3d = "0.23"
 
-[dev-dependencies]
-nalgebra = { version = "0.33", features = ["rand"] }
-futures-test = "0.3"
-serial_test = "3"
-approx = "0.5"
-async-std = { version = "1", features = ["attributes"] }
-bevy = { version = "0.15.0", features = ["shader_format_glsl", "shader_format_spirv", "webgpu"] }
-wgsparkl_testbed3d = { path = "../wgsparkl-testbed3d" }
\ No newline at end of file
+# for loaders
+# TODO: slim bevy features
+bevy = { version = "0.15.0", features = [
+    "shader_format_glsl",
+    "shader_format_spirv",
+    "webgpu",
+    "jpeg",
+], optional = true }
+gltf = { version = "1.3", features = [
+    "KHR_materials_pbrSpecularGlossiness",
+], optional = true }
+obj-rs = { version = "0.7", default-features = false, optional = true }
+image = { version = "0.25", optional = true }
+
+[package.metadata.docs.rs]
+cargo-args = ["-Zunstable-options", "-Zrustdoc-scrape-examples"]
diff --git a/crates/wgsparkl3d/examples/sand3.rs b/crates/wgsparkl3d/examples/sand3.rs
deleted file mode 100644
index f4eae09..0000000
--- a/crates/wgsparkl3d/examples/sand3.rs
+++ /dev/null
@@ -1,121 +0,0 @@
-use wgsparkl_testbed3d::{wgsparkl, RapierData};
-
-use bevy::prelude::*;
-use bevy::render::renderer::RenderDevice;
-use nalgebra::{vector, Similarity3, Vector3};
-use rapier3d::prelude::{ColliderBuilder, RigidBodyBuilder};
-use wgebra::GpuSim3;
-use wgparry3d::parry::shape::Cuboid;
-use wgrapier3d::dynamics::BodyDesc;
-use wgsparkl::models::DruckerPrager;
-use wgsparkl::{
-    models::ElasticCoefficients,
-    pipeline::MpmData,
-    solver::{Particle, ParticleDynamics, ParticlePhase, SimulationParams},
-};
-use wgsparkl_testbed3d::{init_testbed, AppState, PhysicsContext, SceneInits};
-
-fn main() {
-    panic!("Run the `testbed3` example instead.");
-}
-
-pub fn sand_demo(
-    mut commands: Commands,
-    device: Res<RenderDevice>,
-    mut app_state: ResMut<AppState>,
-) {
-    let mut rapier_data = RapierData::default();
-    let device = device.wgpu_device();
-
-    let nxz = 45;
-    let cell_width = 1.0;
-    let mut particles = vec![];
-    for i in 0..nxz {
-        for j in 0..100 {
-            for k in 0..nxz {
-                let position = vector![
-                    i as f32 + 0.5 - nxz as f32 / 2.0,
-                    j as f32 + 0.5 + 10.0,
-                    k as f32 + 0.5 - nxz as f32 / 2.0
-                ] * cell_width
-                    / 2.0;
-                let density = 2700.0;
-                let radius = cell_width / 4.0;
-                particles.push(Particle {
-                    position,
-                    dynamics: ParticleDynamics::with_density(radius, density),
-                    model: ElasticCoefficients::from_young_modulus(2_000_000_000.0, 0.2),
-                    plasticity: Some(DruckerPrager::new(2_000_000_000.0, 0.2)),
-                    phase: None,
-                });
-            }
-        }
-    }
-
-    if !app_state.restarting {
-        app_state.num_substeps = 20;
-        app_state.gravity_factor = 1.0;
-    };
-
-    let params = SimulationParams {
-        gravity: vector![0.0, -9.81, 0.0] * app_state.gravity_factor,
-        dt: (1.0 / 60.0) / (app_state.num_substeps as f32),
-    };
-
-    let rb = RigidBodyBuilder::fixed().translation(vector![0.0, -4.0, 0.0]);
-    let rb_handle = rapier_data.bodies.insert(rb);
-    let co = ColliderBuilder::cuboid(100.0, 4.0, 100.0);
-    rapier_data
-        .colliders
-        .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
-
-    let rb = RigidBodyBuilder::fixed().translation(vector![0.0, 5.0, -35.0]);
-    let rb_handle = rapier_data.bodies.insert(rb);
-    let co = ColliderBuilder::cuboid(35.0, 5.0, 0.5);
-    rapier_data
-        .colliders
-        .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
-    let rb = RigidBodyBuilder::fixed().translation(vector![0.0, 5.0, 35.0]);
-    let rb_handle = rapier_data.bodies.insert(rb);
-    let co = ColliderBuilder::cuboid(35.0, 5.0, 0.5);
-    rapier_data
-        .colliders
-        .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
-    let rb = RigidBodyBuilder::fixed().translation(vector![-35.0, 5.0, 0.0]);
-    let rb_handle = rapier_data.bodies.insert(rb);
-    let co = ColliderBuilder::cuboid(0.5, 5.0, 35.0);
-    rapier_data
-        .colliders
-        .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
-    let rb = RigidBodyBuilder::fixed().translation(vector![35.0, 5.0, 0.0]);
-    let rb_handle = rapier_data.bodies.insert(rb);
-    let co = ColliderBuilder::cuboid(0.5, 5.0, 35.0);
-    rapier_data
-        .colliders
-        .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
-
-    let rb = RigidBodyBuilder::kinematic_velocity_based()
-        .translation(vector![0.0, 2.0, 0.0])
-        .rotation(vector![0.0, 0.0, -0.5])
-        .angvel(vector![0.0, -1.0, 0.0]);
-    let rb_handle = rapier_data.bodies.insert(rb);
-    let co = ColliderBuilder::cuboid(0.5, 2.0, 30.0);
-    rapier_data
-        .colliders
-        .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
-
-    let data = MpmData::new(
-        device,
-        params,
-        &particles,
-        &rapier_data.bodies,
-        &rapier_data.colliders,
-        cell_width,
-        60_000,
-    );
-    commands.insert_resource(PhysicsContext {
-        data,
-        rapier_data,
-        particles,
-    });
-}
diff --git a/dist2d/index.html b/dist2d/index.html
index 0bb9815..3794f6a 100644
--- a/dist2d/index.html
+++ b/dist2d/index.html
@@ -3,7 +3,7 @@
 
 <body style="margin: 0px;">
     <script type="module">
-        import init from './elasticity2.js'
+        import init from './wgsparkl-examples-2d.js'
 
         init().catch((error) => {
             if (!error.message.startsWith("Using exceptions for control flow, don't mind me. This isn't actually an error!")) {
diff --git a/dist3d/index.html b/dist3d/index.html
index 51922e3..48feb87 100644
--- a/dist3d/index.html
+++ b/dist3d/index.html
@@ -3,7 +3,7 @@
 
 <body style="margin: 0px;">
     <script type="module">
-        import init from './elasticity3.js'
+        import init from './wgsparkl-examples-3d.js'
 
         init().catch((error) => {
             if (!error.message.startsWith("Using exceptions for control flow, don't mind me. This isn't actually an error!")) {
diff --git a/examples2d/Cargo.toml b/examples2d/Cargo.toml
new file mode 100644
index 0000000..8020513
--- /dev/null
+++ b/examples2d/Cargo.toml
@@ -0,0 +1,24 @@
+[package]
+name = "wgsparkl-examples-2d"
+version = "0.1.0"
+license = "Apache-2.0 OR Custom"
+edition = "2021"
+
+[dependencies]
+nalgebra = { version = "0.33", features = ["rand"] }
+futures-test = "0.3"
+serial_test = "3"
+approx = "0.5"
+async-std = { version = "1", features = ["attributes"] }
+bevy = { version = "0.15.0", features = [
+    "shader_format_glsl",
+    "shader_format_spirv",
+    "webgpu",
+    "jpeg",
+] }
+wgsparkl_testbed2d = { path = "../crates/wgsparkl-testbed2d" }
+wgsparkl2d = { path = "../crates/wgsparkl2d" }
+rapier2d = "0.23"
+
+[lints]
+workspace = true
diff --git a/crates/wgsparkl2d/examples/elastic_cut2.rs b/examples2d/src/elastic_cut2.rs
similarity index 86%
rename from crates/wgsparkl2d/examples/elastic_cut2.rs
rename to examples2d/src/elastic_cut2.rs
index 643b3c5..be78844 100644
--- a/crates/wgsparkl2d/examples/elastic_cut2.rs
+++ b/examples2d/src/elastic_cut2.rs
@@ -1,13 +1,8 @@
-use wgsparkl_testbed2d::{wgsparkl, RapierData};
+use wgsparkl_testbed2d::{wgsparkl, Callbacks, RapierData};
 
-use bevy::prelude::*;
 use bevy::render::renderer::RenderDevice;
-use nalgebra::{point, vector, Similarity2, Vector2};
+use nalgebra::{point, vector, Vector2};
 use rapier2d::prelude::{ColliderBuilder, RigidBodyBuilder};
-use wgebra::GpuSim2;
-use wgparry2d::parry::shape::Cuboid;
-use wgrapier2d::dynamics::{BodyDesc, GpuVelocity};
-use wgsparkl::models::DruckerPrager;
 use wgsparkl::solver::ParticlePhase;
 use wgsparkl::{
     models::ElasticCoefficients,
@@ -15,17 +10,17 @@ use wgsparkl::{
     solver::{Particle, SimulationParams},
 };
 use wgsparkl2d::solver::ParticleDynamics;
-use wgsparkl_testbed2d::{init_testbed, AppState, PhysicsContext, SceneInits};
+use wgsparkl_testbed2d::{AppState, PhysicsContext};
 
 fn main() {
     panic!("Run the `testbed3` example instead.");
 }
 
 pub fn elastic_cut_demo(
-    mut commands: Commands,
-    device: Res<RenderDevice>,
-    mut app_state: ResMut<AppState>,
-) {
+    device: RenderDevice,
+    app_state: &mut AppState,
+    _callbacks: &mut Callbacks,
+) -> PhysicsContext {
     let mut rapier_data = RapierData::default();
     let device = device.wgpu_device();
 
@@ -48,6 +43,7 @@ pub fn elastic_cut_demo(
                     phase: 1.0,
                     max_stretch: f32::MAX,
                 }),
+                color: None,
             });
         }
     }
@@ -118,9 +114,9 @@ pub fn elastic_cut_demo(
         cell_width,
         60_000,
     );
-    commands.insert_resource(PhysicsContext {
+    PhysicsContext {
         data,
         rapier_data,
         particles,
-    });
+    }
 }
diff --git a/crates/wgsparkl2d/examples/elasticity2.rs b/examples2d/src/elasticity2.rs
similarity index 84%
rename from crates/wgsparkl2d/examples/elasticity2.rs
rename to examples2d/src/elasticity2.rs
index 30019a7..d4203fd 100644
--- a/crates/wgsparkl2d/examples/elasticity2.rs
+++ b/examples2d/src/elasticity2.rs
@@ -1,13 +1,8 @@
-use wgsparkl_testbed2d::{wgsparkl, RapierData};
+use wgsparkl_testbed2d::{wgsparkl, Callbacks, RapierData};
 
-use bevy::prelude::*;
 use bevy::render::renderer::RenderDevice;
-use nalgebra::{point, vector, Similarity2, Vector2};
+use nalgebra::{vector, Vector2};
 use rapier2d::prelude::{ColliderBuilder, RigidBodyBuilder};
-use wgebra::GpuSim2;
-use wgparry2d::parry::shape::Cuboid;
-use wgrapier2d::dynamics::{BodyDesc, GpuVelocity};
-use wgsparkl::models::DruckerPrager;
 use wgsparkl::solver::ParticlePhase;
 use wgsparkl::{
     models::ElasticCoefficients,
@@ -15,17 +10,17 @@ use wgsparkl::{
     solver::{Particle, SimulationParams},
 };
 use wgsparkl2d::solver::ParticleDynamics;
-use wgsparkl_testbed2d::{init_testbed, AppState, PhysicsContext, SceneInits};
+use wgsparkl_testbed2d::{AppState, PhysicsContext};
 
 fn main() {
     panic!("Run the `testbed3` example instead.");
 }
 
 pub fn elastic_demo(
-    mut commands: Commands,
-    device: Res<RenderDevice>,
-    mut app_state: ResMut<AppState>,
-) {
+    device: RenderDevice,
+    app_state: &mut AppState,
+    _callbacks: &mut Callbacks,
+) -> PhysicsContext {
     let device = device.wgpu_device();
     let mut rapier_data = RapierData::default();
 
@@ -49,6 +44,7 @@ pub fn elastic_demo(
                     phase: 1.0,
                     max_stretch: f32::MAX,
                 }),
+                color: None,
             });
         }
     }
@@ -98,9 +94,9 @@ pub fn elastic_demo(
         cell_width,
         60_000,
     );
-    commands.insert_resource(PhysicsContext {
+    PhysicsContext {
         data,
         rapier_data,
         particles,
-    });
+    }
 }
diff --git a/examples2d/src/main.rs b/examples2d/src/main.rs
new file mode 100644
index 0000000..a2ce6a1
--- /dev/null
+++ b/examples2d/src/main.rs
@@ -0,0 +1,37 @@
+use wgsparkl_testbed2d::SceneInitFn;
+
+use bevy::prelude::*;
+use wgsparkl_testbed2d::{init_testbed, SceneInits};
+
+mod elastic_cut2;
+mod elasticity2;
+mod sand2;
+
+pub fn main() {
+    let mut app = App::new();
+    init_testbed(&mut app);
+    app.add_systems(
+        Startup,
+        (register_scenes, start_default_scene)
+            .chain()
+            .after(wgsparkl_testbed2d::startup::setup_app),
+    );
+    app.run();
+}
+
+fn register_scenes(world: &mut World) {
+    let scenes: Vec<(String, SceneInitFn)> = vec![
+        ("sand".to_string(), Box::new(sand2::sand_demo)),
+        ("elastic".to_string(), Box::new(elasticity2::elastic_demo)),
+        (
+            "elastic cut".to_string(),
+            Box::new(elastic_cut2::elastic_cut_demo),
+        ),
+    ];
+    let mut inits = world.resource_mut::<SceneInits>();
+    inits.scenes = scenes;
+}
+
+fn start_default_scene(mut commands: Commands, scenes: Res<SceneInits>) {
+    scenes.init_scene(&mut commands, 0);
+}
diff --git a/crates/wgsparkl2d/examples/sand2.rs b/examples2d/src/sand2.rs
similarity index 91%
rename from crates/wgsparkl2d/examples/sand2.rs
rename to examples2d/src/sand2.rs
index 357bf7e..3b73fb0 100644
--- a/crates/wgsparkl2d/examples/sand2.rs
+++ b/examples2d/src/sand2.rs
@@ -1,31 +1,26 @@
-use wgsparkl_testbed2d::{wgsparkl, RapierData};
+use wgsparkl_testbed2d::{wgsparkl, Callbacks, RapierData};
 
-use bevy::prelude::*;
 use bevy::render::renderer::RenderDevice;
-use nalgebra::{point, vector, Similarity2, Vector2};
+use nalgebra::{vector, Vector2};
 use rapier2d::prelude::{ColliderBuilder, RigidBodyBuilder};
-use wgebra::GpuSim2;
-use wgparry2d::parry::shape::Cuboid;
-use wgrapier2d::dynamics::{BodyDesc, GpuVelocity};
 use wgsparkl::models::DruckerPrager;
-use wgsparkl::solver::ParticlePhase;
 use wgsparkl::{
     models::ElasticCoefficients,
     pipeline::MpmData,
     solver::{Particle, SimulationParams},
 };
 use wgsparkl2d::solver::ParticleDynamics;
-use wgsparkl_testbed2d::{init_testbed, AppState, PhysicsContext, SceneInits};
+use wgsparkl_testbed2d::{AppState, PhysicsContext};
 
 fn main() {
     panic!("Run the `testbed3` example instead.");
 }
 
 pub fn sand_demo(
-    mut commands: Commands,
-    device: Res<RenderDevice>,
-    mut app_state: ResMut<AppState>,
-) {
+    device: RenderDevice,
+    app_state: &mut AppState,
+    _callbacks: &mut Callbacks,
+) -> PhysicsContext {
     let mut rapier_data = RapierData::default();
     let device = device.wgpu_device();
 
@@ -45,6 +40,7 @@ pub fn sand_demo(
                 model: ElasticCoefficients::from_young_modulus(10_000_000.0, 0.2),
                 plasticity: Some(DruckerPrager::new(10_000_000.0, 0.2)),
                 phase: None,
+                color: None,
             });
         }
     }
@@ -176,9 +172,9 @@ pub fn sand_demo(
         cell_width,
         60_000,
     );
-    commands.insert_resource(PhysicsContext {
+    PhysicsContext {
         data,
         rapier_data,
         particles,
-    });
+    }
 }
diff --git a/examples3d/Cargo.toml b/examples3d/Cargo.toml
new file mode 100644
index 0000000..4bd9a89
--- /dev/null
+++ b/examples3d/Cargo.toml
@@ -0,0 +1,32 @@
+[package]
+name = "wgsparkl-examples-3d"
+version = "0.1.0"
+license = "Apache-2.0 OR Custom"
+edition = "2021"
+
+[dependencies]
+nalgebra = { version = "0.33", features = ["rand"] }
+futures-test = "0.3"
+serial_test = "3"
+approx = "0.5"
+async-std = { version = "1", features = ["attributes"] }
+bevy = { version = "0.15.0", features = [
+    "shader_format_glsl",
+    "shader_format_spirv",
+    "webgpu",
+    "jpeg",
+] }
+wgsparkl_testbed3d = { path = "../crates/wgsparkl-testbed3d" }
+wgsparkl3d = { path = "../crates/wgsparkl3d", features = [
+    "load_gltf",
+    "load_bevy",
+    "load_obj",
+] }
+obj-rs = { version = "0.7", default-features = false }
+bevy_editor_cam = "*"
+image = "0.25"
+gltf = { version = "1.3", features = ["KHR_materials_pbrSpecularGlossiness"] }
+rapier3d = "0.23"
+
+[lints]
+workspace = true
diff --git a/examples3d/assets/banana.glb b/examples3d/assets/banana.glb
new file mode 100644
index 0000000..104c1d6
Binary files /dev/null and b/examples3d/assets/banana.glb differ
diff --git a/examples3d/assets/chefs_knife_open_blade.glb b/examples3d/assets/chefs_knife_open_blade.glb
new file mode 100644
index 0000000..7374bc9
Binary files /dev/null and b/examples3d/assets/chefs_knife_open_blade.glb differ
diff --git a/examples3d/assets/credits.txt b/examples3d/assets/credits.txt
new file mode 100644
index 0000000..c7edf64
--- /dev/null
+++ b/examples3d/assets/credits.txt
@@ -0,0 +1,8 @@
+# Thanks 
+
+- Shiba from zixisun02 at https://sketchfab.com/3d-models/shiba-faef9fe5ace445e7b2989d1c1ece361c
+- Banana from ReCo_3D at https://sketchfab.com/3d-models/banana-ed0537c2201d4ad7b8ddb8368d844174
+- Chef's knife from Davide Specchi at https://sketchfab.com/3d-models/chefs-knife-85ededce157545e68e51fa87768a08cb
+  - modifications:
+    - set origin to the middle of the blade + rotated so blade is facing Y down.
+    - inverted faces on a blade side to properly cut
diff --git a/examples3d/assets/shiba.glb b/examples3d/assets/shiba.glb
new file mode 100644
index 0000000..1fc2d6b
Binary files /dev/null and b/examples3d/assets/shiba.glb differ
diff --git a/examples3d/examples/glb_to_point_cloud.rs b/examples3d/examples/glb_to_point_cloud.rs
new file mode 100644
index 0000000..a8bd54f
--- /dev/null
+++ b/examples3d/examples/glb_to_point_cloud.rs
@@ -0,0 +1,71 @@
+//! This example demonstrates how to load a GLB file and display it as a point cloud.
+//!
+//! No MPM simulation is involved here.
+
+use bevy::{
+    app::{App, Startup, Update},
+    core_pipeline::core_3d::Camera3d,
+    ecs::system::Commands,
+    math::Vec3,
+    prelude::*,
+    render::camera::Camera,
+    DefaultPlugins,
+};
+use bevy_editor_cam::{prelude::EditorCam, DefaultEditorCamPlugins};
+use nalgebra::{Isometry3, Transform3, UnitQuaternion, Vector3};
+use std::{fs::File, io::Read};
+use wgsparkl3d::load_mesh3d::load_gltf::load_model_with_colors;
+
+fn main() {
+    App::new()
+        .add_plugins((DefaultPlugins, DefaultEditorCamPlugins))
+        .add_systems(Startup, init_scene)
+        .add_systems(Update, display_point_cloud)
+        .run();
+}
+
+#[derive(Component)]
+pub struct PointCloud {
+    pub positions: Vec<(Vec3, Color)>,
+}
+
+pub const SAMPLE_PER_UNIT: f32 = 20.0;
+
+fn init_scene(mut commands: Commands) {
+    commands.spawn((Camera3d::default(), Camera::default(), EditorCam::default()));
+
+    let mut file = File::open("assets/shiba.glb").expect("Failed to open GLB file");
+    let mut buffer = Vec::new();
+    file.read_to_end(&mut buffer).expect("Failed to read file");
+    let pc_grid = load_model_with_colors(
+        &buffer,
+        Transform3::from_matrix_unchecked(
+            Isometry3::from_parts(
+                Vector3::new(0.0, 0.0, 0.0).into(),
+                UnitQuaternion::identity(),
+            )
+            .to_matrix()
+            .scale(3.0),
+        ),
+        None,
+        SAMPLE_PER_UNIT,
+    );
+
+    commands.spawn(PointCloud {
+        positions: pc_grid
+            .iter()
+            .map(|p| {
+                let pos = Vec3::new(p.0.x, p.0.y, p.0.z);
+                (pos, Color::from(Srgba::from_u8_array(p.1)))
+            })
+            .collect::<Vec<_>>(),
+    });
+}
+
+fn display_point_cloud(pcs: Query<&PointCloud>, mut gizmos: Gizmos) {
+    for pc in pcs.iter() {
+        for p in pc.positions.iter() {
+            gizmos.sphere(p.0, 1f32 / SAMPLE_PER_UNIT / 2f32, p.1);
+        }
+    }
+}
diff --git a/examples3d/src/banana3.rs b/examples3d/src/banana3.rs
new file mode 100644
index 0000000..feb0223
--- /dev/null
+++ b/examples3d/src/banana3.rs
@@ -0,0 +1,182 @@
+use std::fs::File;
+use std::io::Read;
+
+use crate::utils::default_scene::{self, SAMPLE_PER_UNIT};
+
+use bevy::color::palettes::css;
+use bevy::{prelude::*, render::renderer::RenderDevice};
+use nalgebra::{vector, Isometry3, Transform3, UnitQuaternion, Vector3};
+use rapier3d::prelude::RigidBodyHandle;
+use rapier3d::prelude::{ColliderBuilder, RigidBodyBuilder, SharedShape, TriMeshFlags};
+use wgsparkl3d::load_mesh3d::load_gltf::{load_model_trimeshes, load_model_with_colors};
+use wgsparkl3d::pipeline::MpmData;
+use wgsparkl3d::solver::SimulationParams;
+use wgsparkl_testbed3d::Callbacks;
+use wgsparkl_testbed3d::{AppState, PhysicsContext};
+use wgsparkl_testbed3d::{Callback, RapierData};
+
+pub fn demo(
+    device: RenderDevice,
+    app_state: &mut AppState,
+    callbacks: &mut Callbacks,
+) -> PhysicsContext {
+    let mut file = File::open("assets/banana.glb").expect("Failed to open GLB file");
+    let mut buffer = Vec::new();
+    file.read_to_end(&mut buffer).expect("Failed to read file");
+
+    let pc_grid = load_model_with_colors(
+        &buffer,
+        Transform3::from_matrix_unchecked(
+            Isometry3::from_parts(
+                Vector3::new(0.0, 2.4, 0.0).into(),
+                UnitQuaternion::from_axis_angle(&Vector3::y_axis(), -90f32.to_radians()),
+            )
+            .to_matrix()
+                * nalgebra::Matrix4::new_scaling(0.3),
+        ),
+        Some(css::BLANCHED_ALMOND.to_u8_array()),
+        SAMPLE_PER_UNIT,
+    )
+    .iter()
+    .map(|p| {
+        (
+            Vec3::new(p.0.x, p.0.y, p.0.z),
+            Color::from(Srgba::from_u8_array(p.1)),
+        )
+    })
+    .collect::<Vec<_>>();
+
+    let params = SimulationParams {
+        gravity: vector![0.0, -9.81, 0.0] * app_state.gravity_factor,
+        dt: (1.0 / 60.0) / (app_state.num_substeps as f32),
+    };
+    let mut rapier_data = RapierData::default();
+    if !app_state.restarting {
+        app_state.num_substeps = 30;
+        app_state.gravity_factor = 1.0;
+    };
+    default_scene::spawn_ground_and_walls(&mut rapier_data);
+
+    let mut particles = vec![];
+    for (pos, color) in pc_grid {
+        let particle =
+            default_scene::create_particle(&pos, Some(color), 1f32 / SAMPLE_PER_UNIT / 2f32);
+        particles.push(particle);
+    }
+    // Slicer
+    let mut file =
+        File::open("assets/chefs_knife_open_blade.glb").expect("Failed to open GLB file");
+    let mut buffer = Vec::new();
+    file.read_to_end(&mut buffer).expect("Failed to read file");
+    let mut slicer_trimeshes = load_model_trimeshes(&buffer);
+    slicer_trimeshes.iter_mut().for_each(|trimesh| {
+        trimesh.0.iter_mut().for_each(|v| {
+            *v *= 10f32;
+        });
+    });
+
+    let rb = RigidBodyBuilder::kinematic_position_based().translation(vector![0.0, 5.0, 0.0]);
+    let parent_handle = rapier_data.bodies.insert(rb);
+    /*
+    // NOTE: Heightfield may be more efficient or more predictable than trimesh for the knife.
+    let heights = nalgebra::DMatrix::zeros(10, 5);
+    let heightfield = rapier3d::prelude::HeightField::new(heights, vector![2.0, 1.0, 10.0]);
+    let (mut vtx, idx) = heightfield.to_trimesh();
+    vtx.iter_mut().for_each(|pt| {
+        *pt = Isometry3::rotation(vector![0f32, 0f32, -90f32.to_radians()]) * *pt
+            + vector![0.0, 1f32, 0f32]
+    });
+    let co = ColliderBuilder::trimesh(vtx, idx).unwrap();
+    rapier_data
+        .colliders
+        .insert_with_parent(co, parent_handle, &mut rapier_data.bodies);
+    */
+    for (vertices, indices) in slicer_trimeshes.iter() {
+        // Insert collider into rapier state.
+
+        let collider = ColliderBuilder::new(
+            SharedShape::trimesh_with_flags(
+                vertices.clone(),
+                indices
+                    .chunks_exact(3)
+                    .map(|i| [i[0] as u32, i[1] as u32, i[2] as u32])
+                    .collect(),
+                TriMeshFlags::FIX_INTERNAL_EDGES,
+            )
+            .unwrap(),
+        );
+        rapier_data
+            .colliders
+            .insert_with_parent(collider, parent_handle, &mut rapier_data.bodies);
+    }
+
+    let data = MpmData::new(
+        device.wgpu_device(),
+        params,
+        &particles,
+        &rapier_data.bodies,
+        &rapier_data.colliders,
+        1f32 / default_scene::SAMPLE_PER_UNIT,
+        60_000,
+    );
+
+    let new_physics = PhysicsContext {
+        data,
+        rapier_data,
+        particles,
+    };
+    //
+    callbacks.0.push(move_knife_function(parent_handle));
+
+    new_physics
+}
+
+fn move_knife_function(body_handle: RigidBodyHandle) -> Callback {
+    Box::new(
+        move |_render, physics: &mut PhysicsContext, _timestamps, app_state: &AppState| {
+            let t = app_state.physics_time_seconds as f32;
+
+            let body = physics.rapier_data.bodies.get_mut(body_handle).unwrap();
+            let length = 1.31;
+            let width = 0.35;
+            let x_pos = 2.8;
+            let y_pos = 1.3;
+            let z_pos = -1.5;
+            let velocity = 0.9;
+            let extended_width = 0.15; // extended width for the horizontal move to the right
+            let period =
+                (2f32 * length + width + extended_width + 2f32 * width + extended_width) / velocity;
+
+            let i = (t / period).floor();
+            let dis = velocity * (t - period * i);
+
+            // Determine the new position of the knife based on the current displacement
+            let new_pos = if dis < length {
+                // Moving downwards
+                vector![x_pos - width * i, y_pos - dis, z_pos]
+            } else if length <= dis && dis < length + width + extended_width {
+                // Moving horizontally to the right with extended width
+                vector![x_pos - width * i + (dis - length), y_pos - length, z_pos]
+            } else if length + width + extended_width <= dis
+                && dis < 2f32 * length + width + extended_width
+            {
+                // Moving upwards
+                vector![
+                    x_pos - width * i + width + extended_width,
+                    y_pos - length + dis - (length + width + extended_width),
+                    z_pos
+                ]
+            } else {
+                // Moving horizontally to the left, back to the starting x position
+                vector![
+                    x_pos - width * i + width + extended_width
+                        - (dis - 2.0 * length - width - extended_width),
+                    y_pos,
+                    z_pos
+                ]
+            };
+
+            body.set_translation(new_pos, true);
+        },
+    )
+}
diff --git a/crates/wgsparkl3d/examples/elastic_cut3.rs b/examples3d/src/elastic_cut3.rs
similarity index 84%
rename from crates/wgsparkl3d/examples/elastic_cut3.rs
rename to examples3d/src/elastic_cut3.rs
index 98b4c1a..ccb0876 100644
--- a/crates/wgsparkl3d/examples/elastic_cut3.rs
+++ b/examples3d/src/elastic_cut3.rs
@@ -1,30 +1,21 @@
-use wgsparkl_testbed3d::{wgsparkl, RapierData};
+use wgsparkl_testbed3d::{wgsparkl, Callbacks, RapierData};
 
-use bevy::prelude::*;
 use bevy::render::renderer::RenderDevice;
-use nalgebra::{vector, DMatrix, Isometry3, Similarity3, Vector3};
+use nalgebra::{vector, DMatrix, Isometry3};
 use rapier3d::geometry::HeightField;
 use rapier3d::prelude::{ColliderBuilder, RigidBodyBuilder};
-use wgebra::GpuSim3;
-use wgparry3d::parry::shape::Cuboid;
-use wgrapier3d::dynamics::BodyDesc;
-use wgsparkl::models::DruckerPrager;
 use wgsparkl::{
     models::ElasticCoefficients,
     pipeline::MpmData,
     solver::{Particle, ParticleDynamics, ParticlePhase, SimulationParams},
 };
-use wgsparkl_testbed3d::{init_testbed, AppState, PhysicsContext, SceneInits};
-
-fn main() {
-    panic!("Run the `testbed3` example instead.");
-}
+use wgsparkl_testbed3d::{AppState, PhysicsContext};
 
 pub fn elastic_cut_demo(
-    mut commands: Commands,
-    device: Res<RenderDevice>,
-    mut app_state: ResMut<AppState>,
-) {
+    device: RenderDevice,
+    app_state: &mut AppState,
+    _callbacks: &mut Callbacks,
+) -> PhysicsContext {
     let mut rapier_data = RapierData::default();
     let device = device.wgpu_device();
 
@@ -51,6 +42,7 @@ pub fn elastic_cut_demo(
                         phase: 1.0,
                         max_stretch: f32::MAX,
                     }),
+                    color: None,
                 });
             }
         }
@@ -100,9 +92,9 @@ pub fn elastic_cut_demo(
         cell_width,
         60_000,
     );
-    commands.insert_resource(PhysicsContext {
+    PhysicsContext {
         data,
         rapier_data,
         particles,
-    });
+    }
 }
diff --git a/examples3d/src/glb_to_point_cloud_color.rs b/examples3d/src/glb_to_point_cloud_color.rs
new file mode 100644
index 0000000..5e03b8a
--- /dev/null
+++ b/examples3d/src/glb_to_point_cloud_color.rs
@@ -0,0 +1,67 @@
+use crate::utils::default_scene::{self, SAMPLE_PER_UNIT};
+
+use bevy::{math::Vec3, prelude::*, render::renderer::RenderDevice};
+use nalgebra::{vector, Isometry3, Transform3, UnitQuaternion, Vector3};
+use std::{fs::File, io::Read};
+use wgsparkl3d::load_mesh3d::load_gltf::load_model_with_colors;
+use wgsparkl3d::{pipeline::MpmData, solver::SimulationParams};
+use wgsparkl_testbed3d::{AppState, Callbacks, PhysicsContext, RapierData};
+
+pub fn elastic_color_model_demo(
+    device: RenderDevice,
+    app_state: &mut AppState,
+    _callbacks: &mut Callbacks,
+) -> PhysicsContext {
+    let mut file = File::open("assets/shiba.glb").expect("Failed to open GLB file");
+    let mut buffer = Vec::new();
+    file.read_to_end(&mut buffer).expect("Failed to read file");
+    let pc_grid = load_model_with_colors(
+        &buffer,
+        Transform3::from_matrix_unchecked(
+            Isometry3::from_parts(
+                Vector3::new(0.0, 6.0, 0.0).into(),
+                UnitQuaternion::identity(),
+            )
+            .to_matrix()
+                * nalgebra::Matrix4::new_scaling(3.0),
+        ),
+        None,
+        SAMPLE_PER_UNIT,
+    );
+    let params = SimulationParams {
+        gravity: vector![0.0, -9.81, 0.0] * app_state.gravity_factor,
+        dt: (1.0 / 60.0) / (app_state.num_substeps as f32),
+    };
+    let mut rapier_data = RapierData::default();
+    if !app_state.restarting {
+        app_state.num_substeps = 20;
+        app_state.gravity_factor = 1.0;
+    };
+    default_scene::spawn_ground_and_walls(&mut rapier_data);
+
+    let mut particles = vec![];
+    for (pos, color) in pc_grid {
+        let particle = default_scene::create_particle(
+            &Vec3::new(pos.x, pos.y, pos.z),
+            Some(Color::from(Srgba::from_u8_array(color))),
+            1f32 / SAMPLE_PER_UNIT / 2f32,
+        );
+        particles.push(particle);
+    }
+    default_scene::spawn_ground_and_walls(&mut rapier_data);
+    let data = MpmData::new(
+        device.wgpu_device(),
+        params,
+        &particles,
+        &rapier_data.bodies,
+        &rapier_data.colliders,
+        1f32 / default_scene::SAMPLE_PER_UNIT,
+        60_000,
+    );
+
+    PhysicsContext {
+        data,
+        rapier_data,
+        particles,
+    }
+}
diff --git a/crates/wgsparkl3d/examples/heightfield3.rs b/examples3d/src/heightfield3.rs
similarity index 79%
rename from crates/wgsparkl3d/examples/heightfield3.rs
rename to examples3d/src/heightfield3.rs
index a0281d5..152d1ca 100644
--- a/crates/wgsparkl3d/examples/heightfield3.rs
+++ b/examples3d/src/heightfield3.rs
@@ -1,30 +1,22 @@
-use wgsparkl_testbed3d::{wgsparkl, RapierData};
+use wgsparkl_testbed3d::{wgsparkl, Callbacks, RapierData};
 
-use bevy::prelude::*;
 use bevy::render::renderer::RenderDevice;
-use nalgebra::{vector, DMatrix, Similarity3, Vector3};
+use nalgebra::{vector, DMatrix};
 use rapier3d::geometry::HeightField;
 use rapier3d::prelude::{ColliderBuilder, RigidBodyBuilder};
-use wgebra::GpuSim3;
-use wgparry3d::parry::shape::Cuboid;
-use wgrapier3d::dynamics::BodyDesc;
 use wgsparkl::models::DruckerPrager;
 use wgsparkl::{
     models::ElasticCoefficients,
     pipeline::MpmData,
-    solver::{Particle, ParticleDynamics, ParticlePhase, SimulationParams},
+    solver::{Particle, ParticleDynamics, SimulationParams},
 };
-use wgsparkl_testbed3d::{init_testbed, AppState, PhysicsContext, SceneInits};
-
-fn main() {
-    panic!("Run the `testbed3` example instead.");
-}
+use wgsparkl_testbed3d::{AppState, PhysicsContext};
 
 pub fn heightfield_demo(
-    mut commands: Commands,
-    device: Res<RenderDevice>,
-    mut app_state: ResMut<AppState>,
-) {
+    device: RenderDevice,
+    app_state: &mut AppState,
+    _callbacks: &mut Callbacks,
+) -> PhysicsContext {
     let mut rapier_data = RapierData::default();
     let device = device.wgpu_device();
 
@@ -48,6 +40,7 @@ pub fn heightfield_demo(
                     model: ElasticCoefficients::from_young_modulus(2_000_000_000.0, 0.2),
                     plasticity: Some(DruckerPrager::new(2_000_000_000.0, 0.2)),
                     phase: None,
+                    color: None,
                 });
             }
         }
@@ -84,9 +77,9 @@ pub fn heightfield_demo(
         cell_width,
         60_000,
     );
-    commands.insert_resource(PhysicsContext {
+    PhysicsContext {
         data,
         rapier_data,
         particles,
-    });
+    }
 }
diff --git a/crates/wgsparkl3d/examples/testbed3.rs b/examples3d/src/main.rs
similarity index 54%
rename from crates/wgsparkl3d/examples/testbed3.rs
rename to examples3d/src/main.rs
index ec79fbf..09158a5 100644
--- a/crates/wgsparkl3d/examples/testbed3.rs
+++ b/examples3d/src/main.rs
@@ -1,17 +1,11 @@
-use wgsparkl_testbed3d::{wgsparkl, RapierData};
-
 use bevy::prelude::*;
-use wgparry3d::parry::shape::Cuboid;
-use wgrapier3d::dynamics::BodyDesc;
-use wgsparkl::models::DruckerPrager;
-use wgsparkl::{
-    models::ElasticCoefficients,
-    pipeline::MpmData,
-    solver::{Particle, ParticlePhase, SimulationParams},
-};
-use wgsparkl_testbed3d::{init_testbed, AppState, PhysicsContext, SceneInits};
+use wgsparkl_testbed3d::{init_testbed, SceneInitFn, SceneInits};
+
+pub mod utils;
 
+mod banana3;
 mod elastic_cut3;
+mod glb_to_point_cloud_color;
 mod heightfield3;
 mod sand3;
 
@@ -28,16 +22,21 @@ pub fn main() {
 }
 
 fn register_scenes(world: &mut World) {
-    let scenes = vec![
-        ("sand".to_string(), world.register_system(sand3::sand_demo)),
+    let scenes: Vec<(String, SceneInitFn)> = vec![
+        ("sand".to_string(), Box::new(sand3::sand_demo)),
         (
             "heightfield".to_string(),
-            world.register_system(heightfield3::heightfield_demo),
+            Box::new(heightfield3::heightfield_demo),
         ),
         (
             "elastic_cut".to_string(),
-            world.register_system(elastic_cut3::elastic_cut_demo),
+            Box::new(elastic_cut3::elastic_cut_demo),
+        ),
+        (
+            "elastic_model_colors".to_string(),
+            Box::new(glb_to_point_cloud_color::elastic_color_model_demo),
         ),
+        ("taichi_banana".to_string(), Box::new(banana3::demo)),
     ];
     let mut inits = world.resource_mut::<SceneInits>();
     inits.scenes = scenes;
diff --git a/examples3d/src/sand3.rs b/examples3d/src/sand3.rs
new file mode 100644
index 0000000..2c7570a
--- /dev/null
+++ b/examples3d/src/sand3.rs
@@ -0,0 +1,101 @@
+use wgsparkl_testbed3d::{wgsparkl, Callbacks, RapierData};
+
+use bevy::render::renderer::RenderDevice;
+use nalgebra::vector;
+use rapier3d::prelude::{ColliderBuilder, RigidBodyBuilder};
+use wgsparkl::models::DruckerPrager;
+use wgsparkl::{
+    models::ElasticCoefficients,
+    pipeline::MpmData,
+    solver::{Particle, ParticleDynamics, SimulationParams},
+};
+use wgsparkl_testbed3d::{AppState, PhysicsContext};
+
+pub fn sand_demo(
+    device: RenderDevice,
+    app_state: &mut AppState,
+    _callbacks: &mut Callbacks,
+) -> PhysicsContext {
+    let mut rapier_data = RapierData::default();
+    let device = device.wgpu_device();
+
+    let nxz = 45;
+    let cell_width = 1.0;
+    let mut particles = vec![];
+    for i in 0..nxz {
+        for j in 0..100 {
+            for k in 0..nxz {
+                let position = vector![
+                    i as f32 + 0.5 - nxz as f32 / 2.0,
+                    j as f32 + 0.5 + 10.0,
+                    k as f32 + 0.5 - nxz as f32 / 2.0
+                ] * cell_width
+                    / 2.0;
+                let density = 2700.0;
+                let radius = cell_width / 4.0;
+                particles.push(Particle {
+                    position,
+                    dynamics: ParticleDynamics::with_density(radius, density),
+                    model: ElasticCoefficients::from_young_modulus(2_000_000_000.0, 0.2),
+                    plasticity: Some(DruckerPrager::new(2_000_000_000.0, 0.2)),
+                    phase: None,
+                    color: None,
+                });
+            }
+        }
+    }
+
+    if !app_state.restarting {
+        app_state.num_substeps = 20;
+        app_state.gravity_factor = 1.0;
+    };
+
+    let params = SimulationParams {
+        gravity: vector![0.0, -9.81, 0.0] * app_state.gravity_factor,
+        dt: (1.0 / 60.0) / (app_state.num_substeps as f32),
+    };
+
+    rapier_data.insert_body_and_collider(
+        RigidBodyBuilder::fixed().translation(vector![0.0, -4.0, 0.0]),
+        ColliderBuilder::cuboid(100.0, 4.0, 100.0),
+    );
+
+    rapier_data.insert_body_and_collider(
+        RigidBodyBuilder::fixed().translation(vector![0.0, 5.0, -35.0]),
+        ColliderBuilder::cuboid(35.0, 5.0, 0.5),
+    );
+    rapier_data.insert_body_and_collider(
+        RigidBodyBuilder::fixed().translation(vector![0.0, 5.0, 35.0]),
+        ColliderBuilder::cuboid(35.0, 5.0, 0.5),
+    );
+    rapier_data.insert_body_and_collider(
+        RigidBodyBuilder::fixed().translation(vector![-35.0, 5.0, 0.0]),
+        ColliderBuilder::cuboid(0.5, 5.0, 35.0),
+    );
+    rapier_data.insert_body_and_collider(
+        RigidBodyBuilder::fixed().translation(vector![35.0, 5.0, 0.0]),
+        ColliderBuilder::cuboid(0.5, 5.0, 35.0),
+    );
+
+    let rb = RigidBodyBuilder::kinematic_velocity_based()
+        .translation(vector![0.0, 1.5, 0.0])
+        .rotation(vector![0.0, 0.0, -0.5])
+        .angvel(vector![0.0, -1.0, 0.0]);
+    let co = ColliderBuilder::cuboid(0.5, 1.5, 30.0);
+    let (_rb, _) = rapier_data.insert_body_and_collider(rb, co);
+
+    let data = MpmData::new(
+        device,
+        params,
+        &particles,
+        &rapier_data.bodies,
+        &rapier_data.colliders,
+        cell_width,
+        60_000,
+    );
+    PhysicsContext {
+        data,
+        rapier_data,
+        particles,
+    }
+}
diff --git a/examples3d/src/utils/default_scene.rs b/examples3d/src/utils/default_scene.rs
new file mode 100644
index 0000000..7add5de
--- /dev/null
+++ b/examples3d/src/utils/default_scene.rs
@@ -0,0 +1,52 @@
+#![allow(unused)]
+
+use bevy::prelude::*;
+use nalgebra::vector;
+use rapier3d::prelude::{ColliderBuilder, RigidBodyBuilder};
+use wgsparkl3d::{
+    models::ElasticCoefficients,
+    solver::{Particle, ParticleDynamics, ParticlePhase},
+};
+use wgsparkl_testbed3d::{AppState, RapierData};
+
+pub const SAMPLE_PER_UNIT: f32 = 10.0;
+
+/// Spawns a ground and 4 walls.
+pub fn spawn_ground_and_walls(rapier_data: &mut RapierData) {
+    rapier_data.insert_body_and_collider(
+        RigidBodyBuilder::fixed().translation(vector![0.0, -4.0, 0.0]),
+        ColliderBuilder::cuboid(100.0, 4.0, 100.0),
+    );
+
+    rapier_data.insert_body_and_collider(
+        RigidBodyBuilder::fixed().translation(vector![0.0, 5.0, -35.0]),
+        ColliderBuilder::cuboid(35.0, 5.0, 0.5),
+    );
+    rapier_data.insert_body_and_collider(
+        RigidBodyBuilder::fixed().translation(vector![0.0, 5.0, 35.0]),
+        ColliderBuilder::cuboid(35.0, 5.0, 0.5),
+    );
+    rapier_data.insert_body_and_collider(
+        RigidBodyBuilder::fixed().translation(vector![-35.0, 5.0, 0.0]),
+        ColliderBuilder::cuboid(0.5, 5.0, 35.0),
+    );
+    rapier_data.insert_body_and_collider(
+        RigidBodyBuilder::fixed().translation(vector![35.0, 5.0, 0.0]),
+        ColliderBuilder::cuboid(0.5, 5.0, 35.0),
+    );
+}
+
+pub fn create_particle(pos: &Vec3, color: Option<Color>, radius: f32) -> Particle {
+    let density = 3700.0;
+    Particle {
+        position: pos.to_array().into(),
+        dynamics: ParticleDynamics::with_density(radius, density),
+        model: ElasticCoefficients::from_young_modulus(10_000_000.0, 0.22),
+        plasticity: None,
+        phase: Some(ParticlePhase {
+            phase: 1f32,
+            max_stretch: f32::MAX,
+        }),
+        color: color.map(|c| c.to_linear().to_u8_array()),
+    }
+}
diff --git a/examples3d/src/utils/mod.rs b/examples3d/src/utils/mod.rs
new file mode 100644
index 0000000..d954233
--- /dev/null
+++ b/examples3d/src/utils/mod.rs
@@ -0,0 +1 @@
+pub mod default_scene;
diff --git a/notes.thierry.md b/notes.thierry.md
new file mode 100644
index 0000000..224b579
--- /dev/null
+++ b/notes.thierry.md
@@ -0,0 +1,6 @@
+# Notes
+
+banana example:
+
+- oneway+dynamic coupling makes the knife going places (bouncing weirdly).
+  - workaround: not use dynamic, it wasn't planned anyway.
diff --git a/src/lib.rs b/src/lib.rs
index 8ec9a0b..680d7f9 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -13,6 +13,8 @@ pub extern crate wgrapier3d as wgrapier;
 
 mod collision;
 pub mod grid;
+#[cfg(feature = "dim3")]
+pub mod load_mesh3d;
 pub mod models;
 pub mod pipeline;
 pub mod solver;
diff --git a/src/load_mesh3d/load_bevy.rs b/src/load_mesh3d/load_bevy.rs
new file mode 100644
index 0000000..5a14e87
--- /dev/null
+++ b/src/load_mesh3d/load_bevy.rs
@@ -0,0 +1,31 @@
+use bevy::render::mesh::{Indices, Mesh, VertexAttributeValues};
+use nalgebra::Point3;
+
+use f32 as Real;
+
+pub fn load_mesh_vertices_indices(mesh: &Mesh) -> Option<(Vec<Point3<Real>>, Vec<[u32; 3]>)> {
+    use rapier3d::na::point;
+
+    let vertices = mesh.attribute(Mesh::ATTRIBUTE_POSITION)?;
+    let indices = mesh.indices()?;
+
+    let vtx: Vec<_> = match vertices {
+        VertexAttributeValues::Float32(vtx) => {
+            Some(vtx.chunks(3).map(|v| point![v[0], v[1], v[2]]).collect())
+        }
+        VertexAttributeValues::Float32x3(vtx) => {
+            Some(vtx.iter().map(|v| point![v[0], v[1], v[2]]).collect())
+        }
+        _ => None,
+    }?;
+
+    let idx = match indices {
+        Indices::U16(idx) => idx
+            .chunks_exact(3)
+            .map(|i| [i[0] as u32, i[1] as u32, i[2] as u32])
+            .collect(),
+        Indices::U32(idx) => idx.chunks_exact(3).map(|i| [i[0], i[1], i[2]]).collect(),
+    };
+
+    Some((vtx, idx))
+}
diff --git a/src/load_mesh3d/load_gltf.rs b/src/load_mesh3d/load_gltf.rs
new file mode 100644
index 0000000..be00229
--- /dev/null
+++ b/src/load_mesh3d/load_gltf.rs
@@ -0,0 +1,282 @@
+use f32 as Real;
+use image::RgbaImage;
+use nalgebra::{distance_squared, Matrix4, Point3, Transform3};
+
+use super::get_point_cloud_from_trimesh;
+
+fn extract_embedded_texture(
+    gltf: &gltf::Document,
+    buffers: &[gltf::buffer::Data],
+) -> Option<RgbaImage> {
+    for image in gltf.images() {
+        dbg!(image.name());
+        if let gltf::image::Source::View { view, .. } = image.source() {
+            let buffer = &buffers[view.buffer().index()];
+            let start = view.offset();
+            let end = start + view.length();
+            let image_data = &buffer.0[start..end];
+
+            // Decode based on MIME type
+            if let Ok(img) = image::load_from_memory(image_data) {
+                return Some(img.to_rgba8());
+            } else {
+                eprintln!("Failed to decode texture.");
+            }
+        }
+    }
+    None
+}
+
+fn sample_texture(texture: &RgbaImage, uv: [f32; 2]) -> [u8; 4] {
+    let (width, height) = texture.dimensions();
+
+    // Convert UV (0.0 - 1.0) to pixel coordinates
+    let x = (uv[0] * width as f32).clamp(0.0, (width - 1) as f32) as u32;
+    let y = ((uv[1]) * height as f32).clamp(0.0, (height - 1) as f32) as u32;
+
+    texture.get_pixel(x, y).0
+}
+
+fn closest_point(
+    target: Point3<f32>,
+    points: &[(Point3<f32>, [u8; 4])],
+) -> Option<&(Point3<f32>, [u8; 4])> {
+    points.iter().min_by(|a, b| {
+        distance_squared(&a.0, &target)
+            .partial_cmp(&distance_squared(&b.0, &target))
+            .unwrap()
+    })
+}
+
+// TODO: transform should not be here, but when we spawn the model.
+// this function should have a "resolution" like the amount of points to sample in x/y/z.
+/// Load a glb or gltf file and return a point cloud and a list of trimeshes.
+///
+/// ```
+/// let mut file = File::open("path_to_file.glb").expect("Failed to open GLB file");
+/// let mut buffer = Vec::new();
+/// file.read_to_end(&mut buffer).expect("Failed to read file");
+///
+/// let _ = load_model_with_colors(&buffer);
+/// ```
+pub fn load_model_with_colors<S>(
+    slice: S,
+    transform: Transform3<Real>,
+    color_inside: Option<[u8; 4]>,
+    sample_per_unit: f32,
+) -> Vec<(Point3<Real>, [u8; 4])>
+where
+    S: AsRef<[u8]>,
+{
+    let mut res = load_model_with_vertex_colors(slice);
+    let mut pc_grid = vec![];
+
+    res.0
+        .iter_mut()
+        .for_each(|p| p.0 = transform.transform_point(&p.0));
+
+    for (t_id, trimesh) in res.1.iter_mut().enumerate() {
+        if t_id != 0 {
+            //continue;
+        }
+        trimesh.0.iter_mut().for_each(|p| {
+            let new_p = transform.transform_point(&Point3::new(p.x, p.y, p.z));
+            *p = Point3::new(new_p.x, new_p.y, new_p.z);
+        });
+        let mut pc = get_point_cloud_from_trimesh(&trimesh.0, &trimesh.1, sample_per_unit)
+            .into_iter()
+            .map(|p| {
+                let closest_color = closest_point(p, &res.0).unwrap();
+                if let Some(color_inside) = color_inside {
+                    let distance = distance_squared(&p, &closest_color.0);
+                    (
+                        p,
+                        if distance <= 0.04 {
+                            closest_color.1
+                        } else {
+                            color_inside
+                        },
+                    )
+                } else {
+                    (p, closest_color.1)
+                }
+            })
+            .collect();
+        pc_grid.append(&mut pc);
+    }
+    pc_grid
+}
+
+fn get_node_transform(node: &gltf::Node, parent_transform: Matrix4<f32>) -> Matrix4<f32> {
+    let matrix = node.transform().matrix();
+    parent_transform * Matrix4::from_column_slice(&matrix.concat())
+}
+
+pub struct SceneIterator<'a, T, F>
+where
+    F: Fn(&gltf::Primitive, &Matrix4<f32>) -> T + 'a,
+{
+    stack: Vec<(gltf::Node<'a>, Matrix4<f32>)>,
+    process_primitive: F,
+}
+
+impl<'a, T, F> SceneIterator<'a, T, F>
+where
+    F: Fn(&gltf::Primitive, &Matrix4<f32>) -> T + 'a,
+{
+    fn new(root_nodes: Vec<gltf::Node<'a>>, process_primitive: F) -> Self {
+        let stack = root_nodes
+            .into_iter()
+            .map(|node| (node, Matrix4::identity()))
+            .collect();
+        SceneIterator {
+            stack,
+            process_primitive,
+        }
+    }
+}
+
+impl<'a, T, F> Iterator for SceneIterator<'a, T, F>
+where
+    F: Fn(&gltf::Primitive, &Matrix4<f32>) -> T + 'a,
+{
+    type Item = T;
+
+    fn next(&mut self) -> Option<Self::Item> {
+        while let Some((node, parent_transform)) = self.stack.pop() {
+            let world_transform = get_node_transform(&node, parent_transform);
+
+            if let Some(mesh) = node.mesh() {
+                if let Some(primitive) = mesh.primitives().next() {
+                    return Some((self.process_primitive)(&primitive, &world_transform));
+                }
+            }
+
+            for child in node.children() {
+                self.stack.push((child, world_transform));
+            }
+        }
+        None
+    }
+}
+
+/// Load a glb or gltf file into vertices and indices and its vertex colors from sampled texture.
+///
+/// ```
+/// let mut file = File::open("path_to_file.glb").expect("Failed to open GLB file");
+/// let mut buffer = Vec::new();
+/// file.read_to_end(&mut buffer).expect("Failed to read file");
+///
+/// let (gltf, buffers, _) = gltf::import_slice(slice).expect("Failed to parse GLB");
+///
+/// let _ = load_model_with_point_cloud(&buffer);
+/// ```
+pub fn load_model_with_vertex_colors<S>(
+    slice: S,
+) -> (
+    Vec<(nalgebra::Point3<f32>, [u8; 4])>,
+    Vec<(Vec<nalgebra::Point3<f32>>, Vec<usize>)>,
+)
+where
+    S: AsRef<[u8]>,
+{
+    let (gltf, buffers, _) = gltf::import_slice(slice).expect("Failed to parse GLB");
+
+    // Extract embedded texture
+    let texture = extract_embedded_texture(&gltf, &buffers).expect("Failed to extract texture");
+
+    let mut pcs = vec![];
+    let mut trimeshes = vec![];
+
+    for scene in gltf.scenes() {
+        let iterator = SceneIterator::new(scene.nodes().collect(), |primitive, world_transform| {
+            let reader = primitive.reader(|buffer| Some(&*buffers[buffer.index()].0));
+
+            let mut point_cloud = vec![];
+            let mut positions = vec![];
+
+            if let (Some(positions_iter), Some(tex_coords)) =
+                (reader.read_positions(), reader.read_tex_coords(0))
+            {
+                for (pos, uv) in positions_iter.zip(tex_coords.into_f32()) {
+                    let color = sample_texture(&texture, uv);
+                    let position = Point3::new(pos[0], pos[1], pos[2]);
+                    let transformed = world_transform.transform_point(&position);
+                    point_cloud.push((transformed, color));
+                    positions.push(transformed);
+                }
+
+                let indices = reader
+                    .read_indices()
+                    .expect("No indices found")
+                    .into_u32()
+                    .map(|i| i as usize)
+                    .collect::<Vec<_>>();
+
+                (point_cloud, (positions, indices))
+            } else {
+                (vec![], (vec![], vec![]))
+            }
+        });
+
+        for (pc, trimesh) in iterator {
+            pcs.extend(pc);
+            trimeshes.push(trimesh);
+        }
+    }
+
+    (pcs, trimeshes)
+}
+
+/// Load a glb or gltf file and return a point cloud and a list of trimeshes.
+///
+/// ```
+/// let mut file = File::open("path_to_file.glb").expect("Failed to open GLB file");
+/// let mut buffer = Vec::new();
+/// file.read_to_end(&mut buffer).expect("Failed to read file");
+///
+/// let (gltf, buffers, _) = gltf::import_slice(slice).expect("Failed to parse GLB");
+///
+/// let _ = load_model_trimeshes(&buffer);
+/// ```
+pub fn load_model_trimeshes<S>(slice: S) -> Vec<(Vec<nalgebra::Point3<f32>>, Vec<usize>)>
+where
+    S: AsRef<[u8]>,
+{
+    let (gltf, buffers, _) = gltf::import_slice(slice).expect("Failed to parse GLB");
+
+    let mut trimeshes = vec![];
+
+    for scene in gltf.scenes() {
+        let iterator = SceneIterator::new(scene.nodes().collect(), |primitive, world_transform| {
+            let reader = primitive.reader(|buffer| Some(&*buffers[buffer.index()].0));
+
+            let mut positions = vec![];
+
+            if let Some(positions_iter) = reader.read_positions() {
+                for pos in positions_iter {
+                    let position = Point3::new(pos[0], pos[1], pos[2]);
+                    let transformed = world_transform.transform_point(&position);
+                    positions.push(transformed);
+                }
+
+                let indices = reader
+                    .read_indices()
+                    .expect("No indices found")
+                    .into_u32()
+                    .map(|i| i as usize)
+                    .collect::<Vec<_>>();
+
+                (positions, indices)
+            } else {
+                (vec![], vec![])
+            }
+        });
+
+        for trimesh in iterator {
+            trimeshes.push(trimesh);
+        }
+    }
+
+    trimeshes
+}
diff --git a/src/load_mesh3d/load_obj.rs b/src/load_mesh3d/load_obj.rs
new file mode 100644
index 0000000..992088c
--- /dev/null
+++ b/src/load_mesh3d/load_obj.rs
@@ -0,0 +1,29 @@
+use std::io::BufRead;
+
+use nalgebra::{point, Point3};
+use obj::{raw::object::Polygon, ObjResult};
+
+use super::recenter_and_scale;
+
+pub fn get_point_cloud<T: BufRead>(
+    input: T,
+    sample_per_unit: f32,
+) -> ObjResult<(Vec<Point3<f32>>, Vec<usize>)> {
+    let model = obj::raw::parse_obj(input)?;
+    let mut vertices: Vec<_> = model
+        .positions
+        .iter()
+        .map(|v| point![v.0, v.1, v.2])
+        .collect();
+    let indices: Vec<_> = model
+        .polygons
+        .into_iter()
+        .flat_map(|p| match p {
+            Polygon::P(idx) => idx,
+            Polygon::PT(idx) | Polygon::PN(idx) => idx.iter().map(|i| i.0).collect(),
+            Polygon::PTN(idx) => idx.iter().map(|i| i.0).collect(),
+        })
+        .collect();
+    recenter_and_scale(&mut vertices, sample_per_unit);
+    Ok((vertices, indices))
+}
diff --git a/src/load_mesh3d/mod.rs b/src/load_mesh3d/mod.rs
new file mode 100644
index 0000000..5930fce
--- /dev/null
+++ b/src/load_mesh3d/mod.rs
@@ -0,0 +1,62 @@
+#[cfg(feature = "load_bevy")]
+pub mod load_bevy;
+#[cfg(feature = "load_gltf")]
+pub mod load_gltf;
+#[cfg(feature = "load_obj")]
+pub mod load_obj;
+
+use nalgebra::{Point3, Vector3};
+
+use crate::rapier;
+use f32 as Real;
+use rapier::parry::bounding_volume;
+use rapier::prelude::{PointQuery, TriMesh, TriMeshFlags};
+
+pub fn recenter_and_scale(
+    vertices: &mut Vec<nalgebra::OPoint<f32, nalgebra::Const<3>>>,
+    sample_per_unit: f32,
+) {
+    // Compute the size of the model, to scale it and have similar size for everything.
+    let aabb =
+        bounding_volume::details::point_cloud_aabb(&rapier::na::Isometry::default(), &*vertices);
+    let center = aabb.center();
+    let diag = (aabb.maxs - aabb.mins).norm();
+    vertices
+        .iter_mut()
+        .for_each(|p| *p = (*p - center.coords) * sample_per_unit / diag);
+}
+
+pub fn get_point_cloud_from_trimesh(
+    vertices: &[nalgebra::OPoint<f32, nalgebra::Const<3>>],
+    indices: &[usize],
+    sample_per_unit: f32,
+) -> Vec<Point3<Real>> {
+    let indices: Vec<_> = indices
+        .chunks(3)
+        .map(|idx| [idx[0] as u32, idx[1] as u32, idx[2] as u32])
+        .collect();
+    let trimesh = TriMesh::with_flags(vertices.to_vec(), indices, TriMeshFlags::ORIENTED)
+        .expect("Invalid mesh");
+    let aabb = bounding_volume::details::point_cloud_aabb(
+        &rapier3d::na::Isometry::default(),
+        trimesh.vertices(),
+    );
+    let mut positions = vec![];
+
+    let aabb_sample = aabb.scaled(&Vector3::new(
+        sample_per_unit,
+        sample_per_unit,
+        sample_per_unit,
+    ));
+    for x in aabb_sample.mins.x as i32..aabb_sample.maxs.x as i32 {
+        for y in aabb_sample.mins.y as i32..aabb_sample.maxs.y as i32 {
+            for z in aabb_sample.mins.z as i32..aabb_sample.maxs.z as i32 {
+                let point = Point3::new(x as f32, y as f32, z as f32) / sample_per_unit;
+                if trimesh.contains_local_point(&point) {
+                    positions.push(point);
+                }
+            }
+        }
+    }
+    positions
+}
diff --git a/src/solver/particle2d.rs b/src/solver/particle2d.rs
index ac3bb81..6501e25 100644
--- a/src/solver/particle2d.rs
+++ b/src/solver/particle2d.rs
@@ -55,6 +55,7 @@ pub struct Particle {
     pub model: ElasticCoefficients,
     pub plasticity: Option<DruckerPrager>,
     pub phase: Option<ParticlePhase>,
+    pub color: Option<[u8; 4]>,
 }
 
 pub struct GpuRigidParticles {
diff --git a/src/solver/particle3d.rs b/src/solver/particle3d.rs
index 18bbe46..565a2e6 100644
--- a/src/solver/particle3d.rs
+++ b/src/solver/particle3d.rs
@@ -56,6 +56,7 @@ pub struct Particle {
     pub model: ElasticCoefficients,
     pub plasticity: Option<DruckerPrager>,
     pub phase: Option<ParticlePhase>,
+    pub color: Option<[u8; 4]>,
 }
 
 #[derive(Copy, Clone, Debug, ShaderType)]
diff --git a/src_testbed/lib.rs b/src_testbed/lib.rs
index 3c87f68..985448c 100644
--- a/src_testbed/lib.rs
+++ b/src_testbed/lib.rs
@@ -3,6 +3,7 @@ pub extern crate wgsparkl2d as wgsparkl;
 #[cfg(feature = "dim3")]
 pub extern crate wgsparkl3d as wgsparkl;
 
+use bevy::render::renderer::RenderDevice;
 #[cfg(feature = "dim2")]
 pub use instancing2d as instancing;
 #[cfg(feature = "dim3")]
@@ -22,7 +23,7 @@ pub mod step;
 pub mod ui;
 
 use bevy::asset::load_internal_asset;
-use bevy::ecs::system::SystemId;
+use bevy::ecs::system::{BoxedSystem, SystemId};
 use bevy::pbr::wireframe::WireframePlugin;
 use bevy::prelude::*;
 use bevy_editor_cam::prelude::DefaultEditorCamPlugins;
@@ -55,12 +56,15 @@ pub fn init_testbed(app: &mut App) {
         .add_plugins(instancing::ParticlesMaterialPlugin)
         .add_plugins(bevy_egui::EguiPlugin)
         .init_resource::<SceneInits>()
+        .init_resource::<Callbacks>()
         .add_systems(Startup, startup::setup_app)
         .add_systems(
             Update,
             (
                 ui::update_ui,
-                step::step_simulation,
+                (step::step_simulation, step::callbacks)
+                    .chain()
+                    .run_if(|state: Res<AppState>| state.run_state != RunState::Paused),
                 rigid_graphics::update_rigid_graphics,
                 hot_reload::handle_hot_reloading,
             )
@@ -99,21 +103,10 @@ pub struct AppState {
     pub selected_scene: usize,
     pub hot_reload: HotReloadState,
     pub show_rigid_particles: bool,
+    pub physics_time_seconds: f64,
 }
 
-#[derive(Default)]
-pub struct RapierData {
-    pub bodies: RigidBodySet,
-    pub colliders: ColliderSet,
-    pub impulse_joints: ImpulseJointSet,
-    pub multibody_joints: MultibodyJointSet,
-    pub params: IntegrationParameters,
-    pub physics_pipeline: PhysicsPipeline,
-    pub narrow_phase: NarrowPhase,
-    pub broad_phase: DefaultBroadPhase,
-    pub ccd_solver: CCDSolver,
-    pub islands: IslandManager,
-}
+pub use wgsparkl::rapier::prelude::PhysicsContext as RapierData;
 
 #[derive(Resource)]
 pub struct PhysicsContext {
@@ -122,6 +115,15 @@ pub struct PhysicsContext {
     pub particles: Vec<Particle>,
 }
 
+#[derive(Resource, Default)]
+pub struct Callbacks(pub Vec<Callback>);
+
+pub type Callback = Box<
+    dyn FnMut(Option<&mut RenderContext>, &mut PhysicsContext, &Timestamps, &AppState)
+        + Send
+        + Sync,
+>;
+
 #[derive(Resource, Default)]
 pub struct RenderContext {
     pub instanced_materials: InstancedMaterials,
@@ -168,22 +170,42 @@ pub enum RunState {
 
 #[derive(Resource)]
 pub struct SceneInits {
-    pub scenes: Vec<(String, SystemId)>,
+    pub scenes: Vec<(String, SceneInitFn)>,
     reset_graphics: SystemId,
+    reset_app_state: SystemId,
 }
+pub type SceneInitFn =
+    Box<dyn FnMut(RenderDevice, &mut AppState, &mut Callbacks) -> PhysicsContext + Send + Sync>;
 
 impl SceneInits {
     pub fn init_scene(&self, commands: &mut Commands, scene_id: usize) {
-        commands.run_system(self.scenes[scene_id].1);
+        commands.run_system(self.reset_app_state);
+        commands.run_system_cached_with(run_scene_init, scene_id);
         commands.run_system(self.reset_graphics);
     }
 }
 
+pub fn run_scene_init(
+    scene_id: In<usize>,
+    mut commands: Commands,
+    mut scenes: ResMut<SceneInits>,
+    device: Res<RenderDevice>,
+    mut app_state: ResMut<AppState>,
+    mut callbacks: ResMut<Callbacks>,
+) {
+    commands.insert_resource(scenes.scenes[scene_id.0].1(
+        device.clone(),
+        &mut app_state,
+        &mut callbacks,
+    ));
+}
+
 impl FromWorld for SceneInits {
     fn from_world(world: &mut World) -> Self {
         Self {
             scenes: vec![],
             reset_graphics: world.register_system(startup::setup_graphics),
+            reset_app_state: world.register_system(startup::setup_app_state),
         }
     }
 }
diff --git a/src_testbed/startup.rs b/src_testbed/startup.rs
index 0f8cc3c..0a32a71 100644
--- a/src_testbed/startup.rs
+++ b/src_testbed/startup.rs
@@ -2,7 +2,7 @@ use crate::instancing::{InstanceBuffer, InstanceData, InstanceMaterialData};
 use crate::prep_vertex_buffer::{GpuRenderConfig, RenderConfig, RenderMode, WgPrepVertexBuffer};
 use crate::rigid_graphics::{BevyMaterial, EntityWithGraphics, InstancedMaterials};
 use crate::step::TimestampChannel;
-use crate::{AppState, PhysicsContext, RenderContext, RunState, Timestamps};
+use crate::{AppState, Callbacks, PhysicsContext, RenderContext, RunState, Timestamps};
 use bevy::asset::Assets;
 use bevy::color::Color;
 use bevy::hierarchy::DespawnRecursiveExt;
@@ -47,6 +47,7 @@ pub fn setup_app(mut commands: Commands, device: Res<RenderDevice>) {
         selected_scene: 0,
         hot_reload,
         show_rigid_particles: false,
+        physics_time_seconds: 0.0,
     });
     commands.init_resource::<RenderContext>();
 
@@ -176,7 +177,10 @@ fn setup_particles_graphics(
      */
     let mut instances = vec![];
     for (rb_id, particle) in physics.particles.iter().enumerate() {
-        let base_color = colors[rb_id % colors.len()].to_linear().to_f32_array();
+        let base_color = particle.color.map_or_else(
+            || colors[rb_id % colors.len()].to_linear().to_f32_array(),
+            |c| c.map(|c| c as f32 / 255f32),
+        );
         instances.push(InstanceData {
             deformation: [Vec4::X, Vec4::Y, Vec4::Z],
             #[cfg(feature = "dim2")]
@@ -285,3 +289,8 @@ pub fn setup_graphics(
         &to_clear,
     );
 }
+
+pub fn setup_app_state(mut app_state: ResMut<AppState>, mut callbacks: ResMut<Callbacks>) {
+    app_state.physics_time_seconds = 0.0;
+    callbacks.0.clear();
+}
diff --git a/src_testbed/step.rs b/src_testbed/step.rs
index 5a8481a..7da5d7f 100644
--- a/src_testbed/step.rs
+++ b/src_testbed/step.rs
@@ -1,19 +1,17 @@
 use crate::instancing::InstanceMaterialData;
 use crate::startup::RigidParticlesTag;
-use crate::{AppState, PhysicsContext, RunState, Timestamps};
+use crate::{AppState, Callbacks, PhysicsContext, RenderContext, RunState, Timestamps};
 use async_channel::{Receiver, Sender};
 use bevy::prelude::*;
-use bevy::render::renderer::{RenderDevice, RenderQueue, WgpuWrapper};
+use bevy::render::renderer::{RenderDevice, RenderQueue};
 use bevy::tasks::ComputeTaskPool;
-use std::time::Instant;
 use wgcore::kernel::KernelInvocationQueue;
 use wgcore::re_exports::encase::StorageBuffer;
 use wgcore::timestamps::GpuTimestamps;
-use wgpu::{Device, Queue};
 use wgsparkl::rapier::math::Vector;
 use wgsparkl::rapier::prelude::RigidBodyPosition;
 use wgsparkl::wgparry::math::GpuSim;
-use wgsparkl::wgrapier::dynamics::{GpuBodySet, GpuVelocity};
+use wgsparkl::wgrapier::dynamics::GpuVelocity;
 
 #[derive(Resource)]
 pub struct TimestampChannel {
@@ -21,6 +19,23 @@ pub struct TimestampChannel {
     pub rcv: Receiver<Timestamps>,
 }
 
+pub fn callbacks(
+    mut render: ResMut<RenderContext>,
+    mut physics: ResMut<PhysicsContext>,
+    app_state: ResMut<AppState>,
+    timings: Res<Timestamps>,
+    mut callbacks: ResMut<Callbacks>,
+) {
+    for to_call in callbacks.0.iter_mut() {
+        to_call(
+            Some(&mut render),
+            &mut physics,
+            timings.as_ref(),
+            &app_state,
+        );
+    }
+}
+
 pub fn step_simulation(
     mut timings: ResMut<Timestamps>,
     render_device: Res<RenderDevice>,
@@ -55,10 +70,6 @@ pub fn step_simulation_legacy(
     rigid_particles: &Query<&InstanceMaterialData, With<RigidParticlesTag>>,
     timings_channel: &TimestampChannel,
 ) {
-    if app_state.run_state == RunState::Paused {
-        return;
-    }
-
     let timings = &mut *timings;
 
     while let Ok(new_timings) = timings_channel.rcv.try_recv() {
@@ -105,7 +116,9 @@ pub fn step_simulation_legacy(
             let rb = &physics.rapier_data.bodies[coupling.body];
             GpuVelocity {
                 linear: *rb.linvel()
-                    + gravity * physics.rapier_data.params.dt * (rb.is_dynamic() as u32 as f32)
+                    + gravity
+                        * physics.rapier_data.integration_parameters.dt
+                        * (rb.is_dynamic() as u32 as f32)
                         / (app_state.num_substeps as f32),
                 angular: rb.angvel().clone(),
             }
@@ -186,7 +199,7 @@ pub fn step_simulation_legacy(
                     next_position: new_poses[i].isometry,
                 };
                 let vel = interpolator.interpolate_velocity(
-                    1.0 / (physics.rapier_data.params.dt / divisor),
+                    1.0 / (physics.rapier_data.integration_parameters.dt / divisor),
                     &rb.mass_properties().local_mprops.local_com,
                 );
                 rb.set_linvel(vel.linvel, true);
@@ -196,12 +209,12 @@ pub fn step_simulation_legacy(
         }
     }
 
-    let mut params = physics.rapier_data.params;
+    let mut params = physics.rapier_data.integration_parameters;
     params.dt = params.dt / divisor;
     physics.rapier_data.physics_pipeline.step(
         &nalgebra::zero(),
         &params, // physics.rapier_data.params,
-        &mut physics.rapier_data.islands,
+        &mut physics.rapier_data.island_manager,
         &mut physics.rapier_data.broad_phase,
         &mut physics.rapier_data.narrow_phase,
         &mut physics.rapier_data.bodies,
@@ -213,6 +226,7 @@ pub fn step_simulation_legacy(
         &(),
         &(),
     );
+    app_state.physics_time_seconds += params.dt as f64;
 
     if let Some(timestamps) = std::mem::take(&mut timings.timestamps) {
         let timings_snd = timings_channel.snd.clone();
diff --git a/src_testbed/ui.rs b/src_testbed/ui.rs
index 695392b..c454015 100644
--- a/src_testbed/ui.rs
+++ b/src_testbed/ui.rs
@@ -8,6 +8,7 @@ use bevy_egui::{egui, EguiContexts};
 use nalgebra::vector;
 use wgsparkl::solver::SimulationParams;
 
+#[allow(clippy::too_many_arguments)]
 pub fn update_ui(
     mut commands: Commands,
     mut ui_context: EguiContexts,