diff --git a/CHANGELOG.md b/CHANGELOG.md
index a7c1efed8..f9ec4d1a8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -14,12 +14,21 @@
 
 ### Modified
 
+- The contact constraints regularization parameters have been changed from `erp/damping_ratio` to
+  `natural_frequency/damping_ratio`. This helps define them in a timestep-length independent way. The new variables
+  are named `IntegrationParameters::contact_natural_frequency` and `IntegrationParameters::contact_damping_ratio`.
+- The `IntegrationParameters::normalized_max_penetration_correction` has been replaced
+  by `::normalized_max_corrective_velocity`
+  to make the parameter more timestep-length independent. It is now set to a non-infinite value to eliminate aggressive
+  "popping effects".
 - The `Multibody::forward_kinematics` method will no longer automatically update the poses of the `RigidBody` associated
   to each joint. Instead `Multibody::update_rigid_bodies` has to be called explicitly.
 - The `Multibody::forward_kinematics` method will automatically adjust the multibody’s degrees of freedom if the root
   rigid-body changed type (between dynamic and non-dynamic). It can also optionally apply the root’s rigid-body pose
   instead of the root link’s pose (useful for example if you modified the root rigid-body pose externally and wanted
   to propagate it to the multibody).
+- Remove an internal special-case for contact constraints on fast contacts. The doesn’t seem necessary with the substep
+  solver.
 
 ## v0.19.0 (05 May 2024)
 
diff --git a/benchmarks3d/all_benchmarks3.rs b/benchmarks3d/all_benchmarks3.rs
index 2f723d1e4..c1e0da28f 100644
--- a/benchmarks3d/all_benchmarks3.rs
+++ b/benchmarks3d/all_benchmarks3.rs
@@ -20,6 +20,7 @@ mod joint_fixed3;
 mod joint_prismatic3;
 mod joint_revolute3;
 mod keva3;
+mod many_pyramids3;
 mod many_sleep3;
 mod many_static3;
 mod pyramid3;
@@ -66,6 +67,7 @@ pub fn main() {
         ("ImpulseJoint fixed", joint_fixed3::init_world),
         ("ImpulseJoint revolute", joint_revolute3::init_world),
         ("ImpulseJoint prismatic", joint_prismatic3::init_world),
+        ("Many pyramids", many_pyramids3::init_world),
         ("Keva tower", keva3::init_world),
     ];
 
diff --git a/benchmarks3d/many_pyramids3.rs b/benchmarks3d/many_pyramids3.rs
new file mode 100644
index 000000000..1562aec91
--- /dev/null
+++ b/benchmarks3d/many_pyramids3.rs
@@ -0,0 +1,80 @@
+use rapier3d::prelude::*;
+use rapier_testbed3d::Testbed;
+
+fn create_pyramid(
+    bodies: &mut RigidBodySet,
+    colliders: &mut ColliderSet,
+    offset: Vector<f32>,
+    stack_height: usize,
+    rad: f32,
+) {
+    let shift = rad * 2.0;
+
+    for i in 0usize..stack_height {
+        for j in i..stack_height {
+            let fj = j as f32;
+            let fi = i as f32;
+            let x = (fi * shift / 2.0) + (fj - fi) * shift;
+            let y = fi * shift;
+
+            // Build the rigid body.
+            let rigid_body = RigidBodyBuilder::dynamic().translation(vector![x, y, 0.0] + offset);
+            let handle = bodies.insert(rigid_body);
+            let collider = ColliderBuilder::cuboid(rad, rad, rad);
+            colliders.insert_with_parent(collider, handle, bodies);
+        }
+    }
+}
+
+pub fn init_world(testbed: &mut Testbed) {
+    /*
+     * World
+     */
+    let mut bodies = RigidBodySet::new();
+    let mut colliders = ColliderSet::new();
+    let impulse_joints = ImpulseJointSet::new();
+    let multibody_joints = MultibodyJointSet::new();
+
+    let rad = 0.5;
+    let pyramid_count = 40;
+    let spacing = 4.0;
+
+    /*
+     * Ground
+     */
+    let ground_size = 50.0;
+    let ground_height = 0.1;
+
+    let rigid_body = RigidBodyBuilder::fixed().translation(vector![0.0, -ground_height, 0.0]);
+    let ground_handle = bodies.insert(rigid_body);
+    let collider = ColliderBuilder::cuboid(
+        ground_size,
+        ground_height,
+        pyramid_count as f32 * spacing / 2.0 + ground_size,
+    );
+    colliders.insert_with_parent(collider, ground_handle, &mut bodies);
+
+    /*
+     * Create the cubes
+     */
+    for pyramid_index in 0..pyramid_count {
+        let bottomy = rad;
+        create_pyramid(
+            &mut bodies,
+            &mut colliders,
+            vector![
+                0.0,
+                bottomy,
+                (pyramid_index as f32 - pyramid_count as f32 / 2.0) * spacing
+            ],
+            20,
+            rad,
+        );
+    }
+
+    /*
+     * Set up the testbed.
+     */
+    testbed.set_world(bodies, colliders, impulse_joints, multibody_joints);
+    testbed.look_at(point![100.0, 100.0, 100.0], Point::origin());
+}
diff --git a/src/dynamics/integration_parameters.rs b/src/dynamics/integration_parameters.rs
index c6621282b..2de58aecb 100644
--- a/src/dynamics/integration_parameters.rs
+++ b/src/dynamics/integration_parameters.rs
@@ -1,4 +1,5 @@
 use crate::math::Real;
+use na::RealField;
 use std::num::NonZeroUsize;
 
 // TODO: enabling the block solver in 3d introduces a lot of jitters in
@@ -9,9 +10,9 @@ pub(crate) static BLOCK_SOLVER_ENABLED: bool = cfg!(feature = "dim2");
 #[derive(Copy, Clone, Debug)]
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
 pub struct IntegrationParameters {
-    /// The timestep length (default: `1.0 / 60.0`)
+    /// The timestep length (default: `1.0 / 60.0`).
     pub dt: Real,
-    /// Minimum timestep size when using CCD with multiple substeps (default `1.0 / 60.0 / 100.0`)
+    /// Minimum timestep size when using CCD with multiple substeps (default: `1.0 / 60.0 / 100.0`).
     ///
     /// When CCD with multiple substeps is enabled, the timestep is subdivided
     /// into smaller pieces. This timestep subdivision won't generate timestep
@@ -23,29 +24,40 @@ pub struct IntegrationParameters {
     /// to numerical instabilities.
     pub min_ccd_dt: Real,
 
-    /// 0-1: multiplier for how much of the constraint violation (e.g. contact penetration)
-    /// will be compensated for during the velocity solve.
-    /// (default `0.1`).
-    pub erp: Real,
-    /// 0-1: the damping ratio used by the springs for Baumgarte constraints stabilization.
-    /// Lower values make the constraints more compliant (more "springy", allowing more visible penetrations
-    /// before stabilization).
-    /// (default `20.0`).
-    pub damping_ratio: Real,
+    /// > 0: the damping ratio used by the springs for contact constraint stabilization.
+    ///
+    /// Larger values make the constraints more compliant (allowing more visible
+    /// penetrations before stabilization).
+    /// (default `5.0`).
+    pub contact_damping_ratio: Real,
 
-    /// 0-1: multiplier for how much of the joint violation
-    /// will be compensated for during the velocity solve.
-    /// (default `1.0`).
-    pub joint_erp: Real,
+    /// > 0: the natural frequency used by the springs for contact constraint regularization.
+    ///
+    /// Increasing this value will make it so that penetrations get fixed more quickly at the
+    /// expense of potential jitter effects due to overshooting. In order to make the simulation
+    /// look stiffer, it is recommended to increase the [`Self::contact_damping_ratio`] instead of this
+    /// value.
+    /// (default: `30.0`).
+    pub contact_natural_frequency: Real,
+
+    /// > 0: the natural frequency used by the springs for joint constraint regularization.
+    ///
+    /// Increasing this value will make it so that penetrations get fixed more quickly.
+    /// (default: `1.0e6`).
+    pub joint_natural_frequency: Real,
 
     /// The fraction of critical damping applied to the joint for constraints regularization.
+    ///
+    /// Larger values make the constraints more compliant (allowing more joint
+    /// drift before stabilization).
     /// (default `1.0`).
     pub joint_damping_ratio: Real,
 
     /// The coefficient in `[0, 1]` applied to warmstart impulses, i.e., impulses that are used as the
     /// initial solution (instead of 0) at the next simulation step.
     ///
-    /// This should generally be set to 1. Can be set to 0 if using a large [`Self::erp`] value.
+    /// This should generally be set to 1.
+    ///
     /// (default `1.0`).
     pub warmstart_coefficient: Real,
 
@@ -53,13 +65,13 @@ pub struct IntegrationParameters {
     ///
     /// This value is used internally to estimate some length-based tolerance. In particular, the
     /// values [`IntegrationParameters::allowed_linear_error`],
-    /// [`IntegrationParameters::max_penetration_correction`],
+    /// [`IntegrationParameters::max_corrective_velocity`],
     /// [`IntegrationParameters::prediction_distance`], [`RigidBodyActivation::linear_threshold`]
     /// are scaled by this value implicitly.
     ///
     /// This value can be understood as the number of units-per-meter in your physical world compared
     /// to a human-sized world in meter. For example, in a 2d game, if your typical object size is 100
-    /// pixels, set the `[`Self::length_unit`]` parameter to 100.0. The physics engine will interpret
+    /// pixels, set the [`Self::length_unit`] parameter to 100.0. The physics engine will interpret
     /// it as if 100 pixels is equivalent to 1 meter in its various internal threshold.
     /// (default `1.0`).
     pub length_unit: Real,
@@ -68,17 +80,17 @@ pub struct IntegrationParameters {
     ///
     /// This value is implicitly scaled by [`IntegrationParameters::length_unit`].
     pub normalized_allowed_linear_error: Real,
-    /// Maximum amount of penetration the solver will attempt to resolve in one timestep.
+    /// Maximum amount of penetration the solver will attempt to resolve in one timestep (default: `10.0`).
     ///
     /// This value is implicitly scaled by [`IntegrationParameters::length_unit`].
-    pub normalized_max_penetration_correction: Real,
+    pub normalized_max_corrective_velocity: Real,
     /// The maximal distance separating two objects that will generate predictive contacts (default: `0.002m`).
     ///
     /// This value is implicitly scaled by [`IntegrationParameters::length_unit`].
     pub normalized_prediction_distance: Real,
     /// The number of solver iterations run by the constraints solver for calculating forces (default: `4`).
     pub num_solver_iterations: NonZeroUsize,
-    /// Number of addition friction resolution iteration run during the last solver sub-step (default: `4`).
+    /// Number of addition friction resolution iteration run during the last solver sub-step (default: `0`).
     pub num_additional_friction_iterations: usize,
     /// Number of internal Project Gauss Seidel (PGS) iterations run at each solver iteration (default: `1`).
     pub num_internal_pgs_iterations: usize,
@@ -123,20 +135,53 @@ impl IntegrationParameters {
         }
     }
 
-    /// The ERP coefficient, multiplied by the inverse timestep length.
-    pub fn erp_inv_dt(&self) -> Real {
-        self.erp * self.inv_dt()
+    /// The contact’s spring angular frequency for constraints regularization.
+    pub fn contact_angular_frequency(&self) -> Real {
+        self.contact_natural_frequency * Real::two_pi()
+    }
+
+    /// The [`Self::contact_erp`] coefficient, multiplied by the inverse timestep length.
+    pub fn contact_erp_inv_dt(&self) -> Real {
+        let ang_freq = self.contact_angular_frequency();
+        ang_freq / (self.dt * ang_freq + 2.0 * self.contact_damping_ratio)
+    }
+
+    /// The effective Error Reduction Parameter applied for calculating regularization forces
+    /// on contacts.
+    ///
+    /// This parameter is computed automatically from [`Self::contact_natural_frequency`],
+    /// [`Self::contact_damping_ratio`] and the substep length.
+    pub fn contact_erp(&self) -> Real {
+        self.dt * self.contact_erp_inv_dt()
     }
 
-    /// The joint ERP coefficient, multiplied by the inverse timestep length.
+    /// The joint’s spring angular frequency for constraint regularization.
+    pub fn joint_angular_frequency(&self) -> Real {
+        self.joint_natural_frequency * Real::two_pi()
+    }
+
+    /// The [`Self::joint_erp`] coefficient, multiplied by the inverse timestep length.
     pub fn joint_erp_inv_dt(&self) -> Real {
-        self.joint_erp * self.inv_dt()
+        let ang_freq = self.joint_angular_frequency();
+        ang_freq / (self.dt * ang_freq + 2.0 * self.joint_damping_ratio)
+    }
+
+    /// The effective Error Reduction Parameter applied for calculating regularization forces
+    /// on joints.
+    ///
+    /// This parameter is computed automatically from [`Self::joint_natural_frequency`],
+    /// [`Self::joint_damping_ratio`] and the substep length.
+    pub fn joint_erp(&self) -> Real {
+        self.dt * self.joint_erp_inv_dt()
     }
 
-    /// The CFM factor to be used in the constraints resolution.
-    pub fn cfm_factor(&self) -> Real {
+    /// The CFM factor to be used in the constraint resolution.
+    ///
+    /// This parameter is computed automatically from [`Self::contact_natural_frequency`],
+    /// [`Self::contact_damping_ratio`] and the substep length.
+    pub fn contact_cfm_factor(&self) -> Real {
         // Compute CFM assuming a critically damped spring multiplied by the damping ratio.
-        let inv_erp_minus_one = 1.0 / self.erp - 1.0;
+        let inv_erp_minus_one = 1.0 / self.contact_erp() - 1.0;
 
         // let stiffness = 4.0 * damping_ratio * damping_ratio * projected_mass
         //     / (dt * dt * inv_erp_minus_one * inv_erp_minus_one);
@@ -145,7 +190,10 @@ impl IntegrationParameters {
         // let cfm = 1.0 / (dt * dt * stiffness + dt * damping);
         // NOTE: This simplifies to cfm = cfm_coeff / projected_mass:
         let cfm_coeff = inv_erp_minus_one * inv_erp_minus_one
-            / ((1.0 + inv_erp_minus_one) * 4.0 * self.damping_ratio * self.damping_ratio);
+            / ((1.0 + inv_erp_minus_one)
+                * 4.0
+                * self.contact_damping_ratio
+                * self.contact_damping_ratio);
 
         // Furthermore, we use this coefficient inside of the impulse resolution.
         // Surprisingly, several simplifications happen there.
@@ -166,11 +214,14 @@ impl IntegrationParameters {
         1.0 / (1.0 + cfm_coeff)
     }
 
-    /// The CFM (constraints force mixing) coefficient applied to all joints for constraints regularization
+    /// The CFM (constraints force mixing) coefficient applied to all joints for constraints regularization.
+    ///
+    /// This parameter is computed automatically from [`Self::joint_natural_frequency`],
+    /// [`Self::joint_damping_ratio`] and the substep length.
     pub fn joint_cfm_coeff(&self) -> Real {
         // Compute CFM assuming a critically damped spring multiplied by the damping ratio.
         // The logic is similar to `Self::cfm_factor`.
-        let inv_erp_minus_one = 1.0 / self.joint_erp - 1.0;
+        let inv_erp_minus_one = 1.0 / self.joint_erp() - 1.0;
         inv_erp_minus_one * inv_erp_minus_one
             / ((1.0 + inv_erp_minus_one)
                 * 4.0
@@ -186,11 +237,11 @@ impl IntegrationParameters {
 
     /// Maximum amount of penetration the solver will attempt to resolve in one timestep.
     ///
-    /// This is equal to [`Self::normalized_max_penetration_correction`] multiplied by
+    /// This is equal to [`Self::normalized_max_corrective_velocity`] multiplied by
     /// [`Self::length_unit`].
-    pub fn max_penetration_correction(&self) -> Real {
-        if self.normalized_max_penetration_correction != Real::MAX {
-            self.normalized_max_penetration_correction * self.length_unit
+    pub fn max_corrective_velocity(&self) -> Real {
+        if self.normalized_max_corrective_velocity != Real::MAX {
+            self.normalized_max_corrective_velocity * self.length_unit
         } else {
             Real::MAX
         }
@@ -210,9 +261,9 @@ impl IntegrationParameters {
         Self {
             dt: 1.0 / 60.0,
             min_ccd_dt: 1.0 / 60.0 / 100.0,
-            erp: 0.1,
-            damping_ratio: 20.0,
-            joint_erp: 1.0,
+            contact_natural_frequency: 30.0,
+            contact_damping_ratio: 5.0,
+            joint_natural_frequency: 1.0e6,
             joint_damping_ratio: 1.0,
             warmstart_coefficient: 1.0,
             num_internal_pgs_iterations: 1,
@@ -226,7 +277,7 @@ impl IntegrationParameters {
             // tons of islands, reducing SIMD parallelism opportunities.
             min_island_size: 128,
             normalized_allowed_linear_error: 0.001,
-            normalized_max_penetration_correction: Real::MAX,
+            normalized_max_corrective_velocity: 10.0,
             normalized_prediction_distance: 0.002,
             max_ccd_substeps: 1,
             length_unit: 1.0,
@@ -240,8 +291,7 @@ impl IntegrationParameters {
     /// warmstarting proves to be undesirable for your use-case.
     pub fn tgs_soft_without_warmstart() -> Self {
         Self {
-            erp: 0.6,
-            damping_ratio: 1.0,
+            contact_damping_ratio: 0.25,
             warmstart_coefficient: 0.0,
             num_additional_friction_iterations: 4,
             ..Self::tgs_soft()
@@ -253,12 +303,9 @@ impl IntegrationParameters {
     /// This exists mainly for testing and comparison purpose.
     pub fn pgs_legacy() -> Self {
         Self {
-            erp: 0.8,
-            damping_ratio: 0.25,
-            warmstart_coefficient: 0.0,
-            num_additional_friction_iterations: 4,
             num_solver_iterations: NonZeroUsize::new(1).unwrap(),
-            ..Self::tgs_soft()
+            num_internal_pgs_iterations: 4,
+            ..Self::tgs_soft_without_warmstart()
         }
     }
 }
diff --git a/src/dynamics/solver/contact_constraint/generic_one_body_constraint.rs b/src/dynamics/solver/contact_constraint/generic_one_body_constraint.rs
index 301ce7148..e25499569 100644
--- a/src/dynamics/solver/contact_constraint/generic_one_body_constraint.rs
+++ b/src/dynamics/solver/contact_constraint/generic_one_body_constraint.rs
@@ -231,13 +231,8 @@ impl GenericOneBodyConstraintBuilder {
             .unwrap()
             .local_to_world;
 
-        self.inner.update_with_positions(
-            params,
-            solved_dt,
-            pos2,
-            self.ccd_thickness,
-            &mut constraint.inner,
-        );
+        self.inner
+            .update_with_positions(params, solved_dt, pos2, &mut constraint.inner);
     }
 }
 
diff --git a/src/dynamics/solver/contact_constraint/generic_two_body_constraint.rs b/src/dynamics/solver/contact_constraint/generic_two_body_constraint.rs
index be839f01e..dd63d1f95 100644
--- a/src/dynamics/solver/contact_constraint/generic_two_body_constraint.rs
+++ b/src/dynamics/solver/contact_constraint/generic_two_body_constraint.rs
@@ -342,14 +342,8 @@ impl GenericTwoBodyConstraintBuilder {
             .map(|m| &multibodies[m.multibody].link(m.id).unwrap().local_to_world)
             .unwrap_or_else(|| &bodies[constraint.inner.solver_vel2].position);
 
-        self.inner.update_with_positions(
-            params,
-            solved_dt,
-            pos1,
-            pos2,
-            self.ccd_thickness,
-            &mut constraint.inner,
-        );
+        self.inner
+            .update_with_positions(params, solved_dt, pos1, pos2, &mut constraint.inner);
     }
 }
 
diff --git a/src/dynamics/solver/contact_constraint/one_body_constraint.rs b/src/dynamics/solver/contact_constraint/one_body_constraint.rs
index 1243d11a4..5e6e86b3c 100644
--- a/src/dynamics/solver/contact_constraint/one_body_constraint.rs
+++ b/src/dynamics/solver/contact_constraint/one_body_constraint.rs
@@ -260,13 +260,7 @@ impl OneBodyConstraintBuilder {
         constraint: &mut OneBodyConstraint,
     ) {
         let rb2 = &bodies[constraint.solver_vel2];
-        self.update_with_positions(
-            params,
-            solved_dt,
-            &rb2.position,
-            rb2.ccd_thickness,
-            constraint,
-        )
+        self.update_with_positions(params, solved_dt, &rb2.position, constraint)
     }
 
     // TODO: this code is SOOOO similar to TwoBodyConstraint::update.
@@ -276,12 +270,11 @@ impl OneBodyConstraintBuilder {
         params: &IntegrationParameters,
         solved_dt: Real,
         rb2_pos: &Isometry<Real>,
-        ccd_thickness: Real,
         constraint: &mut OneBodyConstraint,
     ) {
-        let cfm_factor = params.cfm_factor();
+        let cfm_factor = params.contact_cfm_factor();
         let inv_dt = params.inv_dt();
-        let erp_inv_dt = params.erp_inv_dt();
+        let erp_inv_dt = params.contact_erp_inv_dt();
 
         let all_infos = &self.infos[..constraint.num_contacts as usize];
         let all_elements = &mut constraint.elements[..constraint.num_contacts as usize];
@@ -290,7 +283,6 @@ impl OneBodyConstraintBuilder {
         let new_pos1 = self
             .vels1
             .integrate(solved_dt, &rb1.position, &rb1.local_com);
-        let mut is_fast_contact = false;
 
         #[cfg(feature = "dim2")]
         let tangents1 = constraint.dir1.orthonormal_basis();
@@ -309,11 +301,9 @@ impl OneBodyConstraintBuilder {
             // Normal part.
             {
                 let rhs_wo_bias = info.normal_rhs_wo_bias + dist.max(0.0) * inv_dt;
-                let rhs_bias = erp_inv_dt
-                    * (dist + params.allowed_linear_error())
-                        .clamp(-params.max_penetration_correction(), 0.0);
+                let rhs_bias = (erp_inv_dt * (dist + params.allowed_linear_error()))
+                    .clamp(-params.max_corrective_velocity(), 0.0);
                 let new_rhs = rhs_wo_bias + rhs_bias;
-                is_fast_contact = is_fast_contact || (-new_rhs * params.dt > ccd_thickness * 0.5);
 
                 element.normal_part.rhs_wo_bias = rhs_wo_bias;
                 element.normal_part.rhs = new_rhs;
@@ -333,7 +323,7 @@ impl OneBodyConstraintBuilder {
             }
         }
 
-        constraint.cfm_factor = if is_fast_contact { 1.0 } else { cfm_factor };
+        constraint.cfm_factor = cfm_factor;
     }
 }
 
diff --git a/src/dynamics/solver/contact_constraint/one_body_constraint_simd.rs b/src/dynamics/solver/contact_constraint/one_body_constraint_simd.rs
index 5710bc40b..077dd750c 100644
--- a/src/dynamics/solver/contact_constraint/one_body_constraint_simd.rs
+++ b/src/dynamics/solver/contact_constraint/one_body_constraint_simd.rs
@@ -15,7 +15,6 @@ use crate::math::{
 use crate::utils::SimdBasis;
 use crate::utils::{self, SimdAngularInertia, SimdCross, SimdDot};
 use num::Zero;
-use parry::math::SimdBool;
 use parry::utils::SdpMatrix2;
 use simba::simd::{SimdPartialOrd, SimdValue};
 
@@ -262,16 +261,14 @@ impl SimdOneBodyConstraintBuilder {
         _multibodies: &MultibodyJointSet,
         constraint: &mut OneBodyConstraintSimd,
     ) {
-        let cfm_factor = SimdReal::splat(params.cfm_factor());
-        let dt = SimdReal::splat(params.dt);
+        let cfm_factor = SimdReal::splat(params.contact_cfm_factor());
         let inv_dt = SimdReal::splat(params.inv_dt());
         let allowed_lin_err = SimdReal::splat(params.allowed_linear_error());
-        let erp_inv_dt = SimdReal::splat(params.erp_inv_dt());
-        let max_penetration_correction = SimdReal::splat(params.max_penetration_correction());
+        let erp_inv_dt = SimdReal::splat(params.contact_erp_inv_dt());
+        let max_corrective_velocity = SimdReal::splat(params.max_corrective_velocity());
         let warmstart_coeff = SimdReal::splat(params.warmstart_coefficient);
 
         let rb2 = gather![|ii| &bodies[constraint.solver_vel2[ii]]];
-        let ccd_thickness = SimdReal::from(gather![|ii| rb2[ii].ccd_thickness]);
         let poss2 = Isometry::from(gather![|ii| rb2[ii].position]);
 
         let all_infos = &self.infos[..constraint.num_contacts as usize];
@@ -292,7 +289,6 @@ impl SimdOneBodyConstraintBuilder {
             constraint.dir1.cross(&constraint.tangent1),
         ];
 
-        let mut is_fast_contact = SimdBool::splat(false);
         let solved_dt = SimdReal::splat(solved_dt);
 
         for (info, element) in all_infos.iter().zip(all_elements.iter_mut()) {
@@ -305,12 +301,9 @@ impl SimdOneBodyConstraintBuilder {
             {
                 let rhs_wo_bias =
                     info.normal_rhs_wo_bias + dist.simd_max(SimdReal::zero()) * inv_dt;
-                let rhs_bias = (dist + allowed_lin_err)
-                    .simd_clamp(-max_penetration_correction, SimdReal::zero())
-                    * erp_inv_dt;
+                let rhs_bias = ((dist + allowed_lin_err) * erp_inv_dt)
+                    .simd_clamp(-max_corrective_velocity, SimdReal::zero());
                 let new_rhs = rhs_wo_bias + rhs_bias;
-                is_fast_contact =
-                    is_fast_contact | (-new_rhs * dt).simd_gt(ccd_thickness * SimdReal::splat(0.5));
 
                 element.normal_part.rhs_wo_bias = rhs_wo_bias;
                 element.normal_part.rhs = new_rhs;
@@ -330,7 +323,7 @@ impl SimdOneBodyConstraintBuilder {
             }
         }
 
-        constraint.cfm_factor = SimdReal::splat(1.0).select(is_fast_contact, cfm_factor);
+        constraint.cfm_factor = cfm_factor;
     }
 }
 
diff --git a/src/dynamics/solver/contact_constraint/two_body_constraint.rs b/src/dynamics/solver/contact_constraint/two_body_constraint.rs
index 1511d4f7a..d2a2c4940 100644
--- a/src/dynamics/solver/contact_constraint/two_body_constraint.rs
+++ b/src/dynamics/solver/contact_constraint/two_body_constraint.rs
@@ -361,15 +361,7 @@ impl TwoBodyConstraintBuilder {
     ) {
         let rb1 = &bodies[constraint.solver_vel1];
         let rb2 = &bodies[constraint.solver_vel2];
-        let ccd_thickness = rb1.ccd_thickness + rb2.ccd_thickness;
-        self.update_with_positions(
-            params,
-            solved_dt,
-            &rb1.position,
-            &rb2.position,
-            ccd_thickness,
-            constraint,
-        )
+        self.update_with_positions(params, solved_dt, &rb1.position, &rb2.position, constraint)
     }
 
     // Used by both generic and non-generic builders..
@@ -379,18 +371,15 @@ impl TwoBodyConstraintBuilder {
         solved_dt: Real,
         rb1_pos: &Isometry<Real>,
         rb2_pos: &Isometry<Real>,
-        ccd_thickness: Real,
         constraint: &mut TwoBodyConstraint,
     ) {
-        let cfm_factor = params.cfm_factor();
+        let cfm_factor = params.contact_cfm_factor();
         let inv_dt = params.inv_dt();
-        let erp_inv_dt = params.erp_inv_dt();
+        let erp_inv_dt = params.contact_erp_inv_dt();
 
         let all_infos = &self.infos[..constraint.num_contacts as usize];
         let all_elements = &mut constraint.elements[..constraint.num_contacts as usize];
 
-        let mut is_fast_contact = false;
-
         #[cfg(feature = "dim2")]
         let tangents1 = constraint.dir1.orthonormal_basis();
         #[cfg(feature = "dim3")]
@@ -408,11 +397,9 @@ impl TwoBodyConstraintBuilder {
             // Normal part.
             {
                 let rhs_wo_bias = info.normal_rhs_wo_bias + dist.max(0.0) * inv_dt;
-                let rhs_bias = erp_inv_dt
-                    * (dist + params.allowed_linear_error())
-                        .clamp(-params.max_penetration_correction(), 0.0);
+                let rhs_bias = (erp_inv_dt * (dist + params.allowed_linear_error()))
+                    .clamp(-params.max_corrective_velocity(), 0.0);
                 let new_rhs = rhs_wo_bias + rhs_bias;
-                is_fast_contact = is_fast_contact || (-new_rhs * params.dt > ccd_thickness * 0.5);
 
                 element.normal_part.rhs_wo_bias = rhs_wo_bias;
                 element.normal_part.rhs = new_rhs;
@@ -432,7 +419,7 @@ impl TwoBodyConstraintBuilder {
             }
         }
 
-        constraint.cfm_factor = if is_fast_contact { 1.0 } else { cfm_factor };
+        constraint.cfm_factor = cfm_factor;
     }
 }
 
diff --git a/src/dynamics/solver/contact_constraint/two_body_constraint_simd.rs b/src/dynamics/solver/contact_constraint/two_body_constraint_simd.rs
index 4c4ef5229..2e2c68a55 100644
--- a/src/dynamics/solver/contact_constraint/two_body_constraint_simd.rs
+++ b/src/dynamics/solver/contact_constraint/two_body_constraint_simd.rs
@@ -15,7 +15,6 @@ use crate::math::{
 use crate::utils::SimdBasis;
 use crate::utils::{self, SimdAngularInertia, SimdCross, SimdDot};
 use num::Zero;
-use parry::math::SimdBool;
 use parry::utils::SdpMatrix2;
 use simba::simd::{SimdPartialOrd, SimdValue};
 
@@ -251,20 +250,16 @@ impl TwoBodyConstraintBuilderSimd {
         _multibodies: &MultibodyJointSet,
         constraint: &mut TwoBodyConstraintSimd,
     ) {
-        let cfm_factor = SimdReal::splat(params.cfm_factor());
-        let dt = SimdReal::splat(params.dt);
+        let cfm_factor = SimdReal::splat(params.contact_cfm_factor());
         let inv_dt = SimdReal::splat(params.inv_dt());
         let allowed_lin_err = SimdReal::splat(params.allowed_linear_error());
-        let erp_inv_dt = SimdReal::splat(params.erp_inv_dt());
-        let max_penetration_correction = SimdReal::splat(params.max_penetration_correction());
+        let erp_inv_dt = SimdReal::splat(params.contact_erp_inv_dt());
+        let max_corrective_velocity = SimdReal::splat(params.max_corrective_velocity());
         let warmstart_coeff = SimdReal::splat(params.warmstart_coefficient);
 
         let rb1 = gather![|ii| &bodies[constraint.solver_vel1[ii]]];
         let rb2 = gather![|ii| &bodies[constraint.solver_vel2[ii]]];
 
-        let ccd_thickness = SimdReal::from(gather![|ii| rb1[ii].ccd_thickness])
-            + SimdReal::from(gather![|ii| rb2[ii].ccd_thickness]);
-
         let poss1 = Isometry::from(gather![|ii| rb1[ii].position]);
         let poss2 = Isometry::from(gather![|ii| rb2[ii].position]);
 
@@ -279,7 +274,6 @@ impl TwoBodyConstraintBuilderSimd {
             constraint.dir1.cross(&constraint.tangent1),
         ];
 
-        let mut is_fast_contact = SimdBool::splat(false);
         let solved_dt = SimdReal::splat(solved_dt);
 
         for (info, element) in all_infos.iter().zip(all_elements.iter_mut()) {
@@ -292,12 +286,9 @@ impl TwoBodyConstraintBuilderSimd {
             {
                 let rhs_wo_bias =
                     info.normal_rhs_wo_bias + dist.simd_max(SimdReal::zero()) * inv_dt;
-                let rhs_bias = (dist + allowed_lin_err)
-                    .simd_clamp(-max_penetration_correction, SimdReal::zero())
-                    * erp_inv_dt;
+                let rhs_bias = ((dist + allowed_lin_err) * erp_inv_dt)
+                    .simd_clamp(-max_corrective_velocity, SimdReal::zero());
                 let new_rhs = rhs_wo_bias + rhs_bias;
-                is_fast_contact =
-                    is_fast_contact | (-new_rhs * dt).simd_gt(ccd_thickness * SimdReal::splat(0.5));
 
                 element.normal_part.rhs_wo_bias = rhs_wo_bias;
                 element.normal_part.rhs = new_rhs;
@@ -317,7 +308,7 @@ impl TwoBodyConstraintBuilderSimd {
             }
         }
 
-        constraint.cfm_factor = SimdReal::splat(1.0).select(is_fast_contact, cfm_factor);
+        constraint.cfm_factor = cfm_factor;
     }
 }
 
diff --git a/src/dynamics/solver/island_solver.rs b/src/dynamics/solver/island_solver.rs
index 159bfa7f4..718542b20 100644
--- a/src/dynamics/solver/island_solver.rs
+++ b/src/dynamics/solver/island_solver.rs
@@ -47,8 +47,6 @@ impl IslandSolver {
 
         let mut params = *base_params;
         params.dt /= num_solver_iterations as Real;
-        params.damping_ratio /= num_solver_iterations as Real;
-        // params.joint_damping_ratio /= num_solver_iterations as Real;
 
         /*
          *
diff --git a/src/dynamics/solver/velocity_solver.rs b/src/dynamics/solver/velocity_solver.rs
index 9d5c27980..a4e03b280 100644
--- a/src/dynamics/solver/velocity_solver.rs
+++ b/src/dynamics/solver/velocity_solver.rs
@@ -205,17 +205,19 @@ impl VelocitySolver {
             /*
              * Resolution without bias.
              */
-            for _ in 0..params.num_internal_stabilization_iterations {
-                joint_constraints
-                    .solve_wo_bias(&mut self.solver_vels, &mut self.generic_solver_vels);
-                contact_constraints.solve_restitution_wo_bias(
-                    &mut self.solver_vels,
-                    &mut self.generic_solver_vels,
-                );
-            }
+            if params.num_internal_stabilization_iterations > 0 {
+                for _ in 0..params.num_internal_stabilization_iterations {
+                    joint_constraints
+                        .solve_wo_bias(&mut self.solver_vels, &mut self.generic_solver_vels);
+                    contact_constraints.solve_restitution_wo_bias(
+                        &mut self.solver_vels,
+                        &mut self.generic_solver_vels,
+                    );
+                }
 
-            contact_constraints
-                .solve_friction(&mut self.solver_vels, &mut self.generic_solver_vels);
+                contact_constraints
+                    .solve_friction(&mut self.solver_vels, &mut self.generic_solver_vels);
+            }
         }
     }
 
diff --git a/src/lib.rs b/src/lib.rs
index 32e15411c..5c660ba60 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -15,6 +15,7 @@
 #![allow(clippy::too_many_arguments)]
 #![allow(clippy::needless_range_loop)] // TODO: remove this? I find that in the math code using indices adds clarity.
 #![allow(clippy::module_inception)]
+#![allow(unexpected_cfgs)] // This happens due to the dim2/dim3/f32/f64 cfg.
 
 #[cfg(all(feature = "dim2", feature = "f32"))]
 pub extern crate parry2d as parry;
diff --git a/src_testbed/lib.rs b/src_testbed/lib.rs
index 36c9199ae..2cb002ac3 100644
--- a/src_testbed/lib.rs
+++ b/src_testbed/lib.rs
@@ -1,4 +1,5 @@
 #![allow(clippy::too_many_arguments)]
+#![allow(unexpected_cfgs)] // This happens due to the dim2/dim3/f32/f64 cfg.
 
 extern crate nalgebra as na;
 
diff --git a/src_testbed/ui.rs b/src_testbed/ui.rs
index 7f3bdacab..aa9dc7a71 100644
--- a/src_testbed/ui.rs
+++ b/src_testbed/ui.rs
@@ -162,7 +162,7 @@ pub fn update_ui(
             ui.add(
                 Slider::new(
                     &mut integration_parameters.num_internal_stabilization_iterations,
-                    1..=100,
+                    0..=100,
                 )
                 .text("max internal stabilization iters."),
             );
@@ -170,12 +170,35 @@ pub fn update_ui(
                 Slider::new(&mut integration_parameters.warmstart_coefficient, 0.0..=1.0)
                     .text("warmstart coefficient"),
             );
-            ui.add(Slider::new(&mut integration_parameters.erp, 0.0..=1.0).text("erp"));
+
+            let mut substep_params = *integration_parameters;
+            substep_params.dt /= substep_params.num_solver_iterations.get() as Real;
+            let curr_erp = substep_params.contact_erp();
+            let curr_cfm_factor = substep_params.contact_cfm_factor();
+            ui.add(
+                Slider::new(
+                    &mut integration_parameters.contact_natural_frequency,
+                    0.01..=120.0,
+                )
+                .text(format!("contacts Hz (erp = {:.3})", curr_erp)),
+            );
+            ui.add(
+                Slider::new(
+                    &mut integration_parameters.contact_damping_ratio,
+                    0.01..=20.0,
+                )
+                .text(format!(
+                    "damping ratio (cfm-factor = {:.3})",
+                    curr_cfm_factor
+                )),
+            );
             ui.add(
-                Slider::new(&mut integration_parameters.damping_ratio, 0.0..=20.0)
-                    .text("damping ratio"),
+                Slider::new(
+                    &mut integration_parameters.joint_natural_frequency,
+                    0.0..=1200000.0,
+                )
+                .text("joint erp"),
             );
-            ui.add(Slider::new(&mut integration_parameters.joint_erp, 0.0..=1.0).text("joint erp"));
             ui.add(
                 Slider::new(&mut integration_parameters.joint_damping_ratio, 0.0..=20.0)
                     .text("joint damping ratio"),
@@ -197,6 +220,8 @@ pub fn update_ui(
             Slider::new(&mut integration_parameters.min_island_size, 1..=10_000)
                 .text("min island size"),
         );
+        ui.add(Slider::new(&mut state.nsteps, 1..=100).text("sims. per frame"));
+
         let mut frequency = integration_parameters.inv_dt().round() as u32;
         ui.add(Slider::new(&mut frequency, 0..=240).text("frequency (Hz)"));
         integration_parameters.set_inv_dt(frequency as Real);