1use crate::collision::ContactManifold;
10use crate::components::{RigidBody, Transform, Velocity};
11use gizmo_math::{Mat3, Vec3};
12
13pub struct ConstraintSolver {
18 pub iterations: usize,
20 pub baumgarte: f32,
23 pub slop: f32,
25 pub warm_start_factor: f32,
27 pub restitution_velocity_threshold: f32,
29 pub max_linear_correction: f32,
31 pub split_impulse_enabled: bool,
35 pub split_impulse_erp: f32,
37}
38
39impl Default for ConstraintSolver {
40 fn default() -> Self {
41 Self {
42 iterations: 20,
43 baumgarte: 0.15,
44 slop: 0.005,
45 warm_start_factor: 0.85,
46 restitution_velocity_threshold: 1.0,
47 max_linear_correction: 0.02,
48 split_impulse_enabled: true,
49 split_impulse_erp: 0.1,
50 }
51 }
52}
53
54impl ConstraintSolver {
55 pub fn new(iterations: usize) -> Self {
56 Self {
57 iterations,
58 ..Default::default()
59 }
60 }
61
62 pub fn solve_contacts(
67 &self,
68 manifolds: &mut [ContactManifold],
69 rigid_bodies: &[RigidBody],
70 transforms: &[Transform],
71 velocities: &mut [Velocity],
72 entity_index_map: &std::collections::HashMap<u32, usize>,
73 dt: f32,
74 ) {
75 if manifolds.is_empty() {
76 return;
77 }
78
79 let mut pseudo_vel: Vec<(Vec3, Vec3)> = vec![(Vec3::ZERO, Vec3::ZERO); velocities.len()];
83
84 for mid in 0..manifolds.len() {
86 let entity_a_id = manifolds[mid].entity_a.id();
87 let entity_b_id = manifolds[mid].entity_b.id();
88
89 let idx_a = match entity_index_map.get(&entity_a_id) {
90 Some(&i) => i,
91 None => continue,
92 };
93 let idx_b = match entity_index_map.get(&entity_b_id) {
94 Some(&i) => i,
95 None => continue,
96 };
97
98 let inv_m_a = rigid_bodies[idx_a].inv_mass();
99 let inv_m_b = rigid_bodies[idx_b].inv_mass();
100 let inv_i_a = rigid_bodies[idx_a].inv_world_inertia_tensor(transforms[idx_a].rotation);
101 let inv_i_b = rigid_bodies[idx_b].inv_world_inertia_tensor(transforms[idx_b].rotation);
102 let dyn_a = rigid_bodies[idx_a].is_dynamic();
103 let dyn_b = rigid_bodies[idx_b].is_dynamic();
104
105 let com_a = transforms[idx_a].position
106 + transforms[idx_a]
107 .rotation
108 .mul_vec3(rigid_bodies[idx_a].center_of_mass);
109 let com_b = transforms[idx_b].position
110 + transforms[idx_b]
111 .rotation
112 .mul_vec3(rigid_bodies[idx_b].center_of_mass);
113
114 for contact in &manifolds[mid].contacts {
115 let r_a = contact.point - com_a;
116 let r_b = contact.point - com_b;
117
118 let wn = contact.normal * (contact.normal_impulse * self.warm_start_factor);
119 let wt = contact.tangent_impulse * self.warm_start_factor;
120
121 if dyn_a {
122 velocities[idx_a].linear -= wn * inv_m_a;
123 velocities[idx_a].linear -= wt * inv_m_a;
124 velocities[idx_a].angular -= inv_i_a * (r_a.cross(wn) + r_a.cross(wt));
125 }
126 if dyn_b {
127 velocities[idx_b].linear += wn * inv_m_b;
128 velocities[idx_b].linear += wt * inv_m_b;
129 velocities[idx_b].angular += inv_i_b * (r_b.cross(wn) + r_b.cross(wt));
130 }
131 }
132 }
133
134 let inv_dt = if dt > 0.0 { 1.0 / dt } else { 0.0 };
136
137 for _ in 0..self.iterations {
138 for mid in 0..manifolds.len() {
139 let entity_a_id = manifolds[mid].entity_a.id();
140 let entity_b_id = manifolds[mid].entity_b.id();
141
142 let idx_a = match entity_index_map.get(&entity_a_id) {
143 Some(&i) => i,
144 None => continue,
145 };
146 let idx_b = match entity_index_map.get(&entity_b_id) {
147 Some(&i) => i,
148 None => continue,
149 };
150
151 let friction = manifolds[mid].friction;
152 let restitution = manifolds[mid].restitution;
153
154 for cid in 0..manifolds[mid].contacts.len() {
155 let contact_pt = manifolds[mid].contacts[cid].point;
156 let normal = manifolds[mid].contacts[cid].normal;
157 let penetration = manifolds[mid].contacts[cid].penetration;
158 let acc_n = manifolds[mid].contacts[cid].normal_impulse;
159 let acc_t = manifolds[mid].contacts[cid].tangent_impulse;
160
161 let com_a = transforms[idx_a].position
162 + transforms[idx_a]
163 .rotation
164 .mul_vec3(rigid_bodies[idx_a].center_of_mass);
165 let com_b = transforms[idx_b].position
166 + transforms[idx_b]
167 .rotation
168 .mul_vec3(rigid_bodies[idx_b].center_of_mass);
169 let r_a = contact_pt - com_a;
170 let r_b = contact_pt - com_b;
171
172 let inv_m_a = rigid_bodies[idx_a].inv_mass();
173 let inv_m_b = rigid_bodies[idx_b].inv_mass();
174 let inv_i_a =
175 rigid_bodies[idx_a].inv_world_inertia_tensor(transforms[idx_a].rotation);
176 let inv_i_b =
177 rigid_bodies[idx_b].inv_world_inertia_tensor(transforms[idx_b].rotation);
178 let dyn_a = rigid_bodies[idx_a].is_dynamic();
179 let dyn_b = rigid_bodies[idx_b].is_dynamic();
180
181 if !dyn_a && !dyn_b {
182 continue;
183 }
184
185 let va = velocities[idx_a].linear + velocities[idx_a].angular.cross(r_a);
187 let vb = velocities[idx_b].linear + velocities[idx_b].angular.cross(r_b);
188 let rel_vel = vb - va;
189 let vel_norm = rel_vel.dot(normal);
190
191 let r_a_x_n = r_a.cross(normal);
193 let r_b_x_n = r_b.cross(normal);
194 let k_n = inv_m_a
195 + inv_m_b
196 + (inv_i_a * r_a_x_n).dot(r_a_x_n)
197 + (inv_i_b * r_b_x_n).dot(r_b_x_n);
198
199 if k_n < 1e-8 {
200 continue;
201 }
202
203 let bias = if penetration < 0.0 {
207 penetration * inv_dt
209 } else if self.split_impulse_enabled {
210 0.0
213 } else {
214 let correction = (penetration - self.slop)
216 .max(0.0)
217 .min(self.max_linear_correction);
218 self.baumgarte * inv_dt * correction
219 };
220
221 let e = if -vel_norm > self.restitution_velocity_threshold {
223 restitution
224 } else {
225 0.0
226 };
227
228 let delta_n = (-(1.0 + e) * vel_norm + bias) / k_n;
229 let new_acc_n = (acc_n + delta_n).max(0.0); let actual_n = new_acc_n - acc_n;
231 manifolds[mid].contacts[cid].normal_impulse = new_acc_n;
232
233 let imp_n = normal * actual_n;
234 if dyn_a {
235 velocities[idx_a].linear -= imp_n * inv_m_a;
236 velocities[idx_a].angular -= inv_i_a * r_a.cross(imp_n);
237 }
238 if dyn_b {
239 velocities[idx_b].linear += imp_n * inv_m_b;
240 velocities[idx_b].angular += inv_i_b * r_b.cross(imp_n);
241 }
242
243 let va2 = velocities[idx_a].linear + velocities[idx_a].angular.cross(r_a);
246 let vb2 = velocities[idx_b].linear + velocities[idx_b].angular.cross(r_b);
247 let rel2 = vb2 - va2;
248 let tang_v = rel2 - normal * rel2.dot(normal);
249 let tang_mag = tang_v.length();
250
251 if tang_mag < 1e-8 && acc_t.length_squared() < 1e-8 {
252 continue;
253 }
254
255 let tangent = if acc_t.length_squared() > 1e-8 {
256 acc_t.normalize()
257 } else if tang_mag > 1e-8 {
258 tang_v / tang_mag
259 } else {
260 if normal.x.abs() > 0.9 {
261 gizmo_math::Vec3::new(0.0, 1.0, 0.0)
262 .cross(normal)
263 .normalize()
264 } else {
265 gizmo_math::Vec3::new(1.0, 0.0, 0.0)
266 .cross(normal)
267 .normalize()
268 }
269 };
270
271 let r_a_x_t = r_a.cross(tangent);
272 let r_b_x_t = r_b.cross(tangent);
273 let k_t = inv_m_a
274 + inv_m_b
275 + (inv_i_a * r_a_x_t).dot(r_a_x_t)
276 + (inv_i_b * r_b_x_t).dot(r_b_x_t);
277
278 if k_t < 1e-8 {
279 continue;
280 }
281
282 let rel_t = rel2.dot(tangent);
283 let delta_t = -rel_t / k_t;
284
285 let static_mu = manifolds[mid].static_friction;
287 let dynamic_mu = friction;
288 let max_static = static_mu * new_acc_n.abs();
289 let max_dynamic = dynamic_mu * new_acc_n.abs();
290
291 let old_t_along = acc_t.dot(tangent);
293 let new_t_along = old_t_along + delta_t;
294
295 let clamped_t = if new_t_along.abs() <= max_static {
296 new_t_along } else {
298 new_t_along.signum() * max_dynamic };
300
301 let actual_t = clamped_t - old_t_along;
302 manifolds[mid].contacts[cid].tangent_impulse = tangent * clamped_t;
303
304 let imp_t = tangent * actual_t;
305 if dyn_a {
306 velocities[idx_a].linear -= imp_t * inv_m_a;
307 velocities[idx_a].angular -= inv_i_a * r_a.cross(imp_t);
308 }
309 if dyn_b {
310 velocities[idx_b].linear += imp_t * inv_m_b;
311 velocities[idx_b].angular += inv_i_b * r_b.cross(imp_t);
312 }
313 }
314 }
315 }
316
317 if self.split_impulse_enabled {
322 let mut acc_pseudo: Vec<Vec<f32>> = manifolds
324 .iter()
325 .map(|m| vec![0.0f32; m.contacts.len()])
326 .collect();
327
328 let pos_iterations = (self.iterations / 2).max(4);
329 for _ in 0..pos_iterations {
330 for mid in 0..manifolds.len() {
331 let entity_a_id = manifolds[mid].entity_a.id();
332 let entity_b_id = manifolds[mid].entity_b.id();
333
334 let idx_a = match entity_index_map.get(&entity_a_id) {
335 Some(&i) => i,
336 None => continue,
337 };
338 let idx_b = match entity_index_map.get(&entity_b_id) {
339 Some(&i) => i,
340 None => continue,
341 };
342
343 let inv_m_a = rigid_bodies[idx_a].inv_mass();
344 let inv_m_b = rigid_bodies[idx_b].inv_mass();
345 let inv_i_a =
346 rigid_bodies[idx_a].inv_world_inertia_tensor(transforms[idx_a].rotation);
347 let inv_i_b =
348 rigid_bodies[idx_b].inv_world_inertia_tensor(transforms[idx_b].rotation);
349 let dyn_a = rigid_bodies[idx_a].is_dynamic();
350 let dyn_b = rigid_bodies[idx_b].is_dynamic();
351
352 if !dyn_a && !dyn_b {
353 continue;
354 }
355
356 let com_a = transforms[idx_a].position
357 + transforms[idx_a]
358 .rotation
359 .mul_vec3(rigid_bodies[idx_a].center_of_mass);
360 let com_b = transforms[idx_b].position
361 + transforms[idx_b]
362 .rotation
363 .mul_vec3(rigid_bodies[idx_b].center_of_mass);
364
365 for cid in 0..manifolds[mid].contacts.len() {
366 let contact_pt = manifolds[mid].contacts[cid].point;
367 let normal = manifolds[mid].contacts[cid].normal;
368 let penetration = manifolds[mid].contacts[cid].penetration;
369
370 let correction = (penetration - self.slop)
371 .max(0.0)
372 .min(self.max_linear_correction);
373 if correction < 1e-6 {
374 continue;
375 }
376
377 let r_a = contact_pt - com_a;
378 let r_b = contact_pt - com_b;
379
380 let r_a_x_n = r_a.cross(normal);
381 let r_b_x_n = r_b.cross(normal);
382 let k_n = inv_m_a
383 + inv_m_b
384 + (inv_i_a * r_a_x_n).dot(r_a_x_n)
385 + (inv_i_b * r_b_x_n).dot(r_b_x_n);
386 if k_n < 1e-8 {
387 continue;
388 }
389
390 let pv_a = pseudo_vel[idx_a].0 + pseudo_vel[idx_a].1.cross(r_a);
392 let pv_b = pseudo_vel[idx_b].0 + pseudo_vel[idx_b].1.cross(r_b);
393 let pv_rel = pv_b.dot(normal) - pv_a.dot(normal);
394
395 let bias = self.split_impulse_erp * inv_dt * correction;
396 let delta_p = (-pv_rel + bias) / k_n;
400
401 let old_acc = acc_pseudo[mid][cid];
403 let new_acc = (old_acc + delta_p).max(0.0);
404 let actual_delta = new_acc - old_acc;
405 acc_pseudo[mid][cid] = new_acc;
406
407 let imp_p = normal * actual_delta;
408 if dyn_a {
409 pseudo_vel[idx_a].0 -= imp_p * inv_m_a;
410 pseudo_vel[idx_a].1 -= inv_i_a * r_a.cross(imp_p);
411 }
412 if dyn_b {
413 pseudo_vel[idx_b].0 += imp_p * inv_m_b;
414 pseudo_vel[idx_b].1 += inv_i_b * r_b.cross(imp_p);
415 }
416 }
417 }
418 }
419
420 for i in 0..velocities.len() {
422 velocities[i].linear += pseudo_vel[i].0;
423 velocities[i].angular += pseudo_vel[i].1;
424 }
425 }
426 }
427
428 #[allow(clippy::too_many_arguments)]
433 pub fn solve_contact_constraint(
434 &self,
435 rb_a: &mut RigidBody,
436 transform_a: &Transform,
437 vel_a: &mut Velocity,
438 rb_b: &mut RigidBody,
439 transform_b: &Transform,
440 vel_b: &mut Velocity,
441 contact_point: Vec3,
442 normal: Vec3,
443 penetration: f32,
444 static_friction: f32,
445 dynamic_friction: f32,
446 restitution: f32,
447 dt: f32,
448 ) {
449 if !rb_a.is_dynamic() && !rb_b.is_dynamic() {
450 return;
451 }
452
453 let com_a = transform_a.position + transform_a.rotation.mul_vec3(rb_a.center_of_mass);
454 let com_b = transform_b.position + transform_b.rotation.mul_vec3(rb_b.center_of_mass);
455
456 let r_a = contact_point - com_a;
457 let r_b = contact_point - com_b;
458
459 let va = vel_a.linear + vel_a.angular.cross(r_a);
460 let vb = vel_b.linear + vel_b.angular.cross(r_b);
461 let rel_vel = vb - va;
462 let vel_norm = rel_vel.dot(normal);
463
464 if vel_norm > 0.0 {
465 return;
466 } let inv_m_a = rb_a.inv_mass();
469 let inv_m_b = rb_b.inv_mass();
470 let inv_i_a = rb_a.inv_world_inertia_tensor(transform_a.rotation);
471 let inv_i_b = rb_b.inv_world_inertia_tensor(transform_b.rotation);
472
473 let r_a_x_n = r_a.cross(normal);
474 let r_b_x_n = r_b.cross(normal);
475 let k =
476 inv_m_a + inv_m_b + (inv_i_a * r_a_x_n).dot(r_a_x_n) + (inv_i_b * r_b_x_n).dot(r_b_x_n);
477
478 if k < 1e-8 {
479 return;
480 }
481
482 let inv_dt = if dt > 0.0 { 1.0 / dt } else { 0.0 };
483 let bias = self.baumgarte * inv_dt * (penetration - self.slop).max(0.0);
484 let e = if -vel_norm > self.restitution_velocity_threshold {
485 restitution
486 } else {
487 0.0
488 };
489 let j = ((-(1.0 + e) * vel_norm + bias) / k).max(0.0);
490
491 let impulse = normal * j;
492
493 if rb_a.is_dynamic() {
494 vel_a.linear -= impulse * inv_m_a;
495 vel_a.angular -= inv_i_a * r_a.cross(impulse);
496 }
497 if rb_b.is_dynamic() {
498 vel_b.linear += impulse * inv_m_b;
499 vel_b.angular += inv_i_b * r_b.cross(impulse);
500 }
501
502 self.apply_friction_standalone(
504 rb_a,
505 vel_a,
506 rb_b,
507 vel_b,
508 r_a,
509 r_b,
510 normal,
511 static_friction,
512 dynamic_friction,
513 j,
514 &inv_i_a,
515 &inv_i_b,
516 );
517 }
518
519 #[allow(clippy::too_many_arguments)]
520 fn apply_friction_standalone(
521 &self,
522 rb_a: &RigidBody,
523 vel_a: &mut Velocity,
524 rb_b: &RigidBody,
525 vel_b: &mut Velocity,
526 r_a: Vec3,
527 r_b: Vec3,
528 normal: Vec3,
529 static_friction: f32,
530 dynamic_friction: f32,
531 normal_impulse: f32,
532 inv_i_a: &Mat3,
533 inv_i_b: &Mat3,
534 ) {
535 let va = vel_a.linear + vel_a.angular.cross(r_a);
536 let vb = vel_b.linear + vel_b.angular.cross(r_b);
537 let rel = vb - va;
538 let tang_v = rel - normal * rel.dot(normal);
539 let tang_mag = tang_v.length();
540 if tang_mag < 1e-8 {
541 return;
542 }
543 let tangent = tang_v / tang_mag;
544
545 let inv_m_a = rb_a.inv_mass();
546 let inv_m_b = rb_b.inv_mass();
547
548 let r_a_x_t = r_a.cross(tangent);
549 let r_b_x_t = r_b.cross(tangent);
550 let k = inv_m_a
551 + inv_m_b
552 + (*inv_i_a * r_a_x_t).dot(r_a_x_t)
553 + (*inv_i_b * r_b_x_t).dot(r_b_x_t);
554 if k < 1e-8 {
555 return;
556 }
557
558 let max_static = static_friction * normal_impulse.abs();
559 let max_dynamic = dynamic_friction * normal_impulse.abs();
560
561 let delta_t = -tang_mag / k;
562
563 let jt = if delta_t.abs() <= max_static {
564 delta_t
565 } else {
566 delta_t.signum() * max_dynamic
567 };
568
569 let ft = tangent * jt;
570 if rb_a.is_dynamic() {
571 vel_a.linear -= ft * inv_m_a;
572 vel_a.angular -= *inv_i_a * r_a.cross(ft);
573 }
574 if rb_b.is_dynamic() {
575 vel_b.linear += ft * inv_m_b;
576 vel_b.angular += *inv_i_b * r_b.cross(ft);
577 }
578 }
579}
580
581#[cfg(test)]
586mod tests {
587 use super::*;
588 use gizmo_math::Vec3;
589
590 #[test]
591 fn test_solver_creation() {
592 let solver = ConstraintSolver::new(20);
593 assert_eq!(solver.iterations, 20);
594 }
595
596 #[test]
597 fn test_collision_response() {
598 let mut rb_a = RigidBody::default();
599 let mut rb_b = RigidBody::default();
600 rb_a.wake_up();
601 rb_b.wake_up();
602 let transform_a = Transform::new(Vec3::new(0.0, 0.0, 0.0));
603 let transform_b = Transform::new(Vec3::new(0.0, 2.0, 0.0));
604 let mut vel_a = Velocity::new(Vec3::new(0.0, 1.0, 0.0));
605 let mut vel_b = Velocity::new(Vec3::new(0.0, -1.0, 0.0));
606
607 let solver = ConstraintSolver::default();
608 solver.solve_contact_constraint(
609 &mut rb_a,
610 &transform_a,
611 &mut vel_a,
612 &mut rb_b,
613 &transform_b,
614 &mut vel_b,
615 Vec3::new(0.0, 1.0, 0.0),
616 Vec3::new(0.0, 1.0, 0.0),
617 0.1,
618 0.6,
619 0.5,
620 0.5,
621 0.016,
622 );
623
624 assert!(vel_a.linear.y < 1.0);
625 assert!(vel_b.linear.y > -1.0);
626 }
627
628 #[test]
629 fn test_normal_impulse_non_negative() {
630 let mut rb_a = RigidBody::default();
631 let mut rb_b = RigidBody::default();
632 rb_a.wake_up();
633 rb_b.wake_up();
634 let transform_a = Transform::new(Vec3::ZERO);
635 let transform_b = Transform::new(Vec3::new(0.0, 1.0, 0.0));
636 let mut vel_a = Velocity::new(Vec3::new(0.0, 5.0, 0.0));
637 let mut vel_b = Velocity::new(Vec3::new(0.0, -5.0, 0.0));
638
639 let before_a = vel_a.linear.y;
640 let before_b = vel_b.linear.y;
641
642 let solver = ConstraintSolver::default();
643 solver.solve_contact_constraint(
644 &mut rb_a,
645 &transform_a,
646 &mut vel_a,
647 &mut rb_b,
648 &transform_b,
649 &mut vel_b,
650 Vec3::new(0.0, 0.5, 0.0),
651 Vec3::new(0.0, 1.0, 0.0),
652 0.05,
653 0.4,
654 0.3,
655 0.0,
656 0.016,
657 );
658
659 assert!(vel_a.linear.y < before_a);
660 assert!(vel_b.linear.y > before_b);
661 }
662
663 #[test]
664 fn test_resting_contact_no_bounce() {
665 let mut rb_a = RigidBody::default();
667 let mut rb_b = RigidBody::default();
668 rb_b.body_type = crate::components::rigid_body::BodyType::Static;
669 rb_a.wake_up();
670 let transform_a = Transform::new(Vec3::new(0.0, 1.05, 0.0));
671 let transform_b = Transform::new(Vec3::ZERO);
672 let mut vel_a = Velocity::new(Vec3::new(0.0, -0.1, 0.0)); let mut vel_b = Velocity::default();
674
675 let solver = ConstraintSolver::default();
676 solver.solve_contact_constraint(
681 &mut rb_a,
682 &transform_a,
683 &mut vel_a,
684 &mut rb_b,
685 &transform_b,
686 &mut vel_b,
687 Vec3::new(0.0, 0.525, 0.0), Vec3::new(0.0, -1.0, 0.0), 0.05,
690 0.6,
691 0.5,
692 0.0,
693 0.016, );
695 assert!(
696 vel_a.linear.y >= -0.01,
697 "A should not bounce or sink significantly: {}",
698 vel_a.linear.y
699 );
700 assert!(
701 vel_b.linear.y <= 0.01,
702 "B (static) should remain still: {}",
703 vel_b.linear.y
704 );
705 }
706}