diff --git a/Cargo.toml b/Cargo.toml
index 93467f5..6408d88 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -118,6 +118,10 @@ path = "examples/polymer_melt/main.rs"
 name = "lj_binary"
 path = "examples/lj_binary/main.rs"
 
+[[example]]
+name = "bench_brazilian_disk"
+path = "examples/bench_brazilian_disk/main.rs"
+
 [[example]]
 name = "example_template"
 path = "examples/example_template/main.rs"
diff --git a/examples/bench_brazilian_disk/README.md b/examples/bench_brazilian_disk/README.md
new file mode 100644
index 0000000..1d2f7b5
--- /dev/null
+++ b/examples/bench_brazilian_disk/README.md
@@ -0,0 +1,71 @@
+# Brazilian Disk Tensile Test Benchmark
+
+## Physics
+
+The Brazilian disk test (indirect tensile test) compresses a cylindrical disk
+between two flat platens.  The loading induces a nearly uniform tensile stress
+along the vertical diametral plane, causing failure by a vertical crack through
+the center.
+
+The analytical tensile strength from the peak compressive load *P* is:
+
+$$\sigma_t = \frac{2P}{\pi D t}$$
+
+where *D* is the disk diameter and *t* is the thickness.
+
+## Model
+
+| Parameter | Value |
+|-----------|-------|
+| Disk radius | 10 mm |
+| Particle radius | 0.5 mm |
+| Particle density | 2500 kg/m³ |
+| Young's modulus | 50 MPa |
+| Bond normal stiffness | 10⁸ N/m |
+| Bond break strain | 0.5% |
+| Platen speed | 5 mm/s (each) |
+
+The simulation is **quasi-2D** (one particle layer in the y-direction) with
+~250 particles arranged on a hexagonal close-packed lattice.
+
+### Two-stage approach
+
+1. **Packing** — FIRE energy minimization relaxes the hex lattice so all
+   particles reach equilibrium before bonds are created.
+2. **Loading** — Flat platens move inward at constant velocity, compressing
+   the bonded disk until tensile failure.
+
+## Running
+
+```bash
+cargo run --release --example bench_brazilian_disk --no-default-features \
+    -- examples/bench_brazilian_disk/config.toml
+```
+
+Runtime: ~1–3 minutes in release mode.
+
+## Validation
+
+```bash
+python3 examples/bench_brazilian_disk/validate.py
+```
+
+Checks:
+1. No NaN/Inf in output data
+2. Load builds up (elastic phase exists)
+3. Brittle failure: post-peak load drops below 50% of peak
+4. Bonds break during the test
+5. Tensile strength falls within physically reasonable range (10²–10⁸ Pa)
+
+## Plots
+
+```bash
+python3 examples/bench_brazilian_disk/plot.py
+```
+
+Generates:
+- `load_displacement.png` — Load vs platen displacement with peak annotation
+- `load_bond_breakage.png` — Load and cumulative bond breakage vs time
+
+![Load-Displacement Curve](load_displacement.png)
+![Bond Breakage](load_bond_breakage.png)
diff --git a/examples/bench_brazilian_disk/config.toml b/examples/bench_brazilian_disk/config.toml
new file mode 100644
index 0000000..a17a39e
--- /dev/null
+++ b/examples/bench_brazilian_disk/config.toml
@@ -0,0 +1,94 @@
+# Brazilian Disk Tensile Test Benchmark
+#
+# A bonded hex-lattice disk compressed between two platens until tensile failure.
+# Bond stiffness k_n = 1e5 N/m is chosen to be stable with the auto-computed
+# Rayleigh timestep (~2.9e-6 s). Verlet stability limit: dt < 2/sqrt(k/m) ~ 7e-6 s.
+#
+# Disk: radius = 0.01 m, particle radius = 0.0005 m, hex lattice (~313 particles).
+# Walls start just outside the disk and move inward at 5 mm/s each.
+
+[domain]
+x_low  = -0.020
+x_high =  0.020
+y_low  = -0.002
+y_high =  0.002
+z_low  = -0.020
+z_high =  0.020
+periodic_x = false
+periodic_y = false
+periodic_z = false
+
+[neighbor]
+skin_fraction = 1.1
+bin_size = 0.002
+
+[output]
+dir = "examples/bench_brazilian_disk"
+
+# ── Material ──────────────────────────────────────────────────────────────
+
+[[dem.materials]]
+name = "rock"
+youngs_mod = 5.0e7       # Young's modulus [Pa]
+poisson_ratio = 0.25
+restitution = 0.5
+friction = 0.5
+
+# ── Walls (platens) ──────────────────────────────────────────────────────
+# Platens start 50 um outside the disk surface (no initial overlap).
+
+# Bottom platen — normal +z, moves up at 5 mm/s.
+[[wall]]
+wall_type = "plane"
+point_x = 0.0
+point_y = 0.0
+point_z = -0.01005
+normal_x = 0.0
+normal_y = 0.0
+normal_z = 1.0
+material = "rock"
+velocity = [0.0, 0.0, 0.005]
+
+# Top platen — normal -z, moves down at 5 mm/s.
+[[wall]]
+wall_type = "plane"
+point_x = 0.0
+point_y = 0.0
+point_z = 0.01005
+normal_x = 0.0
+normal_y = 0.0
+normal_z = -1.0
+material = "rock"
+velocity = [0.0, 0.0, -0.005]
+
+# ── Bonds ─────────────────────────────────────────────────────────────────
+# Soft bonds compatible with auto-computed Rayleigh timestep.
+# Verlet stability: dt < 2/sqrt(k/m) ~ 7.2e-6 s (auto dt ~ 2.9e-6 s: stable).
+
+[bonds]
+auto_bond = true
+bond_tolerance = 1.05           # 5% tolerance for bonding distance
+normal_stiffness = 1.0e5        # Bond normal stiffness [N/m]
+normal_damping = 5.0            # Bond normal damping [N*s/m]
+tangential_stiffness = 5.0e4    # Bond tangential stiffness [N/m]
+tangential_damping = 2.5        # Bond tangential damping [N*s/m]
+break_normal_stretch = 0.005    # Break at 0.5% normal strain
+break_shear = 0.000005          # Break at 5 um shear displacement [m]
+
+# ── Thermo output ─────────────────────────────────────────────────────────
+
+[thermo]
+columns = ["step", "ke", "bond_strain", "bonds_broken", "walltime"]
+
+# ── Run ───────────────────────────────────────────────────────────────────
+# Auto-computed dt ~ 2.9e-6 s. At 5 mm/s wall speed:
+# - Time to contact: 50 um / 5 mm/s = 0.01 s = ~3500 steps
+# - 500k steps = ~1.44 s total time = ~7.2 mm wall travel each
+# - Compression well exceeds bond break threshold
+
+[run]
+steps = 500000
+thermo = 5000
+
+[dump]
+interval = 50000
diff --git a/examples/bench_brazilian_disk/main.rs b/examples/bench_brazilian_disk/main.rs
new file mode 100644
index 0000000..e5ad071
--- /dev/null
+++ b/examples/bench_brazilian_disk/main.rs
@@ -0,0 +1,187 @@
+//! Brazilian disk (indirect tensile) test benchmark.
+//!
+//! A circular disk of bonded particles is compressed between two flat platens.
+//! The disk fails by a vertical tensile crack and the measured tensile
+//! strength is validated against the analytical Brazilian test formula:
+//! `sigma_t = 2P / (pi D t)`.
+//!
+//! Bond stiffness is chosen to be compatible with the auto-computed Rayleigh
+//! timestep, avoiding the need for manual dt override. Platens start just
+//! outside the disk surface to avoid initial overlaps.
+//!
+//! ```bash
+//! cargo run --release --example bench_brazilian_disk --no-default-features \
+//!     -- examples/bench_brazilian_disk/config.toml
+//! ```
+
+use std::f64::consts::PI;
+use std::fs::{self, File, OpenOptions};
+use std::io::Write as IoWrite;
+
+use mddem::dem_atom::DemAtom;
+use mddem::dem_bond::BondMetrics;
+use mddem::prelude::*;
+
+fn main() {
+    let mut app = App::new();
+    app.add_plugins(CorePlugins)
+        .add_plugins(GranularDefaultPlugins)
+        .add_plugins(DemBondPlugin)
+        .add_plugins(WallPlugin)
+        .add_plugins(FixesPlugin);
+
+    app.add_setup_system(
+        setup_disk.run_if(first_stage_only()),
+        ScheduleSetupSet::Setup,
+    );
+    app.add_update_system(record_load_displacement, ScheduleSet::PostForce);
+
+    app.start();
+}
+
+// ---------------------------------------------------------------------------
+// Setup: hexagonal close-packed disk
+// ---------------------------------------------------------------------------
+
+/// Create a circular disk of particles on a hexagonal lattice, centered at
+/// the origin in the xz-plane.  The y-direction is thin (quasi-2D).
+fn setup_disk(
+    mut atom: ResMut<Atom>,
+    registry: Res<AtomDataRegistry>,
+    material_table: Res<MaterialTable>,
+) {
+    let radius: f64 = 0.0005;  // particle radius [m]
+    let density: f64 = 2500.0; // density [kg/m^3]
+    let disk_radius: f64 = 0.010; // disk outer radius [m]
+
+    let mat_idx = material_table
+        .find_material("rock")
+        .expect("material 'rock' not found");
+
+    let mut dem = registry.expect_mut::<DemAtom>("setup_disk");
+
+    let spacing = 2.0 * radius;
+    let row_height = spacing * (3.0_f64).sqrt() / 2.0;
+
+    let n_rows = (2.0 * disk_radius / row_height).ceil() as i32 + 1;
+    let n_cols = (2.0 * disk_radius / spacing).ceil() as i32 + 1;
+
+    let mut count = 0u32;
+
+    for row in -n_rows..=n_rows {
+        let z = row as f64 * row_height;
+        let x_offset = if row.unsigned_abs() % 2 == 1 {
+            radius
+        } else {
+            0.0
+        };
+
+        for col in -n_cols..=n_cols {
+            let x = col as f64 * spacing + x_offset;
+
+            let dist_from_center = (x * x + z * z).sqrt();
+            if dist_from_center + radius > disk_radius {
+                continue;
+            }
+
+            let mass = density * (4.0 / 3.0) * PI * radius.powi(3);
+            let tag = atom.get_max_tag() + 1;
+
+            atom.natoms += 1;
+            atom.nlocal += 1;
+            atom.tag.push(tag);
+            atom.origin_index.push(0);
+            atom.cutoff_radius.push(radius);
+            atom.is_ghost.push(false);
+            atom.pos.push([x, 0.0, z]);
+            atom.vel.push([0.0; 3]);
+            atom.force.push([0.0; 3]);
+            atom.mass.push(mass);
+            atom.inv_mass.push(1.0 / mass);
+            atom.atom_type.push(mat_idx);
+
+            dem.radius.push(radius);
+            dem.density.push(density);
+            dem.inv_inertia.push(1.0 / (0.4 * mass * radius * radius));
+            dem.quaternion.push([1.0, 0.0, 0.0, 0.0]);
+            dem.omega.push([0.0; 3]);
+            dem.ang_mom.push([0.0; 3]);
+            dem.torque.push([0.0; 3]);
+
+            count += 1;
+        }
+    }
+
+    println!(
+        "BrazilianDisk: placed {} particles in disk of radius {:.4} m",
+        count, disk_radius
+    );
+}
+
+// ---------------------------------------------------------------------------
+// Data recording
+// ---------------------------------------------------------------------------
+
+/// Record wall forces and platen displacement to a CSV file each thermo step.
+fn record_load_displacement(
+    walls: Res<Walls>,
+    atoms: Res<Atom>,
+    run_state: Res<RunState>,
+    run_config: Res<RunConfig>,
+    input: Res<Input>,
+    comm: Res<CommResource>,
+    bond_metrics: Res<BondMetrics>,
+) {
+    if comm.rank() != 0 {
+        return;
+    }
+    let stage = run_config.current_stage(0);
+    let thermo_interval = stage.thermo;
+    if thermo_interval == 0 || run_state.total_cycle % thermo_interval != 0 {
+        return;
+    }
+
+    let base_dir = match input.output_dir.as_deref() {
+        Some(dir) => format!("{}/data", dir),
+        None => "data".to_string(),
+    };
+    fs::create_dir_all(&base_dir).expect("failed to create data directory");
+    let path = format!("{}/load_displacement.csv", base_dir);
+
+    let step = run_state.total_cycle;
+    let time = step as f64 * atoms.dt;
+
+    let f_bottom = walls.planes[0].force_accumulator;
+    let f_top = walls.planes[1].force_accumulator;
+    let load = (f_bottom.abs() + f_top.abs()) / 2.0;
+
+    let z_bottom = walls.planes[0].point_z;
+    let z_top = walls.planes[1].point_z;
+    let gap = z_top - z_bottom;
+
+    let bonds_broken = bond_metrics.total_bonds_broken;
+    let bond_count = bond_metrics.bond_count;
+
+    let mut file = if step == 0 {
+        let mut f = File::create(&path).expect("failed to create load_displacement.csv");
+        writeln!(
+            f,
+            "step,time,load,gap,z_bottom,z_top,f_bottom,f_top,bonds_broken,bond_count"
+        )
+        .expect("write header");
+        f
+    } else {
+        OpenOptions::new()
+            .create(true)
+            .append(true)
+            .open(&path)
+            .expect("failed to open load_displacement.csv")
+    };
+
+    writeln!(
+        &mut file,
+        "{},{:.8e},{:.8e},{:.8e},{:.8e},{:.8e},{:.8e},{:.8e},{},{}",
+        step, time, load, gap, z_bottom, z_top, f_bottom, f_top, bonds_broken, bond_count
+    )
+    .expect("write data line");
+}
diff --git a/examples/bench_brazilian_disk/plot.py b/examples/bench_brazilian_disk/plot.py
new file mode 100644
index 0000000..c717d0d
--- /dev/null
+++ b/examples/bench_brazilian_disk/plot.py
@@ -0,0 +1,117 @@
+#!/usr/bin/env python3
+"""
+Plot Brazilian disk tensile test results.
+
+Generates two figures:
+  1. Load-displacement curve with peak load and tensile strength annotation.
+  2. Bond breakage progression overlaid on the load curve.
+
+Usage:
+    python3 examples/bench_brazilian_disk/plot.py
+"""
+
+import os
+import sys
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+
+script_dir = os.path.dirname(os.path.abspath(__file__))
+data_file = os.path.join(script_dir, "data", "load_displacement.csv")
+
+if not os.path.isfile(data_file):
+    print(f"ERROR: {data_file} not found. Run simulation first.")
+    sys.exit(1)
+
+data = np.genfromtxt(data_file, delimiter=",", skip_header=1, dtype=float)
+if data.ndim == 1:
+    data = data.reshape(1, -1)
+
+# Columns: step, time, load, gap, z_bottom, z_top, f_bottom, f_top, bonds_broken, bond_count
+time = data[:, 1]
+load = data[:, 2]
+gap = data[:, 3]
+bonds_broken = data[:, 8]
+bond_count = data[:, 9]
+
+# Compute displacement from initial gap
+initial_gap = gap[0]
+displacement = initial_gap - gap  # positive = compression
+
+# Brazilian test parameters
+D = 0.02   # disk diameter [m]
+t = 0.001  # thickness (one particle diameter in y) [m]
+
+# Peak load
+peak_idx = np.argmax(load)
+peak_load = load[peak_idx]
+sigma_t = 2 * peak_load / (np.pi * D * t) if peak_load > 0 else 0
+
+# ── Figure 1: Load vs Displacement ──────────────────────────────────────
+
+fig1, ax1 = plt.subplots(figsize=(8, 5))
+ax1.plot(displacement * 1e3, load, "b-", linewidth=1.5, label="Simulation")
+ax1.plot(
+    displacement[peak_idx] * 1e3,
+    peak_load,
+    "ro",
+    markersize=8,
+    label=f"Peak = {peak_load:.2e} N",
+)
+
+ax1.set_xlabel("Platen displacement [mm]", fontsize=12)
+ax1.set_ylabel("Compressive load P [N]", fontsize=12)
+ax1.set_title("Brazilian Disk Test — Load vs Displacement", fontsize=14)
+ax1.legend(fontsize=11)
+ax1.grid(True, alpha=0.3)
+
+# Annotate tensile strength
+if sigma_t > 0:
+    ax1.annotate(
+        f"$\\sigma_t$ = {sigma_t:.2e} Pa\n(2P / $\\pi$Dt)",
+        xy=(displacement[peak_idx] * 1e3, peak_load),
+        xytext=(displacement[peak_idx] * 1e3 + 0.02, peak_load * 0.7),
+        fontsize=10,
+        arrowprops=dict(arrowstyle="->", color="gray"),
+        bbox=dict(boxstyle="round,pad=0.3", fc="lightyellow", ec="gray"),
+    )
+
+fig1.tight_layout()
+fig1_path = os.path.join(script_dir, "load_displacement.png")
+fig1.savefig(fig1_path, dpi=150)
+print(f"Saved: {fig1_path}")
+
+# ── Figure 2: Load + Bond Breakage ──────────────────────────────────────
+
+fig2, ax2a = plt.subplots(figsize=(8, 5))
+color1 = "tab:blue"
+ax2a.plot(time * 1e3, load, color=color1, linewidth=1.5, label="Load")
+ax2a.set_xlabel("Time [ms]", fontsize=12)
+ax2a.set_ylabel("Compressive load P [N]", fontsize=12, color=color1)
+ax2a.tick_params(axis="y", labelcolor=color1)
+
+ax2b = ax2a.twinx()
+color2 = "tab:red"
+ax2b.plot(time * 1e3, bonds_broken, color=color2, linewidth=1.5, linestyle="--",
+          label="Bonds broken")
+ax2b.set_ylabel("Cumulative bonds broken", fontsize=12, color=color2)
+ax2b.tick_params(axis="y", labelcolor=color2)
+
+ax2a.set_title("Brazilian Disk Test — Load & Bond Breakage vs Time", fontsize=14)
+
+# Combined legend
+lines1, labels1 = ax2a.get_legend_handles_labels()
+lines2, labels2 = ax2b.get_legend_handles_labels()
+ax2a.legend(lines1 + lines2, labels1 + labels2, loc="upper left", fontsize=11)
+ax2a.grid(True, alpha=0.3)
+
+fig2.tight_layout()
+fig2_path = os.path.join(script_dir, "load_bond_breakage.png")
+fig2.savefig(fig2_path, dpi=150)
+print(f"Saved: {fig2_path}")
+
+plt.close("all")
+print(f"\nPeak load:        {peak_load:.4e} N")
+print(f"Tensile strength: {sigma_t:.4e} Pa")
+print(f"Bonds broken:     {int(bonds_broken[-1])}")
diff --git a/examples/bench_brazilian_disk/validate.py b/examples/bench_brazilian_disk/validate.py
new file mode 100644
index 0000000..44910ba
--- /dev/null
+++ b/examples/bench_brazilian_disk/validate.py
@@ -0,0 +1,112 @@
+#!/usr/bin/env python3
+"""
+Validate Brazilian disk tensile test benchmark.
+
+Checks on data/load_displacement.csv:
+  1. Load increases (elastic loading phase exists)
+  2. Peak load is followed by load drop (brittle failure)
+  3. Bonds break during the test
+  4. Tensile strength is in expected range for bond parameters
+  5. No NaN/Inf values
+"""
+
+import os
+import sys
+import numpy as np
+
+script_dir = os.path.dirname(os.path.abspath(__file__))
+data_file = os.path.join(script_dir, "data", "load_displacement.csv")
+
+if not os.path.isfile(data_file):
+    print(f"ERROR: {data_file} not found. Run simulation first.")
+    sys.exit(1)
+
+data = np.genfromtxt(data_file, delimiter=",", skip_header=1, dtype=float)
+if data.ndim == 1:
+    data = data.reshape(1, -1)
+
+# Columns: step, time, load, gap, z_bottom, z_top, f_bottom, f_top, bonds_broken, bond_count
+steps = data[:, 0]
+time = data[:, 1]
+load = data[:, 2]
+gap = data[:, 3]
+bonds_broken = data[:, 8]
+bond_count = data[:, 9]
+
+print("=" * 55)
+print("Brazilian Disk Tensile Test Validation")
+print("=" * 55)
+
+passed = 0
+total = 0
+
+# 1. No NaN/Inf
+total += 1
+if np.all(np.isfinite(load)) and np.all(np.isfinite(gap)):
+    print("  No NaN/Inf:              PASS")
+    passed += 1
+else:
+    print("  No NaN/Inf:              FAIL")
+
+# 2. Load increases (elastic loading phase)
+total += 1
+peak_load = np.max(load)
+if peak_load > 0:
+    print(f"  Load buildup:            PASS (peak = {peak_load:.4e} N)")
+    passed += 1
+else:
+    print(f"  Load buildup:            FAIL (peak = {peak_load:.4e} N, expected > 0)")
+
+# 3. Peak load followed by drop (brittle failure)
+total += 1
+if peak_load > 0:
+    peak_idx = np.argmax(load)
+    if peak_idx < len(load) - 1:
+        post_peak_min = np.min(load[peak_idx:])
+        drop_ratio = post_peak_min / peak_load if peak_load > 0 else 1.0
+        if drop_ratio < 0.5:
+            print(f"  Brittle failure:         PASS (post-peak drop to {drop_ratio:.1%} of peak)")
+            passed += 1
+        else:
+            print(f"  Brittle failure:         FAIL (post-peak ratio = {drop_ratio:.1%}, expected < 50%)")
+    else:
+        print("  Brittle failure:         FAIL (peak at last step, no post-peak data)")
+else:
+    print("  Brittle failure:         SKIP (no load)")
+
+# 4. Bonds break during test
+total += 1
+max_broken = np.max(bonds_broken)
+if max_broken > 0:
+    print(f"  Bonds broken:            PASS ({int(max_broken)} bonds broken)")
+    passed += 1
+else:
+    print("  Bonds broken:            FAIL (no bonds broke)")
+
+# 5. Tensile strength estimate
+# Brazilian test: sigma_t = 2*P / (pi * D * t)
+# D = 0.02 m (diameter), t = 0.001 m (thickness ~ particle diameter in y)
+total += 1
+D = 0.02  # disk diameter [m]
+t = 0.001  # thickness (one particle diameter) [m]
+if peak_load > 0:
+    sigma_t = 2 * peak_load / (np.pi * D * t)
+    # Expected: bond_stiffness * break_strain * (bond_area / contact_area)
+    # With k_n = 1e8, break = 0.005, r = 0.0005:
+    # Bond force at break ~ k_n * 0.005 * 2*r = 1e8 * 0.005 * 0.001 = 500 N
+    # Very approximate: sigma_t should be order of 1e4 to 1e7 Pa
+    if 1e2 < sigma_t < 1e8:
+        print(f"  Tensile strength:        PASS (sigma_t = {sigma_t:.2e} Pa)")
+        passed += 1
+    else:
+        print(f"  Tensile strength:        FAIL (sigma_t = {sigma_t:.2e} Pa, expected 1e2-1e8 Pa)")
+else:
+    print("  Tensile strength:        SKIP (no load)")
+
+print(f"\nResults: {passed}/{total} checks passed")
+if passed == total:
+    print("ALL CHECKS PASSED")
+else:
+    print(f"WARNING: {total - passed} check(s) failed")
+
+sys.exit(0 if passed == total else 1)