From 6a34858ee28d72578035e140f4cfbcc1f36d078f Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Tue, 4 Aug 2026 16:28:43 -0500 Subject: [PATCH 1/5] Fill gaps and validate inputs in the Python potential bindings Several parts of the documented C++ API had no Python equivalent, which made the GCP and convergent-formulation tutorials impossible to follow from Python: - Add SmoothCollisions.compute_adaptive_dhat. Without it, adaptive dhat was unreachable from Python even though build() accepts use_adaptive_dhat=True and requires this to be called first. - Add SmoothContactParameters.adaptive_dhat_ratio property. - Add BarrierPotential.stiffness and .use_physical_barrier properties, mirroring set_stiffness()/set_use_physical_barrier() in C++. Rename the Python SmoothContactPotential class from "SmoothPotential" to match the C++ name. It had no in-tree users and the package is still a 2.0 alpha, so this is a straight rename with no alias. Validate preconditions in the bindings rather than relying on the C++ asserts. BarrierPotential asserts dhat > 0, stiffness > 0, and a non-null barrier, but assert() is compiled out under NDEBUG, so a release build would silently accept a bad value and produce undefined behavior. The bindings now raise ValueError, following the existing py::value_error convention in common.hpp. The new assert_positive helper is written as !(value > 0) so NaN is rejected as well. Co-Authored-By: Claude Opus 5 --- .../collisions/normal/normal_collisions.cpp | 15 +++++ python/src/common.hpp | 28 +++++++++ python/src/potentials/barrier_potential.cpp | 62 ++++++++++++++++--- 3 files changed, 95 insertions(+), 10 deletions(-) diff --git a/python/src/collisions/normal/normal_collisions.cpp b/python/src/collisions/normal/normal_collisions.cpp index 901959ab1..d394cbc98 100644 --- a/python/src/collisions/normal/normal_collisions.cpp +++ b/python/src/collisions/normal/normal_collisions.cpp @@ -19,6 +19,21 @@ void define_smooth_collisions(py::module_& m, const std::string& name) { py::class_(m, name.c_str()) .def(py::init()) + .def( + "compute_adaptive_dhat", &SmoothCollisions::compute_adaptive_dhat, + R"ipc_Qu8mg5v7( + Compute the per-element adaptive dhat from the rest configuration. + + Note: + Must be called before build() when using use_adaptive_dhat=True. + + Parameters: + mesh: The collision mesh. + vertices: Vertices of the collision mesh. + params: SmoothContactParameters. + broad_phase: Broad phase method. + )ipc_Qu8mg5v7", + "mesh"_a, "vertices"_a, "params"_a, "broad_phase"_a = nullptr) .def( "build", py::overload_cast< diff --git a/python/src/common.hpp b/python/src/common.hpp index 2fc728f4c..719cd8816 100644 --- a/python/src/common.hpp +++ b/python/src/common.hpp @@ -15,6 +15,34 @@ using namespace py::literals; #include #include +/// @brief Check that a parameter is strictly positive. +/// @throws py::value_error (Python ValueError) if value is not positive. +/// @note The C++ API only enforces this with an assert(), which is compiled +/// out under NDEBUG. Validating here means Python users get an +/// exception rather than undefined behavior in a release build. +/// @note Written as !(value > 0) so that NaN is rejected too. +inline void assert_positive(const double value, const std::string& name) +{ + if (!(value > 0)) { + throw py::value_error( + "Parameter " + name + " has invalid value: expected " + name + + " > 0 but got " + name + " = " + std::to_string(value)); + } +} + +/// @brief Check that a shared_ptr parameter is not None. +/// @throws py::value_error (Python ValueError) if ptr is null. +template +inline void +assert_not_none(const std::shared_ptr& ptr, const std::string& name) +{ + if (ptr == nullptr) { + throw py::value_error( + "Parameter " + name + " has invalid value: expected " + name + + " to be a valid object but got None"); + } +} + template void assert_2D_or_3D_vector( const Eigen::MatrixBase& v, const std::string& name) diff --git a/python/src/potentials/barrier_potential.cpp b/python/src/potentials/barrier_potential.cpp index 1d8f74561..25224fc9c 100644 --- a/python/src/potentials/barrier_potential.cpp +++ b/python/src/potentials/barrier_potential.cpp @@ -9,7 +9,12 @@ void define_barrier_potential(py::module_& m) { py::class_(m, "BarrierPotential") .def( - py::init(), + py::init([](const double dhat, const double stiffness, + const bool use_physical_barrier) { + assert_positive(dhat, "dhat"); + assert_positive(stiffness, "stiffness"); + return BarrierPotential(dhat, stiffness, use_physical_barrier); + }), R"ipc_Qu8mg5v7( Construct a barrier potential. @@ -17,12 +22,21 @@ void define_barrier_potential(py::module_& m) dhat: The activation distance of the barrier. stiffness: The stiffness of the barrier. use_physical_barrier: Whether to use the physical barrier. + + Raises: + ValueError: If dhat or stiffness is not positive. )ipc_Qu8mg5v7", "dhat"_a, "stiffness"_a, "use_physical_barrier"_a = false) .def( - py::init< - const std::shared_ptr, const double, const double, - const bool>(), + py::init([](std::shared_ptr barrier, const double dhat, + const double stiffness, + const bool use_physical_barrier) { + assert_not_none(barrier, "barrier"); + assert_positive(dhat, "dhat"); + assert_positive(stiffness, "stiffness"); + return BarrierPotential( + std::move(barrier), dhat, stiffness, use_physical_barrier); + }), R"ipc_Qu8mg5v7( Construct a barrier potential. @@ -31,18 +45,40 @@ void define_barrier_potential(py::module_& m) dhat: The activation distance of the barrier. stiffness: The stiffness of the barrier. use_physical_barrier: Whether to use the physical barrier. + + Raises: + ValueError: If barrier is None, or dhat or stiffness is not positive. )ipc_Qu8mg5v7", "barrier"_a, "dhat"_a, "stiffness"_a, "use_physical_barrier"_a = false) .def_property( - "dhat", &BarrierPotential::dhat, &BarrierPotential::set_dhat, - "Barrier activation distance.") + "dhat", &BarrierPotential::dhat, + [](BarrierPotential& self, const double dhat) { + assert_positive(dhat, "dhat"); + self.set_dhat(dhat); + }, + "Barrier activation distance. Must be positive.") + .def_property( + "stiffness", &BarrierPotential::stiffness, + [](BarrierPotential& self, const double stiffness) { + assert_positive(stiffness, "stiffness"); + self.set_stiffness(stiffness); + }, + "Barrier stiffness. Must be positive.") + .def_property( + "use_physical_barrier", &BarrierPotential::use_physical_barrier, + &BarrierPotential::set_use_physical_barrier, + "Whether to use the physical barrier.") .def_property( "barrier", py::cpp_function( &BarrierPotential::barrier, py::return_value_policy::reference), - &BarrierPotential::set_barrier, - "Barrier function used to compute the potential."); + [](BarrierPotential& self, + const std::shared_ptr& barrier) { + assert_not_none(barrier, "barrier"); + self.set_barrier(barrier); + }, + "Barrier function used to compute the potential. Must not be None."); } void define_smooth_potential(py::module_& m) @@ -73,9 +109,15 @@ void define_smooth_potential(py::module_& m) .def_readonly("beta_t", &SmoothContactParameters::beta_t) .def_readonly("alpha_n", &SmoothContactParameters::alpha_n) .def_readonly("beta_n", &SmoothContactParameters::beta_n) - .def_readonly("r", &SmoothContactParameters::r); + .def_readonly("r", &SmoothContactParameters::r) + .def_property( + "adaptive_dhat_ratio", + &SmoothContactParameters::adaptive_dhat_ratio, + &SmoothContactParameters::set_adaptive_dhat_ratio, + "Ratio of the distance to the interaction set in the rest " + "configuration used as the per-element adaptive dhat."); - py::class_(m, "SmoothPotential") + py::class_(m, "SmoothContactPotential") .def( py::init(), R"ipc_Qu8mg5v7( From b2e4c29975f27a9d9f08f3aa3cbcf99b15af2c70 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Tue, 4 Aug 2026 16:29:01 -0500 Subject: [PATCH 2/5] Update tutorials to match the current API Verified every snippet in docs/source/tutorials by extracting the C++ into a compile harness (-fsyntax-only against the real headers) and running the Python against a built ipctk. Both now pass end to end. Removed/renamed API the tutorials still used: - ipc::point_triangle_ccd and the other free narrow-phase functions are now methods on NarrowPhaseCCD subclasses; no longer exists. - The four *_nonlinear_ccd free functions are now NonlinearCCD methods. - CollisionStencil::ccd takes stencil vertices, not (vertices, edges, faces); use dof() to gather them. - TangentialCollisions::build no longer takes barrier_stiffness. In C++ the stale call still compiled, silently binding barrier_stiffness to mu_s and mu to mu_k. Stiffness now comes from the normal potential. - CollisionMesh gained an orient_vertex mask, so the 4-argument construct_is_on_surface form no longer compiles. - ProjectToPSD is now PSDProjectionMethod (and NONE was undocumented). - Candidates::build takes a BroadPhase*, so the C++ call needs &broad_phase. - ipctk.Collisions does not exist; rest_positions is a property, not a method; initial_barrier_stiffness returns max_barrier_stiffness instead of taking it. Also fixed code that never worked: two Python snippets were SyntaxErrors (multi-line assignment without parentheses), a missing semicolon and a stray one, MatrixXd where MatrixXi/MatrixXd was required, filib::Interval qualified as ipc::Interval, and various undefined or misspelled identifiers (mesh vs collision_mesh, collision vs collisions, map_displacement). Corrected the note on conservative CCD. TightInclusionCCD does not scale the returned TOI in the normal path; it inflates the minimum separation the query stops at, capped at 1e-4, and only scales the TOI in the fallback taken when that query returns a TOI below SMALL_TOI. Because the cap usually binds, changing conservative_rescaling often has no effect on the result at all, which the previous wording actively obscured. Co-Authored-By: Claude Opus 5 --- docs/source/tutorials/adhesion.rst | 10 +-- docs/source/tutorials/convergent.rst | 20 +++-- docs/source/tutorials/gcp.rst | 2 +- docs/source/tutorials/getting_started.rst | 100 +++++++++++++++------- docs/source/tutorials/nonlinear_ccd.rst | 24 +++--- docs/source/tutorials/ogc.rst | 32 +++---- docs/source/tutorials/simulation.rst | 66 +++++++++----- 7 files changed, 162 insertions(+), 92 deletions(-) diff --git a/docs/source/tutorials/adhesion.rst b/docs/source/tutorials/adhesion.rst index 29e4901bf..1cca5597d 100644 --- a/docs/source/tutorials/adhesion.rst +++ b/docs/source/tutorials/adhesion.rst @@ -47,7 +47,7 @@ We can build a normal adhesion potential object and compute the adhesion potenti const double Y = 1e3; const double eps_c = 0.5; - const ipc::NormalAdhesionPotential A_n(dhat_p, dhat_a, Y, eps_c) + const ipc::NormalAdhesionPotential A_n(dhat_p, dhat_a, Y, eps_c); double adhesion_potential = A_n(normal_collisions, collision_mesh, vertices); .. md-tab-item:: Python @@ -132,7 +132,7 @@ We can build a tangential adhesion potential object and compute the adhesion pot eps_a = 0.01 A_t = ipctk.TangentialAdhesionPotential(eps_a) - adhesion_potential = A_t(tangential_collisions, collision_mesh, displacement); + adhesion_potential = A_t(tangential_collisions, collision_mesh, displacement) Derivatives ^^^^^^^^^^^ @@ -177,8 +177,7 @@ Similar to the friction model (see, `friction + #include // ... @@ -565,10 +571,10 @@ The following example shows how to use the narrow phase to determine if a point Eigen::Vector3d t2_t1 = t2_t0; // triangle vertex 2 at t=1 double toi; // output time of impact - bool is_colliding = ipc::point_triangle_ccd( + bool is_colliding = ipc::TightInclusionCCD().point_triangle_ccd( p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1, toi); assert(is_colliding); - assert(abs(toi - 0.5) < 1e-8); + assert(toi <= 0.5); // conservative estimate of the exact TOI of 0.5 .. md-tab-item:: Python @@ -591,12 +597,35 @@ The following example shows how to use the narrow phase to determine if a point # returns a boolean indicating if the point is colliding with the triangle # and the time of impact (TOI) - is_colliding, toi = ipctk.point_triangle_ccd( + is_colliding, toi = ipctk.TightInclusionCCD().point_triangle_ccd( p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1) - assert(is_colliding) - assert(abs(toi - 0.5) < 1e-8) + assert is_colliding + assert toi <= 0.5 # conservative estimate of the exact TOI of 0.5 -Alternatively, the ``FaceVertexCandidate`` class contains a ``ccd`` function that can be used to determine if the face-vertex pairing is colliding: +.. note:: + The returned time of impact is a *conservative* under-estimate, so do not test + it for exact equality. For the query above the exact TOI is ``0.5``, but the + returned value is slightly less than that. + + ``TightInclusionCCD`` achieves this primarily by inflating the *distance* at + which the query stops rather than by scaling the resulting TOI. It runs the + narrow-phase query with a minimum separation of + + .. math:: + d_\text{min} + \min\left((1 - r)(d_0 - d_\text{min}),\ 10^{-4}\right), + + where :math:`r` is ``conservative_rescaling`` and :math:`d_0` is the distance + at :math:`t=0`. The TOI therefore comes back early by roughly that separation + divided by the distance travelled. Note the :math:`10^{-4}` cap: for + well-separated primitives it is what binds, so changing + ``conservative_rescaling`` has no effect on the result. + + The TOI itself is multiplied by ``conservative_rescaling`` only in a fallback + path, when the query above returns a TOI below + ``TightInclusionCCD::SMALL_TOI``; the query is then rerun with the true + :math:`d_\text{min}` and the result scaled to keep it away from zero. + +Alternatively, the ``FaceVertexCandidate`` class contains a ``ccd`` function that can be used to determine if the face-vertex pairing is colliding. It takes the *stencil* vertices (the four vertices of the face-vertex pair), which you can gather from the full vertex matrix using ``CollisionStencil::dof``: .. md-tab-set:: @@ -606,9 +635,13 @@ Alternatively, the ``FaceVertexCandidate`` class contains a ``ccd`` function tha ipc::FaceVertexCandidate candidate = ...; // face-vertex candidate + const Eigen::MatrixXi& edges = collision_mesh.edges(); + const Eigen::MatrixXi& faces = collision_mesh.faces(); + double toi; // output time of impact bool is_colliding = candidate.ccd( - vertices_t0, vertices_t1, collision_mesh.edges(), collision_mesh.faces(), toi); + candidate.dof(vertices_t0, edges, faces), + candidate.dof(vertices_t1, edges, faces), toi); .. md-tab-item:: Python @@ -616,12 +649,15 @@ Alternatively, the ``FaceVertexCandidate`` class contains a ``ccd`` function tha candidate = ... # face-vertex candidate + edges, faces = collision_mesh.edges, collision_mesh.faces + # returns a boolean indicating if the point is colliding with the triangle # and the time of impact (TOI) is_colliding, toi = candidate.ccd( - vertices_t0, vertices_t1, collision_mesh.edges, collision_mesh.faces) + candidate.dof(vertices_t0, edges, faces), + candidate.dof(vertices_t1, edges, faces)) -The same can be done for point-edge collisions using the ``point_edge_ccd`` function or ``EdgeVertexCandidate`` class and for edge-edge collisions using the ``edge_edge_ccd`` function or ``EdgeEdgeCandidate`` class. +The same can be done for point-edge collisions using the ``NarrowPhaseCCD::point_edge_ccd`` method or ``EdgeVertexCandidate`` class and for edge-edge collisions using the ``NarrowPhaseCCD::edge_edge_ccd`` method or ``EdgeEdgeCandidate`` class. .. _minimum-separation-ccd: @@ -644,7 +680,8 @@ To do this, we need to set the ``min_distance`` parameter when calling ``is_step Eigen::MatrixXd collision_free_vertices = (vertices_t1 - vertices_t0) * max_step_size + vertices_t0; assert(ipc::is_step_collision_free( - mesh, vertices_t0, collision_free_vertices, /*min_distance=*/1e-4)); + collision_mesh, vertices_t0, collision_free_vertices, + /*min_distance=*/1e-4)); .. md-tab-item:: Python @@ -653,7 +690,8 @@ To do this, we need to set the ``min_distance`` parameter when calling ``is_step max_step_size = ipctk.compute_collision_free_stepsize( collision_mesh, vertices_t0, vertices_t1, min_distance=1e-4) - collision_free_vertices = - (vertices_t1 - vertices_t0) * max_step_size + vertices_t0 - assert(ipctk.is_step_collision_free( - mesh, vertices_t0, collision_free_vertices, min_distance=1e-4)) + collision_free_vertices = ( + (vertices_t1 - vertices_t0) * max_step_size + vertices_t0) + assert ipctk.is_step_collision_free( + collision_mesh, vertices_t0, collision_free_vertices, + min_distance=1e-4) diff --git a/docs/source/tutorials/nonlinear_ccd.rst b/docs/source/tutorials/nonlinear_ccd.rst index 347b4a7e6..2c9b9eb4a 100644 --- a/docs/source/tutorials/nonlinear_ccd.rst +++ b/docs/source/tutorials/nonlinear_ccd.rst @@ -5,27 +5,27 @@ We also implement CCD of nonlinear trajectories (of linear geometry) using the m The method works by transforming the nonlinear trajectories into (adaptive) piecewise linear trajectories with an envelope/minimum separation around each piece, enclosing the nonlinear trajectory. The method then performs CCD on the piecewise linear trajectories to find the earliest time of impact. -We provide the following functions to perform nonlinear CCD: +Nonlinear CCD is provided by the ``NonlinearCCD`` class, which exposes the following methods: .. md-tab-set:: .. md-tab-item:: C++ - * :cpp:func:`ipc::point_point_nonlinear_ccd`, - * :cpp:func:`ipc::point_edge_nonlinear_ccd`, - * :cpp:func:`ipc::edge_edge_nonlinear_ccd`, and - * :cpp:func:`ipc::point_triangle_nonlinear_ccd`. + * :cpp:func:`ipc::NonlinearCCD::point_point_ccd`, + * :cpp:func:`ipc::NonlinearCCD::point_edge_ccd`, + * :cpp:func:`ipc::NonlinearCCD::edge_edge_ccd`, and + * :cpp:func:`ipc::NonlinearCCD::point_triangle_ccd`. - Each of these functions take as input a :cpp:class:`ipc::NonlinearTrajectory` object for the endpoints of the linear geometry. + Each of these methods take as input a :cpp:class:`ipc::NonlinearTrajectory` object for the endpoints of the linear geometry. .. md-tab-item:: Python - * :py:func:`ipctk.point_point_nonlinear_ccd`, - * :py:func:`ipctk.point_edge_nonlinear_ccd`, - * :py:func:`ipctk.edge_edge_nonlinear_ccd`, and - * :py:func:`ipctk.point_triangle_nonlinear_ccd`. + * :py:meth:`ipctk.NonlinearCCD.point_point_ccd`, + * :py:meth:`ipctk.NonlinearCCD.point_edge_ccd`, + * :py:meth:`ipctk.NonlinearCCD.edge_edge_ccd`, and + * :py:meth:`ipctk.NonlinearCCD.point_triangle_ccd`. - Each of these functions take as input a :py:class:`ipctk.NonlinearTrajectory` object for the endpoints of the linear geometry. + Each of these methods take as input a :py:class:`ipctk.NonlinearTrajectory` object for the endpoints of the linear geometry. For example, the following code defines a rigid trajectory in 2D in order to perform nonlinear CCD between a point and edge: @@ -207,7 +207,7 @@ The following code snippet shows an example of how to use interval arithmetic to Eigen::ConstRef center, Eigen::ConstRef point, const double omega, - const Interval& t) + const filib::Interval& t) { // 2×2 matrix of intervals representing the rotation matrix Matrix2I R; diff --git a/docs/source/tutorials/ogc.rst b/docs/source/tutorials/ogc.rst index 995cdfa14..04cd3e1cb 100644 --- a/docs/source/tutorials/ogc.rst +++ b/docs/source/tutorials/ogc.rst @@ -263,18 +263,18 @@ Putting it all together, a single simulation step using OGC looks like this: // 2. Warm start (Predict & Initialize) // x is current position, pred_x is x^t + dt * v^t trust_region.warm_start_time_step( - mesh, x, pred_x, collisions, dhat); + collision_mesh, x, pred_x, collisions, dhat); // 3. Solver Loop for (int i = 0; i < max_iterations; ++i) { // Update trust region if too many vertices hit the bound in previous step - trust_region.update_if_needed(mesh, x, collisions, dhat); + trust_region.update_if_needed(collision_mesh, x, collisions, dhat); // Compute search direction (Solver specific) Eigen::MatrixXd dx = compute_search_direction(x, ...); // Filter the step to respect OGC bounds - trust_region.filter_step(mesh, x, dx); + trust_region.filter_step(collision_mesh, x, dx); // Update positions x += dx; @@ -293,18 +293,18 @@ Putting it all together, a single simulation step using OGC looks like this: # 2. Warm start (Predict & Initialize) trust_region.warm_start_time_step( - mesh, x, pred_x, collisions, dhat) + collision_mesh, x, pred_x, collisions, dhat) # 3. Solver Loop for i in range(max_iterations): # Update trust region if needed - trust_region.update_if_needed(mesh, x, collisions, dhat) + trust_region.update_if_needed(collision_mesh, x, collisions, dhat) # Compute search direction (Solver specific) dx = compute_search_direction(x, ...) # Filter the step to respect OGC bounds - trust_region.filter_step(mesh, x, dx) + trust_region.filter_step(collision_mesh, x, dx) # Update positions x += dx @@ -416,9 +416,9 @@ Using ``planar_filter_step`` .. code-block:: c++ // Inside the solver loop, replace: - // trust_region.filter_step(mesh, x, dx); + // trust_region.filter_step(collision_mesh, x, dx); // with: - trust_region.planar_filter_step(mesh, x, dx); + trust_region.planar_filter_step(collision_mesh, x, dx); // Optionally tune the relaxation ratio via the struct member: // trust_region.relaxed_radius_scaling = 0.9; // default @@ -428,9 +428,9 @@ Using ``planar_filter_step`` .. code-block:: python # Inside the solver loop, replace: - # trust_region.filter_step(mesh, x, dx) + # trust_region.filter_step(collision_mesh, x, dx) # with: - trust_region.planar_filter_step(mesh, x, dx) + trust_region.planar_filter_step(collision_mesh, x, dx) # Optionally tune the relaxation ratio via the struct member: # trust_region.relaxed_radius_scaling = 0.9 # default @@ -451,18 +451,18 @@ Full Optimization Loop with Planar-DAT // 2. Warm start (Predict & Initialize) trust_region.warm_start_time_step( - mesh, x, pred_x, collisions, dhat); + collision_mesh, x, pred_x, collisions, dhat); // 3. Solver Loop for (int i = 0; i < max_iterations; ++i) { // Update trust region centers/radii if needed - trust_region.update_if_needed(mesh, x, collisions, dhat); + trust_region.update_if_needed(collision_mesh, x, collisions, dhat); // Compute search direction (Solver specific) Eigen::MatrixXd dx = compute_search_direction(x, ...); // Filter step using Planar-DAT (direction-aware truncation) - trust_region.planar_filter_step(mesh, x, dx); + trust_region.planar_filter_step(collision_mesh, x, dx); // Update positions x += dx; @@ -481,18 +481,18 @@ Full Optimization Loop with Planar-DAT # 2. Warm start (Predict & Initialize) trust_region.warm_start_time_step( - mesh, x, pred_x, collisions, dhat) + collision_mesh, x, pred_x, collisions, dhat) # 3. Solver Loop for i in range(max_iterations): # Update trust region centers/radii if needed - trust_region.update_if_needed(mesh, x, collisions, dhat) + trust_region.update_if_needed(collision_mesh, x, collisions, dhat) # Compute search direction (Solver specific) dx = compute_search_direction(x, ...) # Filter step using Planar-DAT (direction-aware truncation) - trust_region.planar_filter_step(mesh, x, dx) + trust_region.planar_filter_step(collision_mesh, x, dx) # Update positions x += dx diff --git a/docs/source/tutorials/simulation.rst b/docs/source/tutorials/simulation.rst index 7faf069be..ab2eb7954 100644 --- a/docs/source/tutorials/simulation.rst +++ b/docs/source/tutorials/simulation.rst @@ -38,7 +38,7 @@ From the full (volumetric) mesh vertices and surface edges/faces which index int // TODO: Show how to load a volumetric mesh from a file (e.g., using MshIO) // Faces of the surface mesh with indices into full_rest_positions - Eigen::MatrixXd faces; + Eigen::MatrixXi faces; igl::boundary_facets(tets, faces); // Edges of the surface mesh with indices into full_rest_positions @@ -71,7 +71,7 @@ This ``CollisionMesh`` can then be used just as any other ``CollisionMesh``. How .. code-block:: c++ // Convert full vertices to surface vertices - Eigen::VectorXd vertices = collision_mesh.vertices(full_vertices); + Eigen::MatrixXd vertices = collision_mesh.vertices(full_vertices); // Construct the set of collisions ipc::NormalCollisions collisions; @@ -84,7 +84,7 @@ This ``CollisionMesh`` can then be used just as any other ``CollisionMesh``. How double b = B(collisions, collision_mesh, vertices); // Convert full velocities to surface velocities - Eigen::VectorXd velocities = collision_mesh.map_displacements(full_velocities); + Eigen::MatrixXd velocities = collision_mesh.map_displacements(full_velocities); // Construct the set of friction collisions ipc::TangentialCollisions tangential_collisions; @@ -103,7 +103,7 @@ This ``CollisionMesh`` can then be used just as any other ``CollisionMesh``. How vertices = collision_mesh.vertices(full_vertices) # Construct the set of collisions - collisions = ipctk.Collisions() + collisions = ipctk.NormalCollisions() collisions.build(collision_mesh, vertices, dhat) # Construct a barrier potential @@ -144,12 +144,12 @@ When computing the gradient and Hessian of the potentials, the derivatives will .. code-block:: python - B = BarrierPotential(dhat, stiffness) + B = ipctk.BarrierPotential(dhat, stiffness) - grad = B.gradient(collision, collision_mesh, vertices) + grad = B.gradient(collisions, collision_mesh, vertices) grad_full = collision_mesh.to_full_dof(grad) - hess = B.hessian(collision, collision_mesh, vertices) + hess = B.hessian(collisions, collision_mesh, vertices) hess_full = collision_mesh.to_full_dof(hess) Codimensional Vertices @@ -170,9 +170,14 @@ In some cases, the collision mesh vertices are not the same as the surface verti std::vector is_on_surface = ipc::CollisionMesh::construct_is_on_surface( full_rest_positions.rows(), boundary_edges, codim_vertices); - // Construct the collision mesh from the is_on_surface vector and full mesh data + // is_orient_vertex marks the vertices with a well-defined orientation + // (i.e., those usable for signed distances). Codimensional vertices + // have no orientation, so false is a safe default for all vertices. + std::vector is_orient_vertex(full_rest_positions.rows(), false); + + // Construct the collision mesh from the masks and full mesh data ipc::CollisionMesh collision_mesh( - is_on_surface, full_rest_positions, edges, faces); + is_on_surface, is_orient_vertex, full_rest_positions, edges, faces); .. md-tab-item:: Python @@ -185,9 +190,14 @@ In some cases, the collision mesh vertices are not the same as the surface verti is_on_surface = ipctk.CollisionMesh.construct_is_on_surface( len(full_rest_positions), boundary_edges, codim_vertices) - # Construct the collision mesh from the is_on_surface vector and full mesh data + # is_orient_vertex marks the vertices with a well-defined orientation + # (i.e., those usable for signed distances). Codimensional vertices + # have no orientation, so False is a safe default for all vertices. + is_orient_vertex = [False] * len(full_rest_positions) + + # Construct the collision mesh from the masks and full mesh data collision_mesh = ipctk.CollisionMesh( - is_on_surface, full_rest_positions, edges, faces) + is_on_surface, is_orient_vertex, full_rest_positions, edges, faces) Nonlinear Bases and Curved Meshes ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -209,8 +219,8 @@ While IPC cannot directly handle nonlinear finite element bases and/or curved me Eigen::MatrixXd proxy_rest_positions; Eigen::MatrixXi proxy_edges, proxy_faces; // Load the proxy mesh from a file - igl::read_triangle_mesh("proxy.ply", rest_positions, faces); - igl::edges(faces, edges); + igl::read_triangle_mesh("proxy.ply", proxy_rest_positions, proxy_faces); + igl::edges(proxy_faces, proxy_edges); // Or build it from the volumetric mesh // Linear map from the finite element mesh to the collision proxy @@ -225,8 +235,8 @@ While IPC cannot directly handle nonlinear finite element bases and/or curved me # Finite element mesh fe_mesh = meshio.read("mesh.msh") - fe_rest_positions = mesh.points - tets = mesh.cells_dict["tetra"] + fe_rest_positions = fe_mesh.points + tets = fe_mesh.cells_dict["tetra"] # Collision proxy mesh # Load the proxy mesh from a file @@ -239,10 +249,10 @@ While IPC cannot directly handle nonlinear finite element bases and/or curved me # Linear map from the finite element mesh to the collision proxy displacement_map = ... # build or load the displacement map - collision_mesh = CollisionMesh( + collision_mesh = ipctk.CollisionMesh( proxy_rest_positions, proxy_edges, proxy_faces, displacement_map) -We can then map the displacements using ``collision_mesh.map_displacement(fe_displacements)`` or directly get the displaced proxy mesh vertices using ``collision_mesh.displace_vertices(fe_displacements)``. Similarly, we can map forces/potential gradients using ``collision_mesh.to_full_dof(collision_forces)`` or force Jacobians/potential Hessians using ``collision_mesh.to_full_dof(potential_hessian)``. +We can then map the displacements using ``collision_mesh.map_displacements(fe_displacements)`` or directly get the displaced proxy mesh vertices using ``collision_mesh.displace_vertices(fe_displacements)``. Similarly, we can map forces/potential gradients using ``collision_mesh.to_full_dof(collision_forces)`` or force Jacobians/potential Hessians using ``collision_mesh.to_full_dof(potential_hessian)``. .. warning:: The function ``CollisionMesh::vertices(full_positions)`` should not be used in this case because the rest positions used to construct the ``CollisionMesh`` are not the same as the finite element mesh's rest positions. Instead, use ``CollisionMesh::displace_vertices(fe_displacements)`` where ``fe_displacements`` is already the solution of the PDE or can be computed as ``fe_displacements = fe_positions - fe_rest_positions`` from deformed and rest positions. @@ -257,10 +267,24 @@ To remedy this, we can project the Hessian onto the positive semidefinite (PSD) .. md-tab-item:: C++ - - ``ProjectToPSD::CLAMP``: Clamp the negative eigenvalues of the Hessian to 0. This is the same as used by :cite:t:`Li2020IPC`. - - ``ProjectToPSD::ABS``: Set the negative eigenvalues of the Hessian to their absolute value. This is the method proposed by :cite:t:`Chen2024Stabler`. + - ``PSDProjectionMethod::NONE``: Do not project the Hessian. This is the default. + - ``PSDProjectionMethod::CLAMP``: Clamp the negative eigenvalues of the Hessian to 0. This is the same as used by :cite:t:`Li2020IPC`. + - ``PSDProjectionMethod::ABS``: Set the negative eigenvalues of the Hessian to their absolute value. This is the method proposed by :cite:t:`Chen2024Stabler`. + + .. code-block:: c++ + + Eigen::SparseMatrix hess = B.hessian( + collisions, collision_mesh, vertices, + ipc::PSDProjectionMethod::CLAMP); .. md-tab-item:: Python - - ``ProjectToPSD.CLAMP``: Clamp the negative eigenvalues of the Hessian to 0. This is the same as used by :cite:t:`Li2020IPC`. - - ``ProjectToPSD.ABS``: Set the negative eigenvalues of the Hessian to their absolute value. This is the method proposed by :cite:t:`Chen2024Stabler`. \ No newline at end of file + - ``PSDProjectionMethod.NONE``: Do not project the Hessian. This is the default. + - ``PSDProjectionMethod.CLAMP``: Clamp the negative eigenvalues of the Hessian to 0. This is the same as used by :cite:t:`Li2020IPC`. + - ``PSDProjectionMethod.ABS``: Set the negative eigenvalues of the Hessian to their absolute value. This is the method proposed by :cite:t:`Chen2024Stabler`. + + .. code-block:: python + + hess = B.hessian( + collisions, collision_mesh, vertices, + project_hessian_to_psd=ipctk.PSDProjectionMethod.CLAMP) \ No newline at end of file From f0711e561e618e987b450933ea7e7a9fcf07bc07 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Tue, 4 Aug 2026 17:19:43 -0500 Subject: [PATCH 3/5] Add Python tests for the potential bindings The input validation and newly exposed APIs added in 6a34858 had no automated coverage, so a regression would have been silent. Uses unittest.TestCase rather than plain functions so assertRaises is available under both nose2 (the CI runner) and pytest, without adding a pytest dependency. Covers: - BarrierPotential validation: ctor and setters reject <= 0 and NaN for dhat and stiffness, and None for barrier. Also asserts object state is unchanged after a rejected assignment, and that tiny-but-positive values still pass. - BarrierPotential.stiffness reaches the evaluation path, not just a stored field: tripling it triples the potential and gradient. - BarrierPotential.use_physical_barrier via the property is equivalent to the ctor kwarg for potential, gradient, and Hessian, plus a companion test that the flag changes the result at all so that equivalence is not vacuous. - SmoothCollisions.compute_adaptive_dhat as a differential pair: a baseline test pins that this mesh/dhat combination produces spurious nonzero forces at rest without adaptive dhat, and the adaptive test asserts they are exactly zero with it. The baseline is what keeps the second test meaningful. - SmoothContactParameters.adaptive_dhat_ratio round-trips and actually reaches compute_adaptive_dhat: larger ratios activate monotonically more collisions in a deformed configuration. - The SmoothContactPotential rename, guarded in both directions. Verified the tests bite by mutation testing: reverting the ctor validation, the setter validation, and the adaptive_dhat_ratio setter (to a no-op) turns them red, while the untouched barrier-ctor overload keeps passing, so the failures are specific rather than blanket. Potential/gradient/Hessian comparisons use a relative tolerance rather than exact equality, since the sums are parallel reductions whose operand order is not reproducible. Co-Authored-By: Claude Opus 5 --- python/tests/test_potentials.py | 264 ++++++++++++++++++++++++++++++++ 1 file changed, 264 insertions(+) create mode 100644 python/tests/test_potentials.py diff --git a/python/tests/test_potentials.py b/python/tests/test_potentials.py new file mode 100644 index 000000000..4c13e4900 --- /dev/null +++ b/python/tests/test_potentials.py @@ -0,0 +1,264 @@ +"""Tests for the potential bindings. + +Covers the Python-side input validation (which the C++ API only enforces with +assert(), compiled out under NDEBUG) and the smooth-contact/GCP bindings. + +Uses unittest.TestCase so that assertRaises is available under both nose2 (the +runner used in CI) and pytest, without adding a pytest dependency. +""" + +import unittest + +import numpy as np +from find_ipctk import ipctk +from utils import load_mesh + +# two-cubes-close.ply has a ~0.069 gap between the cubes, so dhat=0.1 activates +# a few hundred collisions while staying fast (~1 ms to build). +DHAT = 0.1 +STIFFNESS = 1.0 + + +def two_cubes(): + V, E, F = load_mesh("two-cubes-close.ply") + mesh = ipctk.CollisionMesh(V, E, F) + return mesh, mesh.rest_positions + + +def normal_collisions(mesh, vertices, dhat=DHAT): + collisions = ipctk.NormalCollisions() + collisions.build(mesh, vertices, dhat) + return collisions + + +class TestBarrierPotentialValidation(unittest.TestCase): + """The C++ ctor/setters assert dhat > 0, stiffness > 0, barrier != nullptr. + + assert() is compiled out under NDEBUG, so the bindings must validate and + raise ValueError instead of admitting undefined behavior in a release build. + """ + + INVALID = [0.0, -1.0, -1e-3, float("nan")] + + def test_ctor_rejects_invalid_dhat(self): + for dhat in self.INVALID: + with self.subTest(dhat=dhat): + with self.assertRaises(ValueError): + ipctk.BarrierPotential(dhat, STIFFNESS) + + def test_ctor_rejects_invalid_stiffness(self): + for stiffness in self.INVALID: + with self.subTest(stiffness=stiffness): + with self.assertRaises(ValueError): + ipctk.BarrierPotential(DHAT, stiffness) + + def test_ctor_rejects_none_barrier(self): + with self.assertRaises(ValueError): + ipctk.BarrierPotential(None, DHAT, STIFFNESS) + + def test_barrier_ctor_rejects_invalid_dhat(self): + for dhat in self.INVALID: + with self.subTest(dhat=dhat): + with self.assertRaises(ValueError): + ipctk.BarrierPotential( + ipctk.ClampedLogBarrier(), dhat, STIFFNESS) + + def test_dhat_setter_rejects_invalid(self): + for dhat in self.INVALID: + with self.subTest(dhat=dhat): + B = ipctk.BarrierPotential(DHAT, STIFFNESS) + with self.assertRaises(ValueError): + B.dhat = dhat + + def test_stiffness_setter_rejects_invalid(self): + for stiffness in self.INVALID: + with self.subTest(stiffness=stiffness): + B = ipctk.BarrierPotential(DHAT, STIFFNESS) + with self.assertRaises(ValueError): + B.stiffness = stiffness + + def test_barrier_setter_rejects_none(self): + B = ipctk.BarrierPotential(DHAT, STIFFNESS) + with self.assertRaises(ValueError): + B.barrier = None + + def test_state_unchanged_after_rejected_set(self): + """A rejected assignment must not partially apply.""" + B = ipctk.BarrierPotential(DHAT, STIFFNESS) + for attr, bad in (("dhat", 0.0), ("stiffness", -1.0)): + with self.subTest(attr=attr): + with self.assertRaises(ValueError): + setattr(B, attr, bad) + self.assertEqual(B.dhat, DHAT) + self.assertEqual(B.stiffness, STIFFNESS) + + def test_valid_values_accepted(self): + B = ipctk.BarrierPotential(DHAT, STIFFNESS) + B.dhat = 2e-3 + self.assertEqual(B.dhat, 2e-3) + B.stiffness = 5.0 + self.assertEqual(B.stiffness, 5.0) + B.barrier = ipctk.ClampedLogBarrier() + # use_physical_barrier has no precondition + B.use_physical_barrier = True + self.assertTrue(B.use_physical_barrier) + B.use_physical_barrier = False + self.assertFalse(B.use_physical_barrier) + # Tiny but positive must be allowed; only <= 0 and NaN are rejected. + ipctk.BarrierPotential(1e-300, 1e-300) + ipctk.BarrierPotential(ipctk.ClampedLogBarrier(), DHAT, STIFFNESS) + + +class TestBarrierPotentialProperties(unittest.TestCase): + """The stiffness/use_physical_barrier properties must reach the evaluation + path, not merely round-trip through a stored field.""" + + @classmethod + def setUpClass(cls): + cls.mesh, cls.vertices = two_cubes() + cls.collisions = normal_collisions(cls.mesh, cls.vertices) + assert len(cls.collisions) > 0, "fixture produced no collisions" + + def test_ctor_roundtrip(self): + B = ipctk.BarrierPotential(DHAT, 2.5, use_physical_barrier=True) + self.assertEqual(B.dhat, DHAT) + self.assertEqual(B.stiffness, 2.5) + self.assertTrue(B.use_physical_barrier) + + def test_stiffness_scales_potential_linearly(self): + """kappa multiplies the barrier potential, so tripling it must triple + the value. Guards against a setter that stores but is never read.""" + B1 = ipctk.BarrierPotential(DHAT, 1.0) + B3 = ipctk.BarrierPotential(DHAT, 1.0) + B3.stiffness = 3.0 + + p1 = B1(self.collisions, self.mesh, self.vertices) + p3 = B3(self.collisions, self.mesh, self.vertices) + self.assertGreater(p1, 0.0) + # rtol, not exact: the sum is a parallel reduction, so the operand + # order (and thus rounding) is not guaranteed to be reproducible. + np.testing.assert_allclose(p3, 3.0 * p1, rtol=1e-9) + + g1 = B1.gradient(self.collisions, self.mesh, self.vertices) + g3 = B3.gradient(self.collisions, self.mesh, self.vertices) + np.testing.assert_allclose(g3, 3.0 * g1, rtol=1e-9) + + def test_use_physical_barrier_setter_matches_ctor(self): + """Setting the property must be equivalent to passing the ctor kwarg.""" + for flag in (False, True): + with self.subTest(use_physical_barrier=flag): + via_ctor = ipctk.BarrierPotential( + DHAT, 2.5, use_physical_barrier=flag) + via_setter = ipctk.BarrierPotential(DHAT, 1.0) + via_setter.stiffness = 2.5 + via_setter.use_physical_barrier = flag + + args = (self.collisions, self.mesh, self.vertices) + np.testing.assert_allclose( + via_setter(*args), via_ctor(*args), rtol=1e-9) + np.testing.assert_allclose( + via_setter.gradient(*args), via_ctor.gradient(*args), + rtol=1e-9) + np.testing.assert_allclose( + via_setter.hessian(*args).todense(), + via_ctor.hessian(*args).todense(), rtol=1e-9) + + def test_use_physical_barrier_changes_result(self): + """The flag must actually alter the potential, otherwise the two + branches of the test above would agree trivially.""" + off = ipctk.BarrierPotential(DHAT, STIFFNESS, False) + on = ipctk.BarrierPotential(DHAT, STIFFNESS, True) + args = (self.collisions, self.mesh, self.vertices) + self.assertNotEqual(off(*args), on(*args)) + + +class TestSmoothContactParameters(unittest.TestCase): + def test_adaptive_dhat_ratio_default(self): + params = ipctk.SmoothContactParameters(DHAT, 0.5, 0.0, 0.1, 0.0, 2) + self.assertEqual(params.adaptive_dhat_ratio, 0.5) + + def test_adaptive_dhat_ratio_roundtrip(self): + params = ipctk.SmoothContactParameters(DHAT, 0.5, 0.0, 0.1, 0.0, 2) + for ratio in (0.1, 0.25, 0.9): + params.adaptive_dhat_ratio = ratio + self.assertEqual(params.adaptive_dhat_ratio, ratio) + + def test_adaptive_dhat_ratio_affects_adaptive_dhat(self): + """A larger ratio yields larger per-element dhat, so a deformed + configuration activates more collisions. Guards against the property + being stored but never reaching compute_adaptive_dhat().""" + mesh, rest = two_cubes() + + # Move the right cube toward the left one so the deformed state is + # close enough to activate, but only for large enough adaptive dhat. + deformed = rest.copy() + deformed[rest[:, 0] > 0.96, 0] -= 0.04 + + counts = [] + for ratio in (0.1, 0.5, 0.9): + params = ipctk.SmoothContactParameters(DHAT, 0.5, 0.0, 0.1, 0.0, 2) + params.adaptive_dhat_ratio = ratio + collisions = ipctk.SmoothCollisions() + collisions.compute_adaptive_dhat(mesh, rest, params) + collisions.build(mesh, deformed, params, True) + counts.append(len(collisions)) + + self.assertEqual(counts[0], 0, f"expected no activation, got {counts}") + self.assertLess(counts[0], counts[1], f"not monotonic: {counts}") + self.assertLess(counts[1], counts[2], f"not monotonic: {counts}") + + +class TestSmoothCollisionsAdaptiveDhat(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.mesh, cls.rest = two_cubes() + cls.params = ipctk.SmoothContactParameters( + DHAT, 0.5, 0.0, 0.1, 0.0, 2) + cls.potential = ipctk.SmoothContactPotential(cls.params) + + def _build(self, use_adaptive_dhat): + collisions = ipctk.SmoothCollisions() + if use_adaptive_dhat: + collisions.compute_adaptive_dhat(self.mesh, self.rest, self.params) + collisions.build(self.mesh, self.rest, self.params, use_adaptive_dhat) + return collisions + + def test_non_adaptive_has_spurious_rest_forces(self): + """Baseline: without adaptive dhat, this mesh/dhat pair produces a + nonzero potential at rest. This is what adaptive dhat exists to fix, so + if it ever becomes zero the test below stops proving anything.""" + collisions = self._build(use_adaptive_dhat=False) + potential = self.potential(collisions, self.mesh, self.rest) + gradient = self.potential.gradient(collisions, self.mesh, self.rest) + self.assertGreater(potential, 0.0) + self.assertGreater(np.abs(gradient).max(), 0.0) + + def test_adaptive_dhat_eliminates_spurious_rest_forces(self): + """GCP's 'no spurious forces' guarantee: with adaptive dhat computed + from the rest configuration, the potential and its gradient are exactly + zero in that configuration.""" + collisions = self._build(use_adaptive_dhat=True) + potential = self.potential(collisions, self.mesh, self.rest) + gradient = self.potential.gradient(collisions, self.mesh, self.rest) + self.assertEqual(potential, 0.0) + np.testing.assert_array_equal(gradient, np.zeros_like(gradient)) + + def test_broad_phase_argument_accepted(self): + collisions = ipctk.SmoothCollisions() + collisions.compute_adaptive_dhat( + self.mesh, self.rest, self.params, ipctk.LBVH()) + + +class TestSmoothContactPotentialNaming(unittest.TestCase): + """The Python class was previously exposed as "SmoothPotential", which did + not match the C++ name. Guard the rename in both directions.""" + + def test_matches_cpp_name(self): + self.assertTrue(hasattr(ipctk, "SmoothContactPotential")) + + def test_old_name_removed(self): + self.assertFalse(hasattr(ipctk, "SmoothPotential")) + + +if __name__ == "__main__": + unittest.main() From 98eb9218014656d0755d63c587eb2ecc54a8735d Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Tue, 4 Aug 2026 21:14:28 -0500 Subject: [PATCH 4/5] Describe the narrow-phase alternatives accurately MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Narrow-Phase section ended with "The alternatives are AdditiveCCD and InexactCCD", which was wrong in three ways: - InexactCCD is behind IPC_TOOLKIT_WITH_INEXACT_CCD, which defaults to OFF, so it does not exist in a default build and is absent from the Python module. Now marked opt-in, matching the wording already used in cpp-api/ccd.rst. - It implied a difference in conservatism policy between the three that does not exist. All three compute their margin as dmin + (1-r)(d0 - dmin); TightInclusionCCD alone caps the second term at 1e-4, which is the entire reason it reports a time of impact closer to the exact one for the same query. - It said nothing about AdditiveCCD's actual trade-off. For AdditiveCCD, lead with the strength (>100x faster, reliable in practice) and keep the theoretical caveat subordinate: it does not account for rounding error in its distance computations, but the default 10% margin is large enough to avoid false negatives, at the cost of a less accurate time of impact and more false positives. The failure mode is shrinking that margin, i.e. pushing conservative_rescaling toward 1.0 — not ordinary use. Also note that the margin is a fraction of the initial separation in excess of dmin rather than of the raw distance. That distinction comes from the identity documented on the gap computation in additive_ccd.cpp, (d - xi) = (d^2 - xi^2) / (d + xi), and it is not cosmetic: for a large minimum separation the two readings differ by an order of magnitude. Normalize d_\text{min} to d_\min, the convention already used elsewhere in the tutorials. Co-Authored-By: Claude Opus 5 --- docs/source/tutorials/getting_started.rst | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/docs/source/tutorials/getting_started.rst b/docs/source/tutorials/getting_started.rst index 32cce876a..efc83a52a 100644 --- a/docs/source/tutorials/getting_started.rst +++ b/docs/source/tutorials/getting_started.rst @@ -544,7 +544,19 @@ Possible values for ``broad_phase`` are: ``BruteForce`` (parallel brute force cu Narrow-Phase ^^^^^^^^^^^^ -The narrow phase computes the time of impact between two primitives (e.g., a point and a triangle or two edges in 3D). Narrow-phase algorithms are implemented as subclasses of ``NarrowPhaseCCD``. The default is ``TightInclusionCCD``, the Tight Inclusion CCD method of :cite:t:`Wang2021TightInclusion`, as it is provably conservative (i.e., never misses collisions), accurate (i.e., rarely reports false positives), and efficient. The alternatives are ``AdditiveCCD`` and ``InexactCCD``. +The narrow phase computes the time of impact between two primitives (e.g., a point and a triangle or two edges in 3D). Narrow-phase algorithms are implemented as subclasses of ``NarrowPhaseCCD``. The default is ``TightInclusionCCD``, the Tight Inclusion CCD method of :cite:t:`Wang2021TightInclusion`, as it is provably conservative (i.e., never misses collisions), accurate (i.e., rarely reports false positives), and efficient. It is what every function taking a ``narrow_phase_ccd`` argument uses unless you pass something else. + +Two other implementations are also available: + +- ``AdditiveCCD``, the method of :cite:t:`Li2021CIPC`, is much faster than Tight Inclusion (>100×) and reliable in practice. It does not account for rounding error in its distance computations, so in theory it can miss collisions, but its default margin (10% of the initial separation) is large enough to avoid this. The cost of that margin is a less accurate time of impact and more false positives. Shrinking it -- pushing ``conservative_rescaling`` toward ``1.0`` -- tightens the time of impact but is what can introduce false negatives. Tight Inclusion accounts for the rounding error, so it can reduce its tolerance without that trade-off :cite:p:`Belgrod2023Time`. +- ``InexactCCD``, the original method from the IPC codebase, is disabled by default. To use it, set the ``IPC_TOOLKIT_WITH_INEXACT_CCD`` CMake option to ``ON``. + +All three use the same expression for their conservative margin, stopping short of contact at a separation of + +.. math:: + d_\min + (1 - r)(d_0 - d_\min), + +where :math:`r` is ``conservative_rescaling`` and :math:`d_0` is the distance at :math:`t=0`. Note that the margin is a fraction of the initial separation *in excess of* :math:`d_\min`, not of the raw distance. ``TightInclusionCCD`` differs only in additionally capping the second term at :math:`10^{-4}` (see the note below), which is why it usually reports a time of impact much closer to the exact one than ``AdditiveCCD`` does for the same query. The following example shows how to use the narrow phase to determine if a point is colliding with a triangle (static in this case). @@ -612,7 +624,7 @@ The following example shows how to use the narrow phase to determine if a point narrow-phase query with a minimum separation of .. math:: - d_\text{min} + \min\left((1 - r)(d_0 - d_\text{min}),\ 10^{-4}\right), + d_\min + \min\left((1 - r)(d_0 - d_\min),\ 10^{-4}\right), where :math:`r` is ``conservative_rescaling`` and :math:`d_0` is the distance at :math:`t=0`. The TOI therefore comes back early by roughly that separation @@ -623,7 +635,7 @@ The following example shows how to use the narrow phase to determine if a point The TOI itself is multiplied by ``conservative_rescaling`` only in a fallback path, when the query above returns a TOI below ``TightInclusionCCD::SMALL_TOI``; the query is then rerun with the true - :math:`d_\text{min}` and the result scaled to keep it away from zero. + :math:`d_\min` and the result scaled to keep it away from zero. Alternatively, the ``FaceVertexCandidate`` class contains a ``ccd`` function that can be used to determine if the face-vertex pairing is colliding. It takes the *stencil* vertices (the four vertices of the face-vertex pair), which you can gather from the full vertex matrix using ``CollisionStencil::dof``: From f286192614d6e785cd886c2ff708c56f318775ce Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Tue, 4 Aug 2026 21:29:54 -0500 Subject: [PATCH 5/5] AdditiveCCD::additive_ccd use template for distance_squared functor parameter --- src/ipc/ccd/additive_ccd.cpp | 4 ++-- src/ipc/ccd/additive_ccd.hpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ipc/ccd/additive_ccd.cpp b/src/ipc/ccd/additive_ccd.cpp index b0e536ccb..0cd0abf21 100644 --- a/src/ipc/ccd/additive_ccd.cpp +++ b/src/ipc/ccd/additive_ccd.cpp @@ -68,11 +68,11 @@ AdditiveCCD::AdditiveCCD( conservative_rescaling = _conservative_rescaling; } +template bool AdditiveCCD::additive_ccd( VectorMax12d x, // mutable copy Eigen::ConstRef dx, - const std::function)>& - distance_squared, + const DistanceSqrFunc& distance_squared, const double max_disp_mag, double& toi, const double min_distance, diff --git a/src/ipc/ccd/additive_ccd.hpp b/src/ipc/ccd/additive_ccd.hpp index 875544c22..412758671 100644 --- a/src/ipc/ccd/additive_ccd.hpp +++ b/src/ipc/ccd/additive_ccd.hpp @@ -137,11 +137,11 @@ class AdditiveCCD : public NarrowPhaseCCD { /// @param min_distance The minimum distance between the objects. /// @param tmax The maximum time to check for collisions. /// @return True if a collision was detected, false otherwise. + template bool additive_ccd( VectorMax12d x, // mutable copy Eigen::ConstRef dx, - const std::function)>& - distance_squared, + const DistanceSqrFunc& distance_squared, const double max_disp_mag, double& toi, const double min_distance = 0.0,