Improve stiffness simulation

- Fix strange behavior with `RopeHandle`s
- Implement bidirectional stiffness computation
- Add more settings to fine-tune stiffness behavior

Author: mphe

demo/addons/ropesim/Rope.gd

@@ -4,6 +4,18 @@ class_name Rope
 
 # TODO: Split line rendering into a separate node
 
+enum StiffnessMethod {
+    ## Computes a single forward pass from front to back.
+    ## Produces decent results but may look unusual without [member fixate_begin] enabled.
+    ## It could also behave unstable in certain scenarios, i.e. jittery or erratic movement. [br]
+    ## This was the default behavior before v1.9.
+    Legacy,
+    ## Computes a forward and backward pass.
+    ## Produces overall better and more stable results than [member StiffnessMethod.Legacy].
+    ## Acts more sensitive to [member stiffness_bend_curve].
+    Bidirectional,
+}
+
 ## Triggered when the rope has been registered at the NativeRopeServer.
 signal on_registered()
 
@@ -37,6 +49,20 @@ signal on_point_count_changed()
 ## Might not produce 100% realistic results with fixed points.
 @export var stiffness: float = 0.0
 
+## Determines which method to use for stiffness simulation.
+@export var stiffness_method: StiffnessMethod = StiffnessMethod.Bidirectional
+
+## (Optional) Scale the stiffness along the rope.
+## If the rope behaves erratic,
+@export var stiffness_curve: Curve
+
+## (Optional) Scale the stiffness individually per segment according to how far each segment is bent.
+## The y value of the curve is multiplied with the stiffness.
+## x = 0 corresponds to no bending and x = 1 corresponds to 180° bending in regard to the previous
+## segment.
+## If unset, defaults to a quadratic curve.
+@export var stiffness_bend_curve: Curve
+
 ## Gravity
 @export var gravity: float = 100
 
@@ -450,7 +476,6 @@ func _set_seg_dist(value: Curve) -> void:
         segment_length_distribution.changed.connect(update_segments)
     update_segments()
 
-
 static func _resize_with_default(array: Variant, new_size: int, default: Variant) -> void:
     @warning_ignore("unsafe_method_access")
     var oldsize: int = array.size()

src/NativeRopeContext.cpp

@@ -10,7 +10,7 @@ const Vector2 VECTOR_ZERO = Vector2();   // NOLINT(cert-err58-cpp)
 const Vector2 VECTOR_DOWN = Vector2(0, 1);   // NOLINT(cert-err58-cpp)
 
 
-static float get_point_perc(int index, const PackedVector2Array& points)
+static float get_point_perc(int64_t index, const PackedVector2Array& points)
 {
     return points.size() > 0 ? static_cast<float>(index) / static_cast<float>(points.size() - 1) : 0;
 }
@@ -54,6 +54,9 @@ void NativeRopeContext::load_context(Node2D* rope)
     wind = rope->get("wind");
     damping = rope->get("damping");
     stiffness = rope->get("stiffness");
+    stiffness_method = static_cast<StiffnessMethod>(static_cast<int>(rope->get("stiffness_method")));
+    stiffness_curve = rope->get("stiffness_curve");
+    stiffness_bend_curve = rope->get("stiffness_bend_curve");
     max_endpoint_distance = rope->get("max_endpoint_distance");
     num_constraint_iterations = rope->get("num_constraint_iterations");
     seg_lengths = rope->call("get_segment_lengths");
@@ -132,7 +135,7 @@ void NativeRopeContext::_simulate_velocities(double delta)
 
     // Apply other forces
     _simulate_wind(&velocities);
-    _simulate_stiffness(&velocities);
+    _simulate_stiffness(delta, &velocities);
 
     // Move points according to velocity
     for (int i = first_idx; i < size; ++i)
@@ -158,46 +161,83 @@ void NativeRopeContext::_simulate_wind(PackedVector2Array* velocities) const
     }
 }
 
-void NativeRopeContext::_simulate_stiffness(PackedVector2Array* velocities) const
+void NativeRopeContext::_simulate_stiffness(double delta, PackedVector2Array* velocities) const
 {
-    // NOTE: oldpoints should not be used here, see comments in simulate_velocities().
-
     if (stiffness <= 0)
         return;
 
-    Vector2 parent_seg_dir = rope->get_global_transform().basis_xform(VECTOR_DOWN).normalized();
-    Vector2 last_stiffness_force;
+    const auto num_points = points.size();
+
+    if (stiffness_method == StiffnessMethod::Legacy)
+    {
+        const Vector2 parent_seg_dir = rope->get_global_transform().basis_xform(VECTOR_DOWN).normalized();
+        _simulate_stiffness_segment(delta, velocities, 1, num_points, parent_seg_dir);  // Start at second and stop at last point
+    }
+    else
+    {
+        const Vector2 end_parent = (points[num_points - 2] - points[num_points - 1]).normalized();
+
+        if (fixate_begin)
+        {
+            // TODO: Make direction configurable, e.g. an enum to select the mode ("transform", "gravity", "custom").
+            const Vector2 parent_seg_dir = rope->get_global_transform().basis_xform(VECTOR_DOWN).normalized();
+            _simulate_stiffness_segment(delta, velocities, 1, num_points, parent_seg_dir);  // Start at second and stop at last point
+            _simulate_stiffness_segment(delta, velocities, num_points - 2, 0, end_parent);  // Start at second last and stop before first point (because it's fixed)
+        }
+        else
+        {
+            const Vector2 begin_parent = (points[1] - points[0]).normalized();
+            _simulate_stiffness_segment(delta, velocities, 1, num_points, begin_parent);  // Start at second and stop at last point
+            _simulate_stiffness_segment(delta, velocities, num_points - 2, -1, end_parent);  // Start at second last and stop at first point
+        }
+    }
+}
+
+void NativeRopeContext::_simulate_stiffness_segment(double  /*delta*/, PackedVector2Array* velocities, int64_t idx_from, int64_t idx_stop, Vector2 last_stiffness_dir) const
+{
+    // NOTE: oldpoints should not be used here, see comments in simulate_velocities().
 
-    for (int i = 1; i < points.size(); ++i)
+    Vector2 force;
+    const bool use_stiffness_curve = stiffness_curve.is_valid() && stiffness_curve->get_point_count() > 0;
+    const bool use_bend_curve = stiffness_bend_curve.is_valid() && stiffness_bend_curve->get_point_count() > 0;
+    const int64_t idx_step = idx_from <= idx_stop ? 1 : -1;
+
+    for (int64_t i = idx_from; i != idx_stop; i += idx_step)
     {
         // NOTE: Asked a physicist to confirm this computation is physically accurate.
         // He mentioned that, while it is technically correct, there is a material-dependent limit
         // how far an object can bend before bending properties (stiffness) changes.
         // E.g. a material might be less inclined to snap back into place at smaller bend angles, or
         // a material might stop bending at some point, i.e. when it breaks.
-        // This implementation should probably suffice in most cases, but a more advanced
-        // implementation would include curves for stiffness in relation to segment position and
-        // for stiffness in relation to bend angle.
         //
         //  |  parent_seg_dir     --->  parent_seg_dir.orthogonal()
         //  |                     \
         //  V                      \   seg_dir
         //  \  seg_dir              V
-        //   \
-        //    V
-        const Vector2 seg_dir = (points[i] - points[i - 1]).normalized();
-        const float angle = seg_dir.angle_to(parent_seg_dir);
+        //   \                     / force_dir
+        //    V                   V
+        //   / force_dir
+        //  V
+        //
+        const auto last_idx = i - idx_step;
+        const Vector2 seg_dir = (points[i] - points[last_idx]).normalized();
+        const float angle = seg_dir.angle_to(last_stiffness_dir);  // angle is signed and can be used to determine the force direction
+        const float bend_factor = Math::absf(angle / 3.1415f);
+        const float scale_curve_pos = get_point_perc(idx_step >= 0 ? i : (points.size() - 1 - i), points);
+        const float scale = use_stiffness_curve ? stiffness_curve->sample_baked(scale_curve_pos) : 1.0f;
+        // If no curve specified, use quadratic scale by default
+        const float bend_scale = (use_bend_curve ? stiffness_bend_curve->sample_baked(bend_factor) : bend_factor * bend_factor);
 
         // The force directs orthogonal to the current segment
-        const Vector2 force_dir = seg_dir.orthogonal();
+        const Vector2 force_dir = Math::sign(-angle) * seg_dir.orthogonal();
+        force += force_dir * stiffness * scale * bend_scale;
 
-        // Scale the force the further the segment bends.
-        // angle is signed and can be used to determine the force direction
-        last_stiffness_force += force_dir * (-angle / 3.1415f) * stiffness;
-        parent_seg_dir = seg_dir;
+        // `simulation_weights` will be applied later in `_simulate_velocities()`
+        (*velocities)[i] += force;
 
-        // Update velocity
-        (*velocities)[i] += last_stiffness_force;
+        // Carry velocity over to the next segment. `simulation_weights` must be applied now to not carry over a wrong value.
+        force *= simulation_weights[i];
+        last_stiffness_dir = seg_dir;
     }
 }
 

src/NativeRopeContext.hpp

@@ -7,6 +7,13 @@
 
 namespace godot {
 
+// Should match Rope.gd
+enum class StiffnessMethod : int8_t
+{
+    Legacy,
+    Bidirectional
+};
+
 // Caches properties of a rope node and performs simulation.
 // TODO: Could be used as base class in the future to manage rope data in C++ and not in GDScript.
 class NativeRopeContext  // NOLINT(cppcoreguidelines-special-member-functions)
@@ -22,7 +29,8 @@ class NativeRopeContext  // NOLINT(cppcoreguidelines-special-member-functions)
 protected:
     void _simulate_velocities(double delta);
     void _simulate_wind(PackedVector2Array* velocities) const;
-    void _simulate_stiffness(PackedVector2Array* velocities) const;
+    void _simulate_stiffness(double delta, PackedVector2Array* velocities) const;
+    void _simulate_stiffness_segment(double delta, PackedVector2Array* velocities, int64_t idx_from, int64_t idx_stop, Vector2 last_stiffness_dir) const;
     void _resolve_collisions(double delta, bool disable_contact_reporting);
     void _constraint(double delta);
     void _writeback();
@@ -45,13 +53,16 @@ class NativeRopeContext  // NOLINT(cppcoreguidelines-special-member-functions)
     float damping = 0.0;
     Ref<Curve> damping_curve;
     float stiffness = 0.0;
+    Ref<Curve> stiffness_curve;
+    Ref<Curve> stiffness_bend_curve;
     float max_endpoint_distance = 0.0;
     int num_constraint_iterations = 0;
 
     float collision_radius = 1.0;
     float collision_damping = 0.0;
     int collision_mask = 0;
 
+    StiffnessMethod stiffness_method = StiffnessMethod::Bidirectional;
     bool fixate_begin = true;
     bool resolve_to_begin = false;
     bool resolve_to_end = false;