Repository: ipc-sim/ipc-toolkit (language: python)

Issue:
Unexpectedly slow CCD
Hi, it's me again :wave: 

I noticed that CCD sometimes takes a very long time relative to other components of my simulator. For example, for a small impact problem, with a block dropping (with some initial velocity) on top of another, CCD can take ~20 seconds or more where factorization takes milliseconds. I can reproduce similar behavior for a very simple example, see [this Gist here](https://gist.github.com/Andlon/34d9d869d5267d28d86fb36f35f72538), where a single tetrahedron is squeezed so that the top vertex gets within `dhat` of the lower face.  Here I'm getting timings like this:

```
Potential duration: 6.1011e-05 seconds.
CCD duration: 2.25247 seconds.
```

Here's the most crucial part of the code (see the above gist for the full setup):
```c++
    // Run this several times just to demonstrate that it's not due to some initialization or similar
    for (int i = 0; i < 5; ++i) {
        const auto potential_begin = std::chrono::steady_clock::now();
        const double potential = constraints.compute_potential(mesh, deformed_vertices, dhat);
        const Eigen::VectorXd grad = constraints.compute_potential_gradient(mesh, deformed_vertices, dhat);
        const Eigen::SparseMatrix<double> hessian = constraints.compute_potential_hessian(mesh, deformed_vertices, dhat,project_spd);
        const auto potential_end = std::chrono::steady_clock::now();
        const auto potential_duration = std::chrono::duration<double>(potential_end - potential_begin).count();

        const auto ccd_begin = std::chrono::steady_clock::now();
        const double alpha = ipc::compute_collision_free_stepsize(mesh, rest_vertices, deformed_vertices);
        const auto ccd_end = std::chrono::steady_clock::now();
        const auto ccd_duration = std::chrono::duration<double>(ccd_end - ccd_begin).count();

        std::cout << "Potential duration: " << potential_duration << " seconds." << std::endl;
        std::cout << "CCD duration: " << ccd_duration << " seconds." << std::endl;
    }

```

Build type is `Release`. I don't have all that much experience with CCD, but this seems slower than what I was expecting. Is there something I should configure that I haven't? I've basically just followed the docs/tutorial here. Any hints as to what I can do here, if anything, would be much appreciated :-)


Discussion / hints:
(I had to choose, so I labeled this as a bug, although I'm not sure it is one)
(Also, for the record, I'm aware that the Hessian computed here would have to be passed through `to_full_dof` in order to be consistent with the call to `build_from_full_mesh`.)
It seems like you have setup up everything correctly, but I did notice one thing. Here
```
const double alpha = ipc::compute_collision_free_stepsize(mesh, rest_vertices, deformed_vertices);
```
you are doing CCD starting at the rest positions and going to the deformed positions. This might be what you intended, but normally we perform CCD between successive deformed positions in the Newton optimization for example.

What usually affects CCD performance is having tiny initial distances and/or large trajectories. The former causes our algorithm to refine further to avoid producing a conservative estimate of 0 for the TOI. The latter can also cause larger running times as it has to refine further to match our length-based convergence criteria.

When I have some free time, I will try to set up a unit test from your gist and see if there are any obvious bottlenecks that can be improved.
Looking at your gist it doesn't seem to be either of those cases, so it should have been able to get you a quick answer. I will definitely look more into it later.
> you are doing CCD starting at the rest positions and going to the deformed positions. This might be what you intended, but normally we perform CCD between successive deformed positions in the Newton optimization for example.

Yeah, that was just out of convenience for this example. In my simulator I'm indeed incorporating CCD into the line search.

Fix patch:
diff --git a/src/ipc/ccd/ccd.cpp b/src/ipc/ccd/ccd.cpp
index 2edc5c249..eaf811caa 100644
--- a/src/ipc/ccd/ccd.cpp
+++ b/src/ipc/ccd/ccd.cpp
@@ -9,6 +9,7 @@
 
 #ifdef IPC_TOOLKIT_WITH_CORRECT_CCD
 #include <tight_inclusion/ccd.hpp>
+#include <tight_inclusion/interval_root_finder.hpp>
 #else
 #include <CTCD.h>
 #endif
@@ -32,25 +33,27 @@ bool ccd_strategy(
         double /*min_distance*/,
         bool /*no_zero_toi*/,
         double& /*toi*/)>& ccd,
-    const double max_iterations,
+    const long max_iterations,
     const double min_distance,
     const double initial_distance,
     const double conservative_rescaling,
     double& toi)
 {
-
-    if (initial_distance == 0) {
-        logger().warn("Initial distance is 0, returning toi=0!");
+    if (initial_distance <= min_distance) {
+        logger().warn(
+            "Initial distance {} ≤ d_min={}, returning toi=0!",
+            initial_distance, min_distance);
         toi = 0;
         return true;
     }
 
-    const double min_effective_distance =
-        min_distance + (1.0 - conservative_rescaling) * initial_distance;
-    // #ifdef IPC_TOOLKIT_WITH_CORRECT_CCD
+    double min_effective_distance =
+        (1.0 - conservative_rescaling) * initial_distance;
+#ifdef IPC_TOOLKIT_WITH_CORRECT_CCD
     // Tight Inclusion performs better when the minimum separation is small
-    // min_distance = std::min(min_distance, 1e-4);
-    // #endif
+    min_effective_distance = std::min(min_effective_distance, 1e-4);
+#endif
+    min_effective_distance += min_distance;
 
     assert(min_effective_distance < initial_distance);
 


Test patch:
diff --git a/tests/ccd/test_ccd.cpp b/tests/ccd/test_ccd.cpp
index c61af3d19..bdeadb388 100644
--- a/tests/ccd/test_ccd.cpp
+++ b/tests/ccd/test_ccd.cpp
@@ -637,3 +637,43 @@ TEST_CASE("Slow EE CCD 2", "[ccd][edge-edge][slow][thisone]")
     CAPTURE(toi);
     CHECK(is_impacting);
 }
+
+TEST_CASE("Squash Tet", "[ccd]")
+{
+    const double dhat = 1e-3;
+
+    Eigen::MatrixXd rest_vertices(4, 3);
+    rest_vertices << 0.0, 0.0, 0.0, //
+        1.0, 0.0, 0.0,              //
+        0.0, 1.0, 0.0,              //
+        0.0, 0.0, 1.0;
+
+    Eigen::MatrixXd deformed_vertices = rest_vertices;
+    deformed_vertices(3, 0) += 0.1;
+    deformed_vertices(3, 1) -= 0.1;
+    deformed_vertices(3, 2) = -0.5 * dhat;
+
+    Eigen::MatrixXi edges(6, 2);
+    edges << 0, 1, //
+        0, 2,      //
+        0, 3,      //
+        1, 2,      //
+        1, 3,      //
+        2, 3;
+    Eigen::MatrixXi faces(4, 3);
+    faces << 0, 2, 1, //
+        0, 1, 3,      //
+        0, 3, 2,      //
+        1, 2, 3;
+
+    ipc::CollisionMesh mesh =
+        ipc::CollisionMesh::build_from_full_mesh(rest_vertices, edges, faces);
+
+    // BENCHMARK("compute toi")
+    // {
+    logger().debug(
+        "toi={}",
+        ipc::compute_collision_free_stepsize(
+            mesh, rest_vertices, deformed_vertices));
+    // };
+}
\ No newline at end of file