source/blender/geometry/CMakeLists.txt

@@ -17,6 +17,7 @@ set(INC
 set(SRC
   intern/add_curves_on_mesh.cc
   intern/constraint_solver.cc
+  intern/curve_constraint_solver.cc
   intern/fillet_curves.cc
   intern/mesh_merge_by_distance.cc
   intern/mesh_primitive_cuboid.cc
@@ -34,6 +35,7 @@ set(SRC
 
   GEO_add_curves_on_mesh.hh
   GEO_constraint_solver.hh
+  GEO_curve_constraint_solver.hh
   GEO_fillet_curves.hh
   GEO_mesh_merge_by_distance.hh
   GEO_mesh_primitive_cuboid.hh

source/blender/geometry/GEO_constraint_solver.hh

@@ -24,12 +24,6 @@ class ConstraintSolver {
      * For details see "Fast Simulation of Inextensible Hair and Fur"
      * by Matthias Mueller and Tae Yong Kim. */
     Sequential,
-    /* Position Based Dynamics (PBD) solves constraints based on relative mass.
-     * The solver requires multiple iterations per step. This is generally slower
-     * than the FTL method, but leads to physically correct movement.
-     *
-     * Based on "XPBD: Position-Based Simulation of Compliant Constrained Dynamics" */
-    PositionBasedDynamics,
   };
 
   struct Params {

source/blender/geometry/GEO_curve_constraint_solver.hh

@@ -0,0 +1,26 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+#include "BKE_curves.hh"
+
+namespace blender::geometry::curve_constraing_solver {
+
+void compute_segment_lengths(const bke::CurvesGeometry &curves,
+                             IndexMask curve_selection,
+                             MutableSpan<float> r_segment_lengths);
+
+void solve_length_constraints(const bke::CurvesGeometry &curves,
+                              IndexMask curve_selection,
+                              Span<float> segment_lenghts,
+                              MutableSpan<float3> positions);
+
+void solve_length_and_collision_constraints(const bke::CurvesGeometry &curves,
+                                            IndexMask curve_selection,
+                                            Span<float> segment_lengths,
+                                            Span<float3> start_positions,
+                                            const Mesh &surface,
+                                            const bke::CurvesSurfaceTransforms &transforms,
+                                            MutableSpan<float3> positions);
+
+}  // namespace blender::geometry::curve_constraing_solver

source/blender/geometry/intern/constraint_solver.cc

@@ -313,8 +313,6 @@ void ConstraintSolver::solve_constraints(CurvesGeometry &curves,
     switch (params_.solver_type) {
       case SolverType::Sequential:
         return 1;
-      case SolverType::PositionBasedDynamics:
-        return params_.max_solver_iterations;
     }
     return 0;
   }();
@@ -349,9 +347,6 @@ void ConstraintSolver::solve_curve_constraints(CurvesGeometry &curves,
       /* The sequential solver only moves the 2nd point of each segment. */
       case SolverType::Sequential:
         return params_.use_root_constraints ? 1 : 0;
-      /* The PBD solver moves both points, except for the pinned root. */
-      case SolverType::PositionBasedDynamics:
-        return params_.use_root_constraints ? 1 : 0;
     }
     return 0;
   }();

source/blender/geometry/intern/curve_constraint_solver.cc

@@ -0,0 +1,62 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "GEO_curve_constraint_solver.hh"
+
+namespace blender::geometry::curve_constraing_solver {
+
+void compute_segment_lengths(const bke::CurvesGeometry &curves,
+                             const IndexMask curve_selection,
+                             MutableSpan<float> r_segment_lengths)
+{
+  BLI_assert(r_segment_lengths.size() == curves.points_num());
+
+  const Span<float3> positions = curves.positions();
+  const OffsetIndices points_by_curve = curves.points_by_curve();
+  threading::parallel_for(curve_selection.index_range(), 256, [&](const IndexRange range) {
+    for (const int curve_i : curve_selection.slice(range)) {
+      const IndexRange points = points_by_curve[curve_i].drop_back(1);
+      for (const int point_i : points) {
+        const float3 &p1 = positions[point_i];
+        const float3 &p2 = positions[point_i + 1];
+        const float length = math::distance(p1, p2);
+        r_segment_lengths[point_i] = length;
+      }
+    }
+  });
+}
+
+void solve_length_constraints(const bke::CurvesGeometry &curves,
+                              const IndexMask curve_selection,
+                              const Span<float> segment_lenghts,
+                              MutableSpan<float3> positions)
+{
+  BLI_assert(segment_lenghts.size() == curves.points_num());
+
+  const OffsetIndices points_by_curve = curves.points_by_curve();
+  threading::parallel_for(curve_selection.index_range(), 256, [&](const IndexRange range) {
+    for (const int curve_i : curve_selection.slice(range)) {
+      const IndexRange points = points_by_curve[curve_i].drop_back(1);
+      for (const int point_i : points) {
+        const float3 &p1 = positions[point_i];
+        float3 &p2 = positions[point_i + 1];
+        const float3 direction = math::normalize(p2 - p1);
+        const float goal_length = segment_lenghts[point_i];
+        p2 = p1 + direction * goal_length;
+      }
+    }
+  });
+}
+
+void solve_length_and_collision_constraints(const bke::CurvesGeometry &curves,
+                                            const IndexMask curve_selection,
+                                            const Span<float> segment_lengths,
+                                            const Span<float3> start_positions,
+                                            const Mesh &surface,
+                                            const bke::CurvesSurfaceTransforms &transforms,
+                                            MutableSpan<float3> positions)
+{
+  UNUSED_VARS(
+      curves, curve_selection, segment_lengths, start_positions, surface, transforms, positions);
+}
+
+}  // namespace blender::geometry::curve_constraing_solver