dirt_wall 0.1.4

General plane wall contact forces with Hertz repulsion for DIRT simulations
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129

//! Plane-wall motion modes and update systems.

use crate::geometry::Walls;
use grass_scheduler::prelude::*;
use soil_core::Atom;

/// Wall motion mode at runtime (plane walls only).
pub enum WallMotion {
    /// Wall is stationary.
    Static,
    /// Wall translates at a constant velocity each timestep.
    ConstantVelocity,
    /// Wall oscillates sinusoidally: `displacement = amplitude * sin(2π * frequency * t)`.
    Oscillate {
        /// Peak displacement (meters).
        amplitude: f64,
        /// Oscillation frequency (Hz).
        frequency: f64,
    },
    /// Proportional servo: velocity = clamp(gain × (target_force − measured_force), ±max_velocity).
    Servo {
        /// Desired total contact force (N).
        target_force: f64,
        /// Maximum wall speed (m/s).
        max_velocity: f64,
        /// Proportional gain (m/s per N).
        gain: f64,
    },
}

/// Runtime representation of an infinite plane wall.
///
/// Contact geometry for a plane wall:
/// ```text
///          ·  particle center (px, py, pz)
///         /|
///        / |  distance = dot(P - W, n̂)
///       /  |
///   n̂ ↑   |  delta (overlap) = radius - distance
///      |   |
///  ────W───┴──── wall surface
///      (point_x, point_y, point_z)
/// ```
///
/// The signed distance from the particle center to the wall plane is
/// `dot(pos - point, normal)`. Overlap is `delta = radius - distance`.

// ── Systems ─────────────────────────────────────────────────────────────────

/// Update wall positions and velocities according to their motion mode.
///
/// Runs in [`ParticleSimScheduleSet::PreInitialIntegration`] so walls are moved
/// *before* the integration step each timestep. Advances `walls.time` by `dt`.
pub fn wall_move(mut walls: ResMut<Walls>, atoms: Res<Atom>) {
    let dt = atoms.dt;
    let time = walls.time;

    let nplanes = walls.planes.len();
    for idx in 0..nplanes {
        if !walls.active[idx] {
            continue;
        }

        let wall = &mut walls.planes[idx];
        match wall.motion {
            WallMotion::Static => {}
            WallMotion::ConstantVelocity => {
                wall.point_x += wall.velocity[0] * dt;
                wall.point_y += wall.velocity[1] * dt;
                wall.point_z += wall.velocity[2] * dt;
            }
            WallMotion::Oscillate {
                amplitude,
                frequency,
            } => {
                let phase = 2.0 * std::f64::consts::PI * frequency * (time + dt);
                let disp = amplitude * phase.sin();
                wall.point_x = wall.origin[0] + disp * wall.normal_x;
                wall.point_y = wall.origin[1] + disp * wall.normal_y;
                wall.point_z = wall.origin[2] + disp * wall.normal_z;
                // Velocity = d(disp)/dt = amplitude * 2*pi*freq * cos(phase)
                let vel_mag = amplitude * 2.0 * std::f64::consts::PI * frequency * phase.cos();
                wall.velocity = [
                    vel_mag * wall.normal_x,
                    vel_mag * wall.normal_y,
                    vel_mag * wall.normal_z,
                ];
            }
            WallMotion::Servo {
                target_force,
                max_velocity,
                gain,
            } => {
                let error = target_force - wall.force_accumulator;
                let vel_mag = (gain * error).clamp(-max_velocity, max_velocity);
                wall.velocity = [
                    vel_mag * wall.normal_x,
                    vel_mag * wall.normal_y,
                    vel_mag * wall.normal_z,
                ];
                wall.point_x += wall.velocity[0] * dt;
                wall.point_y += wall.velocity[1] * dt;
                wall.point_z += wall.velocity[2] * dt;
            }
        }
    }

    walls.time += dt;
}

/// Zero all per-wall force accumulators before the force computation pass.
///
/// Runs in [`ParticleSimScheduleSet::PreForce`]. The accumulators are summed during
/// [`wall_contact_force`] and read by servo controllers in the next
/// [`wall_move`] call.
pub fn wall_zero_force_accumulators(mut walls: ResMut<Walls>) {
    for wall in &mut walls.planes {
        wall.force_accumulator = 0.0;
    }
    for wall in &mut walls.cylinders {
        wall.force_accumulator = 0.0;
    }
    for wall in &mut walls.spheres {
        wall.force_accumulator = 0.0;
    }
    for wall in &mut walls.regions {
        wall.force_accumulator = 0.0;
    }
}