|
1 | | -use crate::math::Real; |
2 | | -use na::RealField; |
3 | | - |
4 | 1 | #[cfg(doc)] |
5 | 2 | use super::RigidBodyActivation; |
| 3 | +use crate::math::Real; |
| 4 | +use simba::simd::SimdRealField; |
6 | 5 |
|
7 | 6 | // TODO: enabling the block solver in 3d introduces a lot of jitters in |
8 | 7 | // the 3D domino demo. So for now we dont enable it in 3D. |
@@ -30,6 +29,114 @@ pub enum FrictionModel { |
30 | 29 | Coulomb, |
31 | 30 | } |
32 | 31 |
|
| 32 | +#[derive(Copy, Clone, Debug, PartialEq)] |
| 33 | +#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] |
| 34 | +// TODO: we should be able to combine this with MotorModel. |
| 35 | +/// Coefficients for a spring, typically used for configuring constraint softness for contacts and |
| 36 | +/// joints. |
| 37 | +pub struct SpringCoefficients<N> { |
| 38 | + /// Sets the natural frequency (Hz) of the spring-like constraint. |
| 39 | + /// |
| 40 | + /// Higher values make the constraint stiffer and resolve constraint violations more quickly. |
| 41 | + pub natural_frequency: N, |
| 42 | + /// Sets the damping ratio for the spring-like constraint. |
| 43 | + /// |
| 44 | + /// Larger values make the joint more compliant (allowing more drift before stabilization). |
| 45 | + pub damping_ratio: N, |
| 46 | +} |
| 47 | + |
| 48 | +impl<N: SimdRealField<Element = Real> + Copy> SpringCoefficients<N> { |
| 49 | + /// Initializes spring coefficients from the spring natural frequency and damping ratio. |
| 50 | + pub fn new(natural_frequency: N, damping_ratio: N) -> Self { |
| 51 | + Self { |
| 52 | + natural_frequency, |
| 53 | + damping_ratio, |
| 54 | + } |
| 55 | + } |
| 56 | + |
| 57 | + /// Default softness coefficients for contacts. |
| 58 | + pub fn contact_defaults() -> Self { |
| 59 | + Self { |
| 60 | + natural_frequency: N::splat(30.0), |
| 61 | + damping_ratio: N::splat(5.0), |
| 62 | + } |
| 63 | + } |
| 64 | + |
| 65 | + /// Default softness coefficients for joints. |
| 66 | + pub fn joint_defaults() -> Self { |
| 67 | + Self { |
| 68 | + natural_frequency: N::splat(1.0e6), |
| 69 | + damping_ratio: N::splat(1.0), |
| 70 | + } |
| 71 | + } |
| 72 | + |
| 73 | + /// The contact’s spring angular frequency for constraints regularization. |
| 74 | + pub fn angular_frequency(&self) -> N { |
| 75 | + self.natural_frequency * N::simd_two_pi() |
| 76 | + } |
| 77 | + |
| 78 | + /// The [`Self::erp`] coefficient, multiplied by the inverse timestep length. |
| 79 | + pub fn erp_inv_dt(&self, dt: N) -> N { |
| 80 | + let ang_freq = self.angular_frequency(); |
| 81 | + ang_freq / (dt * ang_freq + N::splat(2.0) * self.damping_ratio) |
| 82 | + } |
| 83 | + |
| 84 | + /// The effective Error Reduction Parameter applied for calculating regularization forces. |
| 85 | + /// |
| 86 | + /// This parameter is computed automatically from [`Self::natural_frequency`], |
| 87 | + /// [`Self::damping_ratio`] and the substep length. |
| 88 | + pub fn erp(&self, dt: N) -> N { |
| 89 | + dt * self.erp_inv_dt(dt) |
| 90 | + } |
| 91 | + |
| 92 | + /// Compute CFM assuming a critically damped spring multiplied by the damping ratio. |
| 93 | + /// |
| 94 | + /// This coefficient softens the impulse applied at each solver iteration. |
| 95 | + pub fn cfm_coeff(&self, dt: N) -> N { |
| 96 | + let one = N::one(); |
| 97 | + let erp = self.erp(dt); |
| 98 | + let erp_is_not_zero = erp.simd_ne(N::zero()); |
| 99 | + let inv_erp_minus_one = one / erp - one; |
| 100 | + |
| 101 | + // let stiffness = 4.0 * damping_ratio * damping_ratio * projected_mass |
| 102 | + // / (dt * dt * inv_erp_minus_one * inv_erp_minus_one); |
| 103 | + // let damping = 4.0 * damping_ratio * damping_ratio * projected_mass |
| 104 | + // / (dt * inv_erp_minus_one); |
| 105 | + // let cfm = 1.0 / (dt * dt * stiffness + dt * damping); |
| 106 | + // NOTE: This simplifies to cfm = cfm_coeff / projected_mass: |
| 107 | + let result = inv_erp_minus_one * inv_erp_minus_one |
| 108 | + / ((one + inv_erp_minus_one) * N::splat(4.0) * self.damping_ratio * self.damping_ratio); |
| 109 | + result.select(erp_is_not_zero, N::zero()) |
| 110 | + } |
| 111 | + |
| 112 | + /// The CFM factor to be used in the constraint resolution. |
| 113 | + /// |
| 114 | + /// This parameter is computed automatically from [`Self::natural_frequency`], |
| 115 | + /// [`Self::damping_ratio`] and the substep length. |
| 116 | + pub fn cfm_factor(&self, dt: N) -> N { |
| 117 | + let one = N::one(); |
| 118 | + let cfm_coeff = self.cfm_coeff(dt); |
| 119 | + |
| 120 | + // We use this coefficient inside the impulse resolution. |
| 121 | + // Surprisingly, several simplifications happen there. |
| 122 | + // Let `m` the projected mass of the constraint. |
| 123 | + // Let `m’` the projected mass that includes CFM: `m’ = 1 / (1 / m + cfm_coeff / m) = m / (1 + cfm_coeff)` |
| 124 | + // We have: |
| 125 | + // new_impulse = old_impulse - m’ (delta_vel - cfm * old_impulse) |
| 126 | + // = old_impulse - m / (1 + cfm_coeff) * (delta_vel - cfm_coeff / m * old_impulse) |
| 127 | + // = old_impulse * (1 - cfm_coeff / (1 + cfm_coeff)) - m / (1 + cfm_coeff) * delta_vel |
| 128 | + // = old_impulse / (1 + cfm_coeff) - m * delta_vel / (1 + cfm_coeff) |
| 129 | + // = 1 / (1 + cfm_coeff) * (old_impulse - m * delta_vel) |
| 130 | + // So, setting cfm_factor = 1 / (1 + cfm_coeff). |
| 131 | + // We obtain: |
| 132 | + // new_impulse = cfm_factor * (old_impulse - m * delta_vel) |
| 133 | + // |
| 134 | + // The value returned by this function is this cfm_factor that can be used directly |
| 135 | + // in the constraint solver. |
| 136 | + one / (one + cfm_coeff) |
| 137 | + } |
| 138 | +} |
| 139 | + |
33 | 140 | /// Configuration parameters that control the physics simulation quality and behavior. |
34 | 141 | /// |
35 | 142 | /// These parameters affect how the physics engine advances time, resolves collisions, and |
@@ -80,34 +187,8 @@ pub struct IntegrationParameters { |
80 | 187 | /// to numerical instabilities. |
81 | 188 | pub min_ccd_dt: Real, |
82 | 189 |
|
83 | | - /// > 0: the damping ratio used by the springs for contact constraint stabilization. |
84 | | - /// |
85 | | - /// Larger values make the constraints more compliant (allowing more visible |
86 | | - /// penetrations before stabilization). |
87 | | - /// (default `5.0`). |
88 | | - pub contact_damping_ratio: Real, |
89 | | - |
90 | | - /// > 0: the natural frequency used by the springs for contact constraint regularization. |
91 | | - /// |
92 | | - /// Increasing this value will make it so that penetrations get fixed more quickly at the |
93 | | - /// expense of potential jitter effects due to overshooting. In order to make the simulation |
94 | | - /// look stiffer, it is recommended to increase the [`Self::contact_damping_ratio`] instead of this |
95 | | - /// value. |
96 | | - /// (default: `30.0`). |
97 | | - pub contact_natural_frequency: Real, |
98 | | - |
99 | | - /// > 0: the natural frequency used by the springs for joint constraint regularization. |
100 | | - /// |
101 | | - /// Increasing this value will make it so that penetrations get fixed more quickly. |
102 | | - /// (default: `1.0e6`). |
103 | | - pub joint_natural_frequency: Real, |
104 | | - |
105 | | - /// The fraction of critical damping applied to the joint for constraints regularization. |
106 | | - /// |
107 | | - /// Larger values make the constraints more compliant (allowing more joint |
108 | | - /// drift before stabilization). |
109 | | - /// (default `1.0`). |
110 | | - pub joint_damping_ratio: Real, |
| 190 | + /// Softness coefficients for contact constraints. |
| 191 | + pub contact_softness: SpringCoefficients<Real>, |
111 | 192 |
|
112 | 193 | /// The coefficient in `[0, 1]` applied to warmstart impulses, i.e., impulses that are used as the |
113 | 194 | /// initial solution (instead of 0) at the next simulation step. |
@@ -193,110 +274,7 @@ impl IntegrationParameters { |
193 | 274 | } |
194 | 275 | } |
195 | 276 |
|
196 | | - /// The contact’s spring angular frequency for constraints regularization. |
197 | | - pub fn contact_angular_frequency(&self) -> Real { |
198 | | - self.contact_natural_frequency * Real::two_pi() |
199 | | - } |
200 | | - |
201 | | - /// The [`Self::contact_erp`] coefficient, multiplied by the inverse timestep length. |
202 | | - pub fn contact_erp_inv_dt(&self) -> Real { |
203 | | - let ang_freq = self.contact_angular_frequency(); |
204 | | - ang_freq / (self.dt * ang_freq + 2.0 * self.contact_damping_ratio) |
205 | | - } |
206 | | - |
207 | | - /// The effective Error Reduction Parameter applied for calculating regularization forces |
208 | | - /// on contacts. |
209 | | - /// |
210 | | - /// This parameter is computed automatically from [`Self::contact_natural_frequency`], |
211 | | - /// [`Self::contact_damping_ratio`] and the substep length. |
212 | | - pub fn contact_erp(&self) -> Real { |
213 | | - self.dt * self.contact_erp_inv_dt() |
214 | | - } |
215 | | - |
216 | | - /// The joint’s spring angular frequency for constraint regularization. |
217 | | - pub fn joint_angular_frequency(&self) -> Real { |
218 | | - self.joint_natural_frequency * Real::two_pi() |
219 | | - } |
220 | | - |
221 | | - /// The [`Self::joint_erp`] coefficient, multiplied by the inverse timestep length. |
222 | | - pub fn joint_erp_inv_dt(&self) -> Real { |
223 | | - let ang_freq = self.joint_angular_frequency(); |
224 | | - ang_freq / (self.dt * ang_freq + 2.0 * self.joint_damping_ratio) |
225 | | - } |
226 | | - |
227 | | - /// The effective Error Reduction Parameter applied for calculating regularization forces |
228 | | - /// on joints. |
229 | | - /// |
230 | | - /// This parameter is computed automatically from [`Self::joint_natural_frequency`], |
231 | | - /// [`Self::joint_damping_ratio`] and the substep length. |
232 | | - pub fn joint_erp(&self) -> Real { |
233 | | - self.dt * self.joint_erp_inv_dt() |
234 | | - } |
235 | | - |
236 | | - /// The CFM factor to be used in the constraint resolution. |
237 | | - /// |
238 | | - /// This parameter is computed automatically from [`Self::contact_natural_frequency`], |
239 | | - /// [`Self::contact_damping_ratio`] and the substep length. |
240 | | - pub fn contact_cfm_factor(&self) -> Real { |
241 | | - // Compute CFM assuming a critically damped spring multiplied by the damping ratio. |
242 | | - // The logic is similar to [`Self::joint_cfm_coeff`]. |
243 | | - let contact_erp = self.contact_erp(); |
244 | | - if contact_erp == 0.0 { |
245 | | - return 0.0; |
246 | | - } |
247 | | - let inv_erp_minus_one = 1.0 / contact_erp - 1.0; |
248 | | - |
249 | | - // let stiffness = 4.0 * damping_ratio * damping_ratio * projected_mass |
250 | | - // / (dt * dt * inv_erp_minus_one * inv_erp_minus_one); |
251 | | - // let damping = 4.0 * damping_ratio * damping_ratio * projected_mass |
252 | | - // / (dt * inv_erp_minus_one); |
253 | | - // let cfm = 1.0 / (dt * dt * stiffness + dt * damping); |
254 | | - // NOTE: This simplifies to cfm = cfm_coeff / projected_mass: |
255 | | - let cfm_coeff = inv_erp_minus_one * inv_erp_minus_one |
256 | | - / ((1.0 + inv_erp_minus_one) |
257 | | - * 4.0 |
258 | | - * self.contact_damping_ratio |
259 | | - * self.contact_damping_ratio); |
260 | | - |
261 | | - // Furthermore, we use this coefficient inside of the impulse resolution. |
262 | | - // Surprisingly, several simplifications happen there. |
263 | | - // Let `m` the projected mass of the constraint. |
264 | | - // Let `m’` the projected mass that includes CFM: `m’ = 1 / (1 / m + cfm_coeff / m) = m / (1 + cfm_coeff)` |
265 | | - // We have: |
266 | | - // new_impulse = old_impulse - m’ (delta_vel - cfm * old_impulse) |
267 | | - // = old_impulse - m / (1 + cfm_coeff) * (delta_vel - cfm_coeff / m * old_impulse) |
268 | | - // = old_impulse * (1 - cfm_coeff / (1 + cfm_coeff)) - m / (1 + cfm_coeff) * delta_vel |
269 | | - // = old_impulse / (1 + cfm_coeff) - m * delta_vel / (1 + cfm_coeff) |
270 | | - // = 1 / (1 + cfm_coeff) * (old_impulse - m * delta_vel) |
271 | | - // So, setting cfm_factor = 1 / (1 + cfm_coeff). |
272 | | - // We obtain: |
273 | | - // new_impulse = cfm_factor * (old_impulse - m * delta_vel) |
274 | | - // |
275 | | - // The value returned by this function is this cfm_factor that can be used directly |
276 | | - // in the constraint solver. |
277 | | - 1.0 / (1.0 + cfm_coeff) |
278 | | - } |
279 | | - |
280 | | - /// The CFM (constraints force mixing) coefficient applied to all joints for constraints regularization. |
281 | | - /// |
282 | | - /// This parameter is computed automatically from [`Self::joint_natural_frequency`], |
283 | | - /// [`Self::joint_damping_ratio`] and the substep length. |
284 | | - pub fn joint_cfm_coeff(&self) -> Real { |
285 | | - // Compute CFM assuming a critically damped spring multiplied by the damping ratio. |
286 | | - // The logic is similar to `Self::contact_cfm_factor`. |
287 | | - let joint_erp = self.joint_erp(); |
288 | | - if joint_erp == 0.0 { |
289 | | - return 0.0; |
290 | | - } |
291 | | - let inv_erp_minus_one = 1.0 / joint_erp - 1.0; |
292 | | - inv_erp_minus_one * inv_erp_minus_one |
293 | | - / ((1.0 + inv_erp_minus_one) |
294 | | - * 4.0 |
295 | | - * self.joint_damping_ratio |
296 | | - * self.joint_damping_ratio) |
297 | | - } |
298 | | - |
299 | | - /// Amount of penetration the engine won’t attempt to correct (default: `0.001` multiplied by |
| 277 | + /// Amount of penetration the engine won't attempt to correct (default: `0.001` multiplied by |
300 | 278 | /// [`Self::length_unit`]). |
301 | 279 | pub fn allowed_linear_error(&self) -> Real { |
302 | 280 | self.normalized_allowed_linear_error * self.length_unit |
@@ -326,10 +304,7 @@ impl Default for IntegrationParameters { |
326 | 304 | Self { |
327 | 305 | dt: 1.0 / 60.0, |
328 | 306 | min_ccd_dt: 1.0 / 60.0 / 100.0, |
329 | | - contact_natural_frequency: 30.0, |
330 | | - contact_damping_ratio: 5.0, |
331 | | - joint_natural_frequency: 1.0e6, |
332 | | - joint_damping_ratio: 1.0, |
| 307 | + contact_softness: SpringCoefficients::contact_defaults(), |
333 | 308 | warmstart_coefficient: 1.0, |
334 | 309 | num_internal_pgs_iterations: 1, |
335 | 310 | num_internal_stabilization_iterations: 1, |
|
0 commit comments