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

Issue:
Co-dimensional simulations
Quick question: Does this toolkit include the ability to handle co-dimensional simulations similar to Co-IPC? 



Discussion / hints:
I am not sure if you mean can the toolkit completely perform codimensional FEM simulations like in the CIPC paper, or if the toolkit has the necessary components to implement CIPC?

This code base includes the necessary components in order to add the IPC's frictional contact model to a simulation pipeline. This includes computing the barrier potentials, smooth friction dissipative potential, and performing CCD. What the toolkit does not have is the FEM components necessary to implement codimensional FEM (FEM basis, elasticity models, and optimization framework for time-stepping). 

That being said, this toolkit contains most of the components shown in the CIPC paper. 

## Implemented

These features from the CIPC paper are already implemented.

* Support for codimensional triangle-triangle, triangle-edge, triangle-point, edge-edge, and edge-point (in 2D) frictional contact.
    * This is sufficient for standard IPC and a lot of CIPC, but misses some types of contact (see below)
* The IPC thickness model (i.e., $\xi$ offset in the barrier potential)
* CCD between codimensional triangle-triangle, triangle-edge, triangle-point, edge-edge, and edge-point (in 2D)
    * This is using the [Tight-Inclusion](https://github.com/Continuous-Collision-Detection/Tight-Inclusion) CCD method for guaranteed conservative CCD under floating-point 
* Thickness offset in the CCD 
    * This is handled through Tight-Inclusion's minimum separation CCD.

## Missing Features

The following components are missing, but could be easily added (I am happy to add them).
* Support for codimensional edge-point (in 3D) and point-point contact. 
    * These are important for contacts between codim. edges with codim. points and for contacts between particles (like in Fig. 22 "Granules on cloth" of the CIPC paper)
* Strain limiting
    * This is more on the side of FEM but could be added to the toolkit. My only concern is it is a little out of scope.
* CCD for contacts between edge-point (in 3D) and point-point contact
    * I believe we can add these using Tight-Inclusion CCD
* Additive CCD
    *  Instead we choose to use Tight-Inclusion because of its guarantees

Thanks so much. My main interest is having a more robust discrete elastic rod framework than is currently implemented in Co-IPC. Eventually, I would want to be able to simulate rods interacting with a deformable solid geometry. Co-IPC gets 75% of the way there as currently written, but the implementation of rod physics is limited (only considers bending and rods that are straight in their stress-free configuration). I'm not a developer, so I'm just trying to determine the path of least resistance given my skills. Any plans to bring discrete elastic rods to polyfem (hint, hint). ;-) Seriously though, your work is amazingly impressive and appreciated! Thank you for bringing these tools to the public.
Thank you for your nice comments. Adding shell and rod models to PolyFEM is something I have wanted to play with for a while, but I don't have a research project that motivates implementing it. I think it will eventually be implemented, but no timeline right now.
Small update: #22 added the thickness offset in the CCD.

Fix patch:
diff --git a/.github/workflows/continuous.yml b/.github/workflows/continuous.yml
index a98802c70..f2076fce6 100644
--- a/.github/workflows/continuous.yml
+++ b/.github/workflows/continuous.yml
@@ -16,7 +16,7 @@ env:
   CTEST_PARALLEL_LEVEL: 2
 
 jobs:
-  Unix:
+  Build:
     name: ${{ matrix.name }} (${{ matrix.config }})
     runs-on: ${{ matrix.os }}
     strategy:
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 31cf2d68f..09d852ff4 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -66,6 +66,5 @@ jobs:
       - name: Deploy
         uses: JamesIves/github-pages-deploy-action@v4.4.1
         with:
-          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-          BRANCH: gh-pages # The branch the action should deploy to.
-          FOLDER: docs/build/html # The folder the action should deploy.
+          branch: gh-pages # The branch the action should deploy to.
+          folder: docs/build/html # The folder the action should deploy.
diff --git a/.gitignore b/.gitignore
index deaadbac9..fe7342623 100644
--- a/.gitignore
+++ b/.gitignore
@@ -243,6 +243,7 @@ core*
 docs/source/doxyoutput
 docs/source/api
 !docs/Makefile
+tests/data/ccd-queries
 tests/data/slow-broadphase-ccd
 tests/data/ccd-failure
 tests/data/friction/chain
diff --git a/CMakeLists.txt b/CMakeLists.txt
index f70473ffe..106e85440 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -70,6 +70,13 @@ option(IPC_TOOLKIT_WITH_RATIONAL_INTERSECTION "Use rational edge-triangle inters
 option(IPC_TOOLKIT_WITH_ROBIN_MAP             "Use Tessil's robin-map rather than std maps"    ON)
 option(IPC_TOOLKIT_WITH_ABSEIL                "Use Abseil's hash functions"                    ON)
 
+include(CMakeDependentOption)
+cmake_dependent_option(IPC_TOOLKIT_TEST_CCD_BENCHMARK "Enable CCD benchmark test" OFF "IPC_TOOLKIT_BUILD_TESTS" OFF)
+if(IPC_TOOLKIT_TEST_CCD_BENCHMARK)
+  set(IPC_TOOLKIT_CCD_BENCHMARK_DIR "" CACHE PATH "Path to the CCD benchmark directory")
+  set(IPC_TOOLKIT_CCD_NEW_BENCHMARK_DIR "" CACHE PATH "Path to the new CCD benchmark directory")
+endif()
+
 # Set default minimum C++ standard
 if(IPC_TOOLKIT_TOPLEVEL_PROJECT)
   set(CMAKE_CXX_STANDARD 17)
@@ -150,9 +157,10 @@ endif()
 include(spdlog)
 target_link_libraries(ipc_toolkit PUBLIC spdlog::spdlog)
 
+# rational-cpp (requires GMP)
 if(IPC_TOOLKIT_WITH_RATIONAL_INTERSECTION)
-  include(gmp)
-  target_link_libraries(ipc_toolkit PUBLIC GMP::GMP)
+  include(rational_cpp)
+  target_link_libraries(ipc_toolkit PUBLIC rational::rational)
 endif()
 
 # Sweep and Tiniest Queue and CCD
diff --git a/README.md b/README.md
index c5cde8305..b641b4df3 100644
--- a/README.md
+++ b/README.md
@@ -72,9 +72,9 @@ The following libraries are used in this project:
 * [Abseil](https://abseil.io/): hashing utilities
     * Enable by using the CMake option `IPC_TOOLKIT_WITH_ABSEIL`
     * Enabled by default
-* [GMP](https://gmplib.org/): rational arithmetic used for exact intersection checks
+* [rational-cpp](https://github.io/zfergus/rational-cpp): rational arithmetic used for exact intersection checks
     * Enable by using the CMake option `IPC_TOOLKIT_WITH_RATIONAL_INTERSECTION`
-    * GMP must be installed at a system level
+    * Requires [GMP](https://gmplib.org/) to be installed at a system level
 * [Etienne Vouga's Collision Detection Library](https://github.com/evouga/collisiondetection): inexact CCD
     * Included for comparison with the original IPC library
     * Enable by disabling the CMake option `IPC_TOOLKIT_WITH_CORRECT_CCD`
@@ -145,4 +145,4 @@ Additionally, you can cite the original IPC paper:
 
 ## License
 
-MIT License © 2020, the IPC-Sim organization (See <a href="https://github.com/ipc-sim/ipc-toolkit/blob/main/LICENSE"><code>LICENSE</code></a> for details)
+MIT License © 2020, the IPC-Sim organization (See <a href="https://github.com/ipc-sim/ipc-toolkit/blob/main/LICENSE"><code>LICENSE</code></a> for details).
diff --git a/cmake/find/FindGMP.cmake b/cmake/find/FindGMP.cmake
deleted file mode 100644
index b928326ad..000000000
--- a/cmake/find/FindGMP.cmake
+++ /dev/null
@@ -1,20 +0,0 @@
-# Try to find the GMP librairies
-# GMP_FOUND - system has GMP lib
-# GMP_INCLUDE_DIRS - the GMP include directory
-# GMP_LIBRARIES - Libraries needed to use GMP
-
-if (GMP_INCLUDE_DIRS AND GMP_LIBRARIES)
-  # Already in cache, be silent
-  set(GMP_FIND_QUIETLY TRUE)
-endif (GMP_INCLUDE_DIRS AND GMP_LIBRARIES)
-
-find_path(GMP_INCLUDE_DIRS NAMES gmp.h PATHS $ENV{GMP_INC} ${GMP_WINDOWS_PATH})
-find_library(GMP_LIBRARIES NAMES gmp libgmp PATHS $ENV{GMP_LIB} ${GMP_WINDOWS_PATH})
-find_library(GMPXX_LIBRARIES NAMES gmpxx libgmpxx PATHS $ENV{GMP_LIB} ${GMP_WINDOWS_PATH})
-#MESSAGE(STATUS "GMP libs: " ${GMP_LIBRARIES} " " ${GMPXX_LIBRARIES} )
-
-include(FindPackageHandleStandardArgs)
-FIND_PACKAGE_HANDLE_STANDARD_ARGS(GMP DEFAULT_MSG GMP_INCLUDE_DIRS GMP_LIBRARIES)
-
-mark_as_advanced(GMP_INCLUDE_DIRS GMP_LIBRARIES)
-message(STATUS "GMP libs: " ${GMP_LIBRARIES} " " ${GMP_INCLUDE_DIRS})
diff --git a/cmake/recipes/abseil.cmake b/cmake/recipes/abseil.cmake
index 5dabe8d5d..10113ed79 100644
--- a/cmake/recipes/abseil.cmake
+++ b/cmake/recipes/abseil.cmake
@@ -1,3 +1,5 @@
+# Abseil (https://github.com/abseil/abseil-cpp)
+# License: Apache 2.0
 if(TARGET absl::flat_hash_map)
     return()
 endif()
diff --git a/cmake/recipes/catch2.cmake b/cmake/recipes/catch2.cmake
index 1cfe6caf3..4efb1cec2 100644
--- a/cmake/recipes/catch2.cmake
+++ b/cmake/recipes/catch2.cmake
@@ -1,14 +1,5 @@
-#
-# Copyright 2020 Adobe. All rights reserved.
-# This file is licensed to you under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License. You may obtain a copy
-# of the License at http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software distributed under
-# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS
-# OF ANY KIND, either express or implied. See the License for the specific language
-# governing permissions and limitations under the License.
-#
+# Catch2 (https://github.com/catchorg/Catch2)
+# License: BSL-1.0
 if(TARGET Catch2::Catch2)
     return()
 endif()
diff --git a/cmake/recipes/ccd_query_io.cmake b/cmake/recipes/ccd_query_io.cmake
new file mode 100644
index 000000000..18fa410ff
--- /dev/null
+++ b/cmake/recipes/ccd_query_io.cmake
@@ -0,0 +1,13 @@
+# CCD Query IO (https://github.com/Continuous-Collision-Detection/CCD-Query-IO)
+# License: MIT
+if(TARGET ccd_io::ccd_io)
+    return()
+endif()
+
+message(STATUS "Third-party: creating target 'ccd_io::ccd_io'")
+
+set(CCD_IO_DOWNLOAD_SAMPLE_QUERIES ON CACHE BOOL "Download sample CCD queries" FORCE)
+set(CCD_IO_SAMPLE_QUERIES_DIR "${PROJECT_SOURCE_DIR}/tests/data/ccd-queries/" CACHE PATH "Where should we download sample queries?")
+
+include(CPM)
+CPMAddPackage("gh:Continuous-Collision-Detection/CCD-Query-IO#36f6093af81a65acc27d9f05ad32d6b5729e8d15")
\ No newline at end of file
diff --git a/cmake/recipes/eigen.cmake b/cmake/recipes/eigen.cmake
index 0bcd2abc6..30992dee2 100644
--- a/cmake/recipes/eigen.cmake
+++ b/cmake/recipes/eigen.cmake
@@ -1,14 +1,5 @@
-#
-# Copyright 2019 Adobe. All rights reserved.
-# This file is licensed to you under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License. You may obtain a copy
-# of the License at http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software distributed under
-# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS
-# OF ANY KIND, either express or implied. See the License for the specific language
-# governing permissions and limitations under the License.
-#
+# Eigen (https://gitlab.com/libeigen/eigen)
+# License: MPL 2.0
 if(TARGET Eigen3::Eigen)
     return()
 endif()
diff --git a/cmake/recipes/evouga_ccd.cmake b/cmake/recipes/evouga_ccd.cmake
index 91b857139..69a1c760e 100644
--- a/cmake/recipes/evouga_ccd.cmake
+++ b/cmake/recipes/evouga_ccd.cmake
@@ -1,4 +1,5 @@
-# Etienne Vouga's CCD Library
+# Etienne Vouga's CCD Library (https://github.com/evouga/collisiondetection.git)
+# License: ???
 if(TARGET evouga::ccd)
     return()
 endif()
@@ -15,11 +16,11 @@ CPMAddPackage(
 
 # file(GLOB EVOUGA_CCD_SOURCE_FILES "${evccd_SOURCE_DIR}/src/*.cpp")
 add_library(evouga_ccd
-    "${evccd_SOURCE_DIR}/src/CTCD.cpp"
+    "${collisiondetection_SOURCE_DIR}/src/CTCD.cpp"
 )
 add_library(evouga::ccd ALIAS evouga_ccd)
 
-target_include_directories(evouga_ccd PUBLIC "${evccd_SOURCE_DIR}/include")
+target_include_directories(evouga_ccd PUBLIC "${collisiondetection_SOURCE_DIR}/include")
 
 include(eigen)
 target_link_libraries(evouga_ccd PUBLIC Eigen3::Eigen)
diff --git a/cmake/recipes/finite_diff.cmake b/cmake/recipes/finite_diff.cmake
index 5211de439..be4ba5455 100644
--- a/cmake/recipes/finite_diff.cmake
+++ b/cmake/recipes/finite_diff.cmake
@@ -1,3 +1,5 @@
+# finite-diff (https://github.com/zfergus/finite-diff)
+# License: MIT
 if(TARGET finitediff::finitediff)
     return()
 endif()
diff --git a/cmake/recipes/gmp.cmake b/cmake/recipes/gmp.cmake
deleted file mode 100644
index 6c9068280..000000000
--- a/cmake/recipes/gmp.cmake
+++ /dev/null
@@ -1,18 +0,0 @@
-if(TARGET GMP::GMP)
-    return()
-endif()
-
-message(STATUS "Third-party: creating targets 'GMP::GMP'")
-
-# We do not have a build recipe for this, so find it as a system installed library.
-find_package(GMP REQUIRED)
-
-if(NOT ${GMP_FOUND})
-  MESSAGE(FATAL_ERROR "Unable to find GMP")
-endif()
-
-add_library(GMP_GMP INTERFACE)
-add_library(GMP::GMP ALIAS GMP_GMP)
-
-target_include_directories(GMP_GMP INTERFACE ${GMP_INCLUDE_DIRS})
-target_link_libraries(GMP_GMP INTERFACE ${GMP_LIBRARIES})
\ No newline at end of file
diff --git a/cmake/recipes/gpu_ccd.cmake b/cmake/recipes/gpu_ccd.cmake
index 81259d641..b38d20d74 100644
--- a/cmake/recipes/gpu_ccd.cmake
+++ b/cmake/recipes/gpu_ccd.cmake
@@ -1,3 +1,5 @@
+# CCD-GPU (https://github.com/dbelgrod/CCD-GPU)
+# License: MIT
 if(TARGET gpu_ccd::gpu_ccd)
     return()
 endif()
diff --git a/cmake/recipes/json.cmake b/cmake/recipes/json.cmake
index 56cc1e45f..6ad2c7e7e 100644
--- a/cmake/recipes/json.cmake
+++ b/cmake/recipes/json.cmake
@@ -1,14 +1,5 @@
-#
-# Copyright 2020 Adobe. All rights reserved.
-# This file is licensed to you under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License. You may obtain a copy
-# of the License at http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software distributed under
-# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS
-# OF ANY KIND, either express or implied. See the License for the specific language
-# governing permissions and limitations under the License.
-#
+# Nlohmann's JSON (https://github.com/nlohmann/json)
+# License: MIT
 if(TARGET nlohmann_json::nlohmann_json)
     return()
 endif()
diff --git a/cmake/recipes/libigl.cmake b/cmake/recipes/libigl.cmake
index 93b4b5789..d1b85259e 100644
--- a/cmake/recipes/libigl.cmake
+++ b/cmake/recipes/libigl.cmake
@@ -1,3 +1,5 @@
+# libigl (https://github.com/libigl/libigl)
+# License: MPL-2.0
 if(TARGET igl::core)
     return()
 endif()
@@ -9,4 +11,4 @@ set(LIBIGL_PREDICATES ON CACHE BOOL "Use exact predicates" FORCE)
 include(eigen)
 
 include(CPM)
-CPMAddPackage("gh:libigl/libigl@2.4.0")
+CPMAddPackage("gh:libigl/libigl@2.5.0")
\ No newline at end of file
diff --git a/cmake/recipes/onetbb.cmake b/cmake/recipes/onetbb.cmake
index 979f6dc19..318d6fdc2 100644
--- a/cmake/recipes/onetbb.cmake
+++ b/cmake/recipes/onetbb.cmake
@@ -1,14 +1,5 @@
-#
-# Copyright 2021 Adobe. All rights reserved.
-# This file is licensed to you under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License. You may obtain a copy
-# of the License at http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software distributed under
-# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS
-# OF ANY KIND, either express or implied. See the License for the specific language
-# governing permissions and limitations under the License.
-#
+# oneTBB (https://github.com/oneapi-src/oneTBB)
+# License: Apache-2.0
 if(TARGET TBB::tbb)
     return()
 endif()
diff --git a/cmake/recipes/pybind11.cmake b/cmake/recipes/pybind11.cmake
index e078f6e8e..d25797a13 100644
--- a/cmake/recipes/pybind11.cmake
+++ b/cmake/recipes/pybind11.cmake
@@ -1,14 +1,5 @@
-#
-# Copyright 2020 Adobe. All rights reserved.
-# This file is licensed to you under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License. You may obtain a copy
-# of the License at http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software distributed under
-# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS
-# OF ANY KIND, either express or implied. See the License for the specific language
-# governing permissions and limitations under the License.
-#
+# pybind11 (https://github.com/pybind/pybind11)
+# License: BSD-style
 if(TARGET pybind11::pybind11)
     return()
 endif()
diff --git a/cmake/recipes/rational_cpp.cmake b/cmake/recipes/rational_cpp.cmake
new file mode 100644
index 000000000..43e02d6a1
--- /dev/null
+++ b/cmake/recipes/rational_cpp.cmake
@@ -0,0 +1,10 @@
+# rational-cpp (https://github.com/zfergus/rational-cpp)
+# License: MIT
+if(TARGET rational::rational)
+    return()
+endif()
+
+message(STATUS "Third-party: creating target 'rational::rational'")
+
+include(CPM)
+CPMAddPackage("gh:zfergus/rational-cpp#687d4ea3436ada7231b8920f3cd5b02b438c21aa")
\ No newline at end of file
diff --git a/cmake/recipes/robin_map.cmake b/cmake/recipes/robin_map.cmake
index fa83d597d..713963384 100644
--- a/cmake/recipes/robin_map.cmake
+++ b/cmake/recipes/robin_map.cmake
@@ -1,3 +1,5 @@
+# robin-map (https://github.com/Tessil/robin-map)
+# License: MIT
 if(TARGET tsl::robin_map)
     return()
 endif()
diff --git a/cmake/recipes/spdlog.cmake b/cmake/recipes/spdlog.cmake
index 90ca20dd5..beac4ce3c 100644
--- a/cmake/recipes/spdlog.cmake
+++ b/cmake/recipes/spdlog.cmake
@@ -1,14 +1,5 @@
-#
-# Copyright 2020 Adobe. All rights reserved.
-# This file is licensed to you under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License. You may obtain a copy
-# of the License at http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software distributed under
-# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS
-# OF ANY KIND, either express or implied. See the License for the specific language
-# governing permissions and limitations under the License.
-#
+# spdlog (https://github.com/gabime/spdlog)
+# License: MIT
 if(TARGET spdlog::spdlog)
     return()
 endif()
diff --git a/cmake/recipes/sweep_and_tiniest_queue.cmake b/cmake/recipes/sweep_and_tiniest_queue.cmake
index 190ea60d6..75c8361bc 100644
--- a/cmake/recipes/sweep_and_tiniest_queue.cmake
+++ b/cmake/recipes/sweep_and_tiniest_queue.cmake
@@ -1,3 +1,5 @@
+# Broadphase GPU (https://github.com/dbelgrod/broadphase-gpu)
+# License: MIT
 if(TARGET STQ::CPU)
     return()
 endif()
diff --git a/cmake/recipes/tight_inclusion.cmake b/cmake/recipes/tight_inclusion.cmake
index d01b0c81b..7d24599d0 100644
--- a/cmake/recipes/tight_inclusion.cmake
+++ b/cmake/recipes/tight_inclusion.cmake
@@ -1,3 +1,5 @@
+# Tight Inclusion (https://github.com/Continuous-Collision-Detection/Tight-Inclusion)
+# License: MIT
 if(TARGET tight_inclusion::tight_inclusion)
     return()
 endif()
diff --git a/docs/source/cpp.rst b/docs/source/cpp.rst
index ad6a56c5a..43d3ba0dd 100644
--- a/docs/source/cpp.rst
+++ b/docs/source/cpp.rst
@@ -73,10 +73,10 @@ Optional
   * Enable by using the CMake option :cmake:`IPC_TOOLKIT_WITH_ABSEIL`
   * Enabled by default
 
-* `GMP <https://gmplib.org/>`__: rational arithmetic used for exact intersection checks
+* `rational-cpp <https://github.io/zfergus/rational-cpp>`__: rational arithmetic used for exact intersection checks
 
-  * Enable by using the CMake option :cmake:`IPC_TOOLKIT_WITH_RATIONAL_INTERSECTION`
-  * GMP must be installed at a system level
+    * Enable by using the CMake option :cmake:`IPC_TOOLKIT_WITH_RATIONAL_INTERSECTION`
+    * Requires `GMP <https://gmplib.org/>`__ to be installed at a system level
 
 * `Etienne Vouga's Collision Detection Library <https://github.com/evouga/collisiondetection>`__: inexact CCD
 
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index 2f3565dd0..4240f3fac 100755
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -27,6 +27,7 @@ pybind11_add_module(ipctk
 
   src/ccd/aabb.cpp
   src/ccd/ccd.cpp
+  src/ccd/additive_ccd.cpp
   src/ccd/inexact_point_edge.cpp
   src/ccd/point_static_plane.cpp
 
diff --git a/python/src/bindings.cpp b/python/src/bindings.cpp
index 74896f211..910e3fb67 100644
--- a/python/src/bindings.cpp
+++ b/python/src/bindings.cpp
@@ -36,6 +36,7 @@ PYBIND11_MODULE(ipctk, m)
     // ccd
     define_ccd_aabb(m);
     define_ccd(m);
+    define_additive_ccd(m);
     define_inexact_point_edge(m);
     define_point_static_plane(m);
 
diff --git a/python/src/ccd/additive_ccd.cpp b/python/src/ccd/additive_ccd.cpp
new file mode 100644
index 000000000..885c85227
--- /dev/null
+++ b/python/src/ccd/additive_ccd.cpp
@@ -0,0 +1,169 @@
+#include <common.hpp>
+
+#include <ipc/ccd/additive_ccd.hpp>
+
+namespace py = pybind11;
+using namespace ipc;
+
+void define_additive_ccd(py::module_& m)
+{
+    using namespace ipc::additive_ccd;
+
+    auto m_accd = m.def_submodule(
+        "additive_ccd", "Additive CCD method of [Li et al. 2021].");
+
+    m_accd.def(
+        "point_point_ccd",
+        [](const VectorMax3d& p0_t0, const VectorMax3d& p1_t0,
+           const VectorMax3d& p0_t1, const VectorMax3d& p1_t1,
+           const double min_distance = 0.0, const double tmax = 1.0,
+           const double conservative_rescaling =
+               DEFAULT_CCD_CONSERVATIVE_RESCALING) {
+            double toi;
+            bool r = point_point_ccd(
+                p0_t0, p1_t0, p0_t1, p1_t1, toi, min_distance, tmax,
+                conservative_rescaling);
+            return std::make_tuple(r, toi);
+        },
+        R"ipc_Qu8mg5v7(
+        Computes the time of impact between two points using continuous collision detection.
+
+        Parameters:
+            p0_t0: The initial position of the first point.
+            p1_t0: The initial position of the second point.
+            p0_t1: The final position of the first point.
+            p1_t1: The final position of the second point.
+            min_distance: The minimum distance between two objects.
+            tmax: The maximum time to check for collisions.
+            conservative_rescaling: The conservative rescaling of the time of impact.
+
+        Returns:
+            Tuple of:
+            True if a collision was detected, false otherwise.
+            The time of impact between the two points.
+        )ipc_Qu8mg5v7",
+        py::arg("p0_t0"), py::arg("p1_t0"), py::arg("p0_t1"), py::arg("p1_t1"),
+        py::arg("min_distance") = 0.0, py::arg("tmax") = 1.0,
+        py::arg("conservative_rescaling") = DEFAULT_CCD_CONSERVATIVE_RESCALING);
+
+    m_accd.def(
+        "point_edge_ccd",
+        [](const VectorMax3d& p_t0, const VectorMax3d& e0_t0,
+           const VectorMax3d& e1_t0, const VectorMax3d& p_t1,
+           const VectorMax3d& e0_t1, const VectorMax3d& e1_t1,
+           const double min_distance = 0.0, const double tmax = 1.0,
+           const double conservative_rescaling =
+               DEFAULT_CCD_CONSERVATIVE_RESCALING) {
+            double toi;
+            bool r = point_edge_ccd(
+                p_t0, e0_t0, e1_t0, p_t1, e0_t1, e1_t1, toi, min_distance, tmax,
+                conservative_rescaling);
+            return std::make_tuple(r, toi);
+        },
+        R"ipc_Qu8mg5v7(
+        Computes the time of impact between a point and an edge using continuous collision detection.
+
+        Parameters:
+            p_t0: The initial position of the point.
+            e0_t0: The initial position of the first endpoint of the edge.
+            e1_t0: The initial position of the second endpoint of the edge.
+            p_t1: The final position of the point.
+            e0_t1: The final position of the first endpoint of the edge.
+            e1_t1: The final position of the second endpoint of the edge.
+            min_distance: The minimum distance between two objects.
+            tmax: The maximum time to check for collisions.
+            conservative_rescaling: The conservative rescaling of the time of impact.
+
+        Returns:
+            Tuple of:
+            True if a collision was detected, false otherwise.
+            The time of impact between the point and the edge.
+        )ipc_Qu8mg5v7",
+        py::arg("p_t0"), py::arg("e0_t0"), py::arg("e1_t0"), py::arg("p_t1"),
+        py::arg("e0_t1"), py::arg("e1_t1"), py::arg("min_distance") = 0.0,
+        py::arg("tmax") = 1.0,
+        py::arg("conservative_rescaling") = DEFAULT_CCD_CONSERVATIVE_RESCALING);
+
+    m_accd.def(
+        "point_triangle_ccd",
+        [](const Eigen::Vector3d& p_t0, const Eigen::Vector3d& t0_t0,
+           const Eigen::Vector3d& t1_t0, const Eigen::Vector3d& t2_t0,
+           const Eigen::Vector3d& p_t1, const Eigen::Vector3d& t0_t1,
+           const Eigen::Vector3d& t1_t1, const Eigen::Vector3d& t2_t1,
+           const double min_distance = 0.0, const double tmax = 1.0,
+           const double conservative_rescaling =
+               DEFAULT_CCD_CONSERVATIVE_RESCALING) {
+            double toi;
+            bool r = point_triangle_ccd(
+                p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1, toi,
+                min_distance, tmax, conservative_rescaling);
+            return std::make_tuple(r, toi);
+        },
+        R"ipc_Qu8mg5v7(
+        Computes the time of impact between a point and a triangle using continuous collision detection.
+
+        Parameters:
+            p_t0: The initial position of the point.
+            t0_t0: The initial position of the first vertex of the triangle.
+            t1_t0: The initial position of the second vertex of the triangle.
+            t2_t0: The initial position of the third vertex of the triangle.
+            p_t1: The final position of the point.
+            t0_t1: The final position of the first vertex of the triangle.
+            t1_t1: The final position of the second vertex of the triangle.
+            t2_t1: The final position of the third vertex of the triangle.
+            min_distance: The minimum distance between two objects.
+            tmax: The maximum time to check for collisions.
+            conservative_rescaling: The conservative rescaling of the time of impact.
+
+        Returns:
+            Tuple of:
+            True if a collision was detected, false otherwise.
+            The time of impact between the point and the triangle.
+        )ipc_Qu8mg5v7",
+        py::arg("p_t0"), py::arg("t0_t0"), py::arg("t1_t0"), py::arg("t2_t0"),
+        py::arg("p_t1"), py::arg("t0_t1"), py::arg("t1_t1"), py::arg("t2_t1"),
+        py::arg("min_distance") = 0.0, py::arg("tmax") = 1.0,
+        py::arg("conservative_rescaling") = DEFAULT_CCD_CONSERVATIVE_RESCALING);
+
+    m_accd.def(
+        "edge_edge_ccd",
+        [](const Eigen::Vector3d& ea0_t0, const Eigen::Vector3d& ea1_t0,
+           const Eigen::Vector3d& eb0_t0, const Eigen::Vector3d& eb1_t0,
+           const Eigen::Vector3d& ea0_t1, const Eigen::Vector3d& ea1_t1,
+           const Eigen::Vector3d& eb0_t1, const Eigen::Vector3d& eb1_t1,
+           const double min_distance = 0.0, const double tmax = 1.0,
+           const double conservative_rescaling =
+               DEFAULT_CCD_CONSERVATIVE_RESCALING) {
+            double toi;
+            bool r = edge_edge_ccd(
+                ea0_t0, ea1_t0, eb0_t0, eb1_t0, ea0_t1, ea1_t1, eb0_t1, eb1_t1,
+                toi, min_distance, tmax, conservative_rescaling);
+            return std::make_tuple(r, toi);
+        },
+        R"ipc_Qu8mg5v7(
+        Computes the time of impact between two edges using continuous collision detection.
+
+        Parameters:
+            ea0_t0: The initial position of the first endpoint of the first edge.
+            ea1_t0: The initial position of the second endpoint of the first edge.
+            eb0_t0: The initial position of the first endpoint of the second edge.
+            eb1_t0: The initial position of the second endpoint of the second edge.
+            ea0_t1: The final position of the first endpoint of the first edge.
+            ea1_t1: The final position of the second endpoint of the first edge.
+            eb0_t1: The final position of the first endpoint of the second edge.
+            eb1_t1: The final position of the second endpoint of the second edge.
+            min_distance: The minimum distance between two objects.
+            tmax: The maximum time to check for collisions.
+            conservative_rescaling: The conservative rescaling of the time of impact.
+
+        Returns:
+            Tuple of:
+            True if a collision was detected, false otherwise.
+            The time of impact between the two edges.
+        )ipc_Qu8mg5v7",
+        py::arg("ea0_t0"), py::arg("ea1_t0"), py::arg("eb0_t0"),
+        py::arg("eb1_t0"), py::arg("ea0_t1"), py::arg("ea1_t1"),
+        py::arg("eb0_t1"), py::arg("eb1_t1"), py::arg("min_distance") = 0.0,
+        py::arg("tmax") = 1.0,
+        py::arg("conservative_rescaling") = DEFAULT_CCD_CONSERVATIVE_RESCALING);
+}
\ No newline at end of file
diff --git a/python/src/ccd/bindings.hpp b/python/src/ccd/bindings.hpp
index 17445284c..fb2398c57 100644
--- a/python/src/ccd/bindings.hpp
+++ b/python/src/ccd/bindings.hpp
@@ -5,5 +5,6 @@ namespace py = pybind11;
 
 void define_ccd_aabb(py::module_& m);
 void define_ccd(py::module_& m);
+void define_additive_ccd(py::module_& m);
 void define_inexact_point_edge(py::module_& m);
 void define_point_static_plane(py::module_& m);
\ No newline at end of file
diff --git a/python/src/ccd/ccd.cpp b/python/src/ccd/ccd.cpp
index 61a98a419..95b6163b0 100644
--- a/python/src/ccd/ccd.cpp
+++ b/python/src/ccd/ccd.cpp
@@ -23,10 +23,30 @@ void define_ccd(py::module_& m)
                 tolerance, max_iterations, conservative_rescaling);
             return std::make_tuple(r, toi);
         },
-        "", py::arg("p_t0"), py::arg("e0_t0"), py::arg("e1_t0"),
-        py::arg("p_t1"), py::arg("e0_t1"), py::arg("e1_t1"),
-        py::arg("min_distance") = 0.0, py::arg("tmax") = 1.0,
-        py::arg("tolerance") = DEFAULT_CCD_TOLERANCE,
+        R"ipc_Qu8mg5v7(
+        Computes the time of impact between a point and an edge in 2D using continuous collision detection.
+
+        Parameters:
+            p_t0: The initial position of the point.
+            e0_t0: The initial position of the first endpoint of the edge.
+            e1_t0: The initial position of the second endpoint of the edge.
+            p_t1: The final position of the point.
+            e0_t1: The final position of the first endpoint of the edge.
+            e1_t1: The final position of the second endpoint of the edge.
+            min_distance: The minimum distance between the objects.
+            tmax: The maximum time to check for collisions.
+            tolerance: The error tolerance for the time of impact.
+            max_iterations: The maximum number of iterations to perform.
+            conservative_rescaling: The conservative rescaling of the time of impact.
+
+        Returns:
+            Tuple of:
+            True if a collision was detected, false otherwise.
+            The time of impact between the point and the edge.
+        )ipc_Qu8mg5v7",
+        py::arg("p_t0"), py::arg("e0_t0"), py::arg("e1_t0"), py::arg("p_t1"),
+        py::arg("e0_t1"), py::arg("e1_t1"), py::arg("min_distance") = 0.0,
+        py::arg("tmax") = 1.0, py::arg("tolerance") = DEFAULT_CCD_TOLERANCE,
         py::arg("max_iterations") = DEFAULT_CCD_MAX_ITERATIONS,
         py::arg("conservative_rescaling") = DEFAULT_CCD_CONSERVATIVE_RESCALING);
 
@@ -45,8 +65,27 @@ void define_ccd(py::module_& m)
                 max_iterations, conservative_rescaling);
             return std::make_tuple(r, toi);
         },
-        "", py::arg("p0_t0"), py::arg("p1_t0"), py::arg("p0_t1"),
-        py::arg("p1_t1"), py::arg("min_distance") = 0.0, py::arg("tmax") = 1.0,
+        R"ipc_Qu8mg5v7(
+        Computes the time of impact between two points in 3D using continuous collision detection.
+
+        Parameters:
+            p0_t0: The initial position of the first point.
+            p1_t0: The initial position of the second point.
+            p0_t1: The final position of the first point.
+            p1_t1: The final position of the second point.
+            min_distance: The minimum distance between the objects.
+            tmax: The maximum time to check for collisions.
+            tolerance: The error tolerance for the time of impact.
+            max_iterations: The maximum number of iterations to perform.
+            conservative_rescaling: The conservative rescaling of the time of impact.
+
+        Returns:
+            Tuple of:
+            True if a collision was detected, false otherwise.
+            The time of impact between the two points.
+        )ipc_Qu8mg5v7",
+        py::arg("p0_t0"), py::arg("p1_t0"), py::arg("p0_t1"), py::arg("p1_t1"),
+        py::arg("min_distance") = 0.0, py::arg("tmax") = 1.0,
         py::arg("tolerance") = DEFAULT_CCD_TOLERANCE,
         py::arg("max_iterations") = DEFAULT_CCD_MAX_ITERATIONS,
         py::arg("conservative_rescaling") = DEFAULT_CCD_CONSERVATIVE_RESCALING);
@@ -67,10 +106,30 @@ void define_ccd(py::module_& m)
                 tolerance, max_iterations, conservative_rescaling);
             return std::make_tuple(r, toi);
         },
-        "", py::arg("p_t0"), py::arg("e0_t0"), py::arg("e1_t0"),
-        py::arg("p_t1"), py::arg("e0_t1"), py::arg("e1_t1"),
-        py::arg("min_distance") = 0.0, py::arg("tmax") = 1.0,
-        py::arg("tolerance") = DEFAULT_CCD_TOLERANCE,
+        R"ipc_Qu8mg5v7(
+        Computes the time of impact between a point and an edge in 3D using continuous collision detection.
+
+        Parameters:
+            p_t0: The initial position of the point.
+            e0_t0: The initial position of the first endpoint of the edge.
+            e1_t0: The initial position of the second endpoint of the edge.
+            p_t1: The final position of the point.
+            e0_t1: The final position of the first endpoint of the edge.
+            e1_t1: The final position of the second endpoint of the edge.
+            min_distance: The minimum distance between the objects.
+            tmax: The maximum time to check for collisions.
+            tolerance: The error tolerance for the time of impact.
+            max_iterations: The maximum number of iterations to perform.
+            conservative_rescaling: The conservative rescaling of the time of impact.
+
+        Returns:
+            Tuple of:
+            True if a collision was detected, false otherwise.
+            The time of impact between the point and the edge.
+        )ipc_Qu8mg5v7",
+        py::arg("p_t0"), py::arg("e0_t0"), py::arg("e1_t0"), py::arg("p_t1"),
+        py::arg("e0_t1"), py::arg("e1_t1"), py::arg("min_distance") = 0.0,
+        py::arg("tmax") = 1.0, py::arg("tolerance") = DEFAULT_CCD_TOLERANCE,
         py::arg("max_iterations") = DEFAULT_CCD_MAX_ITERATIONS,
         py::arg("conservative_rescaling") = DEFAULT_CCD_CONSERVATIVE_RESCALING);
 
@@ -92,9 +151,32 @@ void define_ccd(py::module_& m)
                 conservative_rescaling);
             return std::make_tuple(r, toi);
         },
-        "", py::arg("p_t0"), py::arg("t0_t0"), py::arg("t1_t0"),
-        py::arg("t2_t0"), py::arg("p_t1"), py::arg("t0_t1"), py::arg("t1_t1"),
-        py::arg("t2_t1"), py::arg("min_distance") = 0.0, py::arg("tmax") = 1.0,
+        R"ipc_Qu8mg5v7(
+        Computes the time of impact between a point and a triangle in 3D using continuous collision detection.
+
+        Parameters:
+            p_t0: The initial position of the point.
+            t0_t0: The initial position of the first vertex of the triangle.
+            t1_t0: The initial position of the second vertex of the triangle.
+            t2_t0: The initial position of the third vertex of the triangle.
+            p_t1: The final position of the point.
+            t0_t1: The final position of the first vertex of the triangle.
+            t1_t1: The final position of the second vertex of the triangle.
+            t2_t1: The final position of the third vertex of the triangle.
+            min_distance: The minimum distance between the objects.
+            tmax: The maximum time to check for collisions.
+            tolerance: The error tolerance for the time of impact.
+            max_iterations: The maximum number of iterations to perform.
+            conservative_rescaling: The conservative rescaling of the time of impact.
+
+        Returns:
+            Tuple of:
+            True if a collision was detected, false otherwise.
+            The time of impact between the point and the triangle.
+        )ipc_Qu8mg5v7",
+        py::arg("p_t0"), py::arg("t0_t0"), py::arg("t1_t0"), py::arg("t2_t0"),
+        py::arg("p_t1"), py::arg("t0_t1"), py::arg("t1_t1"), py::arg("t2_t1"),
+        py::arg("min_distance") = 0.0, py::arg("tmax") = 1.0,
         py::arg("tolerance") = DEFAULT_CCD_TOLERANCE,
         py::arg("max_iterations") = DEFAULT_CCD_MAX_ITERATIONS,
         py::arg("conservative_rescaling") = DEFAULT_CCD_CONSERVATIVE_RESCALING);
@@ -117,7 +199,30 @@ void define_ccd(py::module_& m)
                 conservative_rescaling);
             return std::make_tuple(r, toi);
         },
-        "", py::arg("ea0_t0"), py::arg("ea1_t0"), py::arg("eb0_t0"),
+        R"ipc_Qu8mg5v7(
+        Computes the time of impact between two edges in 3D using continuous collision detection.
+
+        Parameters:
+            ea0_t0: The initial position of the first endpoint of the first edge.
+            ea1_t0: The initial position of the second endpoint of the first edge.
+            eb0_t0: The initial position of the first endpoint of the second edge.
+            eb1_t0: The initial position of the second endpoint of the second edge.
+            ea0_t1: The final position of the first endpoint of the first edge.
+            ea1_t1: The final position of the second endpoint of the first edge.
+            eb0_t1: The final position of the first endpoint of the second edge.
+            eb1_t1: The final position of the second endpoint of the second edge.
+            min_distance: The minimum distance between the objects.
+            tmax: The maximum time to check for collisions.
+            tolerance: The error tolerance for the time of impact.
+            max_iterations: The maximum number of iterations to perform.
+            conservative_rescaling: The conservative rescaling of the time of impact.
+
+        Returns:
+            Tuple of:
+            True if a collision was detected, false otherwise.
+            The time of impact between the two edges.
+        )ipc_Qu8mg5v7",
+        py::arg("ea0_t0"), py::arg("ea1_t0"), py::arg("eb0_t0"),
         py::arg("eb1_t0"), py::arg("ea0_t1"), py::arg("ea1_t1"),
         py::arg("eb0_t1"), py::arg("eb1_t1"), py::arg("min_distance") = 0.0,
         py::arg("tmax") = 1.0, py::arg("tolerance") = DEFAULT_CCD_TOLERANCE,
@@ -140,10 +245,30 @@ void define_ccd(py::module_& m)
                 tolerance, max_iterations, conservative_rescaling);
             return std::make_tuple(r, toi);
         },
-        "", py::arg("p_t0"), py::arg("e0_t0"), py::arg("e1_t0"),
-        py::arg("p_t1"), py::arg("e0_t1"), py::arg("e1_t1"),
-        py::arg("min_distance") = 0.0, py::arg("tmax") = 1.0,
-        py::arg("tolerance") = DEFAULT_CCD_TOLERANCE,
+        R"ipc_Qu8mg5v7(
+        Computes the time of impact between a point and an edge in 2D or 3D using continuous collision detection.
+
+        Parameters:
+            p_t0: The initial position of the point.
+            e0_t0: The initial position of the first endpoint of the edge.
+            e1_t0: The initial position of the second endpoint of the edge.
+            p_t1: The final position of the point.
+            e0_t1: The final position of the first endpoint of the edge.
+            e1_t1: The final position of the second endpoint of the edge.
+            min_distance: The minimum distance between the objects.
+            tmax: The maximum time to check for collisions.
+            tolerance: The error tolerance for the time of impact.
+            max_iterations: The maximum number of iterations to perform.
+            conservative_rescaling: The conservative rescaling of the time of impact.
+
+        Returns:
+            Tuple of:
+            True if a collision was detected, false otherwise.
+            The time of impact between the point and the edge.
+        )ipc_Qu8mg5v7",
+        py::arg("p_t0"), py::arg("e0_t0"), py::arg("e1_t0"), py::arg("p_t1"),
+        py::arg("e0_t1"), py::arg("e1_t1"), py::arg("min_distance") = 0.0,
+        py::arg("tmax") = 1.0, py::arg("tolerance") = DEFAULT_CCD_TOLERANCE,
         py::arg("max_iterations") = DEFAULT_CCD_MAX_ITERATIONS,
         py::arg("conservative_rescaling") = DEFAULT_CCD_CONSERVATIVE_RESCALING);
 }
diff --git a/src/ipc/ccd/CMakeLists.txt b/src/ipc/ccd/CMakeLists.txt
index ab49a7630..3a6f0d9c5 100644
--- a/src/ipc/ccd/CMakeLists.txt
+++ b/src/ipc/ccd/CMakeLists.txt
@@ -1,6 +1,8 @@
 set(SOURCES
   aabb.cpp
   aabb.hpp
+  additive_ccd.cpp
+  additive_ccd.hpp
   ccd.cpp
   ccd.hpp
   inexact_point_edge.cpp
diff --git a/src/ipc/ccd/additive_ccd.cpp b/src/ipc/ccd/additive_ccd.cpp
new file mode 100644
index 000000000..39af9dde9
--- /dev/null
+++ b/src/ipc/ccd/additive_ccd.cpp
@@ -0,0 +1,316 @@
+//
+// Source modified from https://github.com/ipc-sim/Codim-IPC
+// under Appache-2.0 License.
+//
+// Modifications:
+//  • remove broad phase functions
+//  • refactor code to use a single implementation of the additive_ccd algorithm
+//  • utilize our distance function rather than porting the Codim-IPC versions
+//  • return true if the initial distance is less than the minimum distance
+//  • add an explicit tmax parameter rather than relying on the initial value of
+//    toi
+//
+// NOTE: These methods are provided for reference comparison with [Li et al.
+// 2021] and is not utilized by the high-level functionality. In compairson to
+// Tight Inclusion CCD, this CCD method is not provably conservative and so can
+// potentially produce false negatives (i.e., miss collisions) due to
+// floating-point rounding error. However, it is much faster than Tight
+// Inclusion CCD (>100×) and very robust due to the gaps and conservative
+// rescaling used.
+//
+
+#include "additive_ccd.hpp"
+
+#include <ipc/distance/point_point.hpp>
+#include <ipc/distance/point_edge.hpp>
+#include <ipc/distance/point_triangle.hpp>
+#include <ipc/distance/edge_edge.hpp>
+
+namespace ipc::additive_ccd {
+
+namespace {
+    template <typename... Args> void subtract_mean(Args&... args)
+    {
+        constexpr double n = sizeof...(args);
+        static_assert(n > 0, "At least one argument is required");
+
+        using T = typename std::tuple_element<0, std::tuple<Args...>>::type;
+        const int dim = std::get<0>(std::tuple<const Args&...>(args...)).size();
+
+        T mean = T::Zero(dim);
+        for (const T& value : { args... }) {
+            mean += value;
+        }
+        mean /= n;
+
+        for (T* value : { &args... }) {
+            (*value) -= mean;
+        }
+    }
+
+    VectorMax12d stack(const VectorMax3d& x) { return x; }
+
+    template <typename... Args>
+    VectorMax12d stack(const VectorMax3d& x0, const Args&... args)
+    {
+        VectorMax12d x(x0.size() * (1 + sizeof...(args)));
+        x.head(x0.size()) = x0;
+        x.tail(x0.size() * sizeof...(args)) = stack(args...);
+        return x;
+    }
+} // namespace
+
+bool additive_ccd(
+    VectorMax12d x,
+    const VectorMax12d& dx,
+    const std::function<double(const VectorMax12d&)>& distance_squared,
+    const double max_disp_mag,
+    double& toi,
+    const double min_distance,
+    const double tmax,
+    const double conservative_rescaling)
+{
+    assert(conservative_rescaling > 0 && conservative_rescaling < 1);
+
+    const double min_distance_sq = min_distance * min_distance;
+
+    double d, d_sq;
+    d = std::sqrt(d_sq = distance_squared(x));
+    assert(d > min_distance);
+
+    double d_func = d_sq - min_distance_sq;
+    assert(d_func > 0);
+    const double gap = // (d - ξ) = (d² - ξ²) / (d + ξ)
+        (1 - conservative_rescaling) * d_func / (d + min_distance);
+
+    toi = 0;
+    while (true) {
+        // tₗ = η ⋅ (d - ξ) / lₚ = η ⋅ (d² - ξ²) / (lₚ ⋅ (d + ξ))
+        const double toi_lower_bound = conservative_rescaling * d_func
+            / ((d + min_distance) * max_disp_mag);
+
+        x += toi_lower_bound * dx;
+
+        d = std::sqrt(d_sq = distance_squared(x));
+
+        d_func = d_sq - min_distance_sq;
+        assert(d_func > 0);
+        if (toi > 0 && d_func / (d + min_distance) < gap) {
+            break; // distance (including thickness) is less than gap
+        }
+
+        toi += toi_lower_bound;
+        if (toi > tmax) {
+            return false; // collision occurs after tmax
+        }
+    }
+
+    return true;
+}
+
+bool point_point_ccd(
+    const VectorMax3d& p0_t0,
+    const VectorMax3d& p1_t0,
+    const VectorMax3d& p0_t1,
+    const VectorMax3d& p1_t1,
+    double& toi,
+    const double min_distance,
+    const double tmax,
+    const double conservative_rescaling)
+{
+    const int dim = p0_t0.size();
+    assert(dim == p1_t0.size() && dim == p0_t1.size() && dim == p1_t1.size());
+
+    const double initial_distance = point_point_distance(p0_t0, p1_t0);
+    if (initial_distance <= min_distance * min_distance) {
+        logger().warn(
+            "Initial distance {} ≤ d_min={}, returning toi=0!",
+            std::sqrt(initial_distance), min_distance);
+        toi = 0;
+        return true;
+    }
+
+    VectorMax3d dp0 = p0_t1 - p0_t0;
+    VectorMax3d dp1 = p1_t1 - p1_t0;
+    subtract_mean(dp0, dp1);
+
+    const double max_disp_mag = dp0.norm() + dp1.norm();
+    if (max_disp_mag == 0) {
+        return false;
+    }
+
+    auto distance_squared = [dim](const VectorMax12d& x) {
+        return point_point_distance(x.head(dim), x.tail(dim));
+    };
+
+    VectorMax12d x = stack(p0_t0, p1_t0);
+    const VectorMax12d dx = stack(dp0, dp1);
+
+    return additive_ccd(
+        x, dx, distance_squared, max_disp_mag, toi, min_distance, tmax,
+        conservative_rescaling);
+}
+
+bool point_edge_ccd(
+    const VectorMax3d& p_t0,
+    const VectorMax3d& e0_t0,
+    const VectorMax3d& e1_t0,
+    const VectorMax3d& p_t1,
+    const VectorMax3d& e0_t1,
+    const VectorMax3d& e1_t1,
+    double& toi,
+    const double min_distance,
+    const double tmax,
+    const double conservative_rescaling)
+{
+    const int dim = p_t0.size();
+    assert(dim == e0_t0.size() && dim == e1_t0.size());
+    assert(dim == p_t1.size() && dim == e0_t1.size() && dim == e1_t1.size());
+
+    const double initial_distance = point_edge_distance(p_t0, e0_t0, e1_t0);
+    if (initial_distance <= min_distance * min_distance) {
+        logger().warn(
+            "Initial distance {} ≤ d_min={}, returning toi=0!",
+            std::sqrt(initial_distance), min_distance);
+        toi = 0;
+        return true;
+    }
+
+    VectorMax3d dp = p_t1 - p_t0;
+    VectorMax3d de0 = e0_t1 - e0_t0;
+    VectorMax3d de1 = e1_t1 - e1_t0;
+    subtract_mean(dp, de0, de1);
+
+    const double max_disp_mag =
+        dp.norm() + std::sqrt(std::max(de0.squaredNorm(), de1.squaredNorm()));
+    if (max_disp_mag == 0) {
+        return false;
+    }
+
+    VectorMax12d x = stack(p_t0, e0_t0, e1_t0);
+    const VectorMax12d dx = stack(dp, de0, de1);
+
+    auto distance_squared = [dim](const VectorMax12d& x) {
+        return point_edge_distance(
+            x.head(dim), x.segment(dim, dim), x.tail(dim));
+    };
+
+    return additive_ccd(
+        x, dx, distance_squared, max_disp_mag, toi, min_distance, tmax,
+        conservative_rescaling);
+}
+
+bool point_triangle_ccd(
+    const Eigen::Vector3d& p_t0,
+    const Eigen::Vector3d& t0_t0,
+    const Eigen::Vector3d& t1_t0,
+    const Eigen::Vector3d& t2_t0,
+    const Eigen::Vector3d& p_t1,
+    const Eigen::Vector3d& t0_t1,
+    const Eigen::Vector3d& t1_t1,
+    const Eigen::Vector3d& t2_t1,
+    double& toi,
+    const double min_distance,
+    const double tmax,
+    const double conservative_rescaling)
+{
+    const double initial_distance =
+        point_triangle_distance(p_t0, t0_t0, t1_t0, t2_t0);
+    if (initial_distance <= min_distance * min_distance) {
+        logger().warn(
+            "Initial distance {} ≤ d_min={}, returning toi=0!",
+            std::sqrt(initial_distance), min_distance);
+        toi = 0;
+        return true;
+    }
+
+    Eigen::Vector3d dp = p_t1 - p_t0;
+    Eigen::Vector3d dt0 = t0_t1 - t0_t0;
+    Eigen::Vector3d dt1 = t1_t1 - t1_t0;
+    Eigen::Vector3d dt2 = t2_t1 - t2_t0;
+    subtract_mean(dp, dt0, dt1, dt2);
+
+    VectorMax12d x = stack(p_t0, t0_t0, t1_t0, t2_t0);
+    const VectorMax12d dx = stack(dp, dt0, dt1, dt2);
+
+    const double max_disp_mag = dp.norm()
+        + std::sqrt(std::max(
+            { dt0.squaredNorm(), dt1.squaredNorm(), dt2.squaredNorm() }));
+    if (max_disp_mag == 0) {
+        return false;
+    }
+
+    auto distance_squared = [](const VectorMax12d& x) {
+        return point_triangle_distance(
+            x.head<3>(), x.segment<3>(3), x.segment<3>(6), x.tail<3>());
+    };
+
+    return additive_ccd(
+        x, dx, distance_squared, max_disp_mag, toi, min_distance, tmax,
+        conservative_rescaling);
+}
+
+bool edge_edge_ccd(
+    const Eigen::Vector3d& ea0_t0,
+    const Eigen::Vector3d& ea1_t0,
+    const Eigen::Vector3d& eb0_t0,
+    const Eigen::Vector3d& eb1_t0,
+    const Eigen::Vector3d& ea0_t1,
+    const Eigen::Vector3d& ea1_t1,
+    const Eigen::Vector3d& eb0_t1,
+    const Eigen::Vector3d& eb1_t1,
+    double& toi,
+    const double min_distance,
+    const double tmax,
+    const double conservative_rescaling)
+{
+    const double initial_distance =
+        edge_edge_distance(ea0_t0, ea1_t0, eb0_t0, eb1_t0);
+    if (initial_distance <= min_distance * min_distance) {
+        logger().warn(
+            "Initial distance {} ≤ d_min={}, returning toi=0!",
+            std::sqrt(initial_distance), min_distance);
+        toi = 0;
+        return true;
+    }
+
+    Eigen::Vector3d dea0 = ea0_t1 - ea0_t0;
+    Eigen::Vector3d dea1 = ea1_t1 - ea1_t0;
+    Eigen::Vector3d deb0 = eb0_t1 - eb0_t0;
+    Eigen::Vector3d deb1 = eb1_t1 - eb1_t0;
+    subtract_mean(dea0, dea1, deb0, deb1);
+
+    const double max_disp_mag =
+        std::sqrt(std::max(dea0.squaredNorm(), dea1.squaredNorm()))
+        + std::sqrt(std::max(deb0.squaredNorm(), deb1.squaredNorm()));
+    if (max_disp_mag == 0) {
+        return false;
+    }
+
+    VectorMax12d x = stack(ea0_t0, ea1_t0, eb0_t0, eb1_t0);
+    const VectorMax12d dx = stack(dea0, dea1, deb0, deb1);
+
+    const double min_distance_sq = min_distance * min_distance;
+    auto distance_squared = [min_distance_sq](const VectorMax12d& x) {
+        const auto& ea0 = x.head<3>();
+        const auto& ea1 = x.segment<3>(3);
+        const auto& eb0 = x.segment<3>(6);
+        const auto& eb1 = x.tail<3>();
+
+        double d_sq = edge_edge_distance(ea0, ea1, eb0, eb1);
+        if (d_sq - min_distance_sq <= 0) {
+            // since we ensured other place that all dist smaller than d̂ are
+            // positive, this must be some far away nearly parallel edges
+            d_sq = std::min(
+                { (ea0 - eb0).squaredNorm(), (ea0 - eb1).squaredNorm(),
+                  (ea1 - eb0).squaredNorm(), (ea1 - eb1).squaredNorm() });
+        }
+        return d_sq;
+    };
+
+    return additive_ccd(
+        x, dx, distance_squared, max_disp_mag, toi, min_distance, tmax,
+        conservative_rescaling);
+}
+
+} // namespace ipc::additive_ccd
\ No newline at end of file
diff --git a/src/ipc/ccd/additive_ccd.hpp b/src/ipc/ccd/additive_ccd.hpp
new file mode 100644
index 000000000..bef785e55
--- /dev/null
+++ b/src/ipc/ccd/additive_ccd.hpp
@@ -0,0 +1,138 @@
+//
+// NOTE: These methods are provided for reference comparison with [Li et al.
+// 2021] and is not utilized by the high-level functionality. In compairson to
+// Tight Inclusion CCD, this CCD method is not provably conservative and so can
+// potentially produce false negatives (i.e., miss collisions) due to
+// floating-point rounding error. However, it is much faster than Tight
+// Inclusion CCD (>100×) and very robust due to the gaps and conservative
+// rescaling used.
+//
+
+#pragma once
+
+#include <ipc/utils/eigen_ext.hpp>
+
+namespace ipc::additive_ccd {
+
+/// The default conservative rescaling value used to avoid taking steps exactly
+/// to impact. Value choosen to based on [Li et al. 2021].
+static constexpr double DEFAULT_CCD_CONSERVATIVE_RESCALING = 0.9;
+
+/// @brief Computes the time of impact between two points using continuous collision detection.
+/// @param p0_t0 The initial position of the first point.
+/// @param p1_t0 The initial position of the second point.
+/// @param p0_t1 The final position of the first point.
+/// @param p1_t1 The final position of the second point.
+/// @param[out] toi The time of impact between the two points.
+/// @param min_distance The minimum distance between two objects.
+/// @param tmax The maximum time to check for collisions.
+/// @param conservative_rescaling The conservative rescaling of the time of impact.
+/// @return True if a collision was detected, false otherwise.
+bool point_point_ccd(
+    const VectorMax3d& p0_t0,
+    const VectorMax3d& p1_t0,
+    const VectorMax3d& p0_t1,
+    const VectorMax3d& p1_t1,
+    double& toi,
+    const double min_distance = 0.0,
+    const double tmax = 1.0,
+    const double conservative_rescaling = DEFAULT_CCD_CONSERVATIVE_RESCALING);
+
+/// @brief Computes the time of impact between a point and an edge using continuous collision detection.
+/// @param p_t0 The initial position of the point.
+/// @param e0_t0 The initial position of the first endpoint of the edge.
+/// @param e1_t0 The initial position of the second endpoint of the edge.
+/// @param p_t1 The final position of the point.
+/// @param e0_t1 The final position of the first endpoint of the edge.
+/// @param e1_t1 The final position of the second endpoint of the edge.
+/// @param[out] toi The time of impact between the point and the edge.
+/// @param min_distance The minimum distance between two objects.
+/// @param tmax The maximum time to check for collisions.
+/// @param conservative_rescaling The conservative rescaling of the time of impact.
+/// @return True if a collision was detected, false otherwise.
+bool point_edge_ccd(
+    const VectorMax3d& p_t0,
+    const VectorMax3d& e0_t0,
+    const VectorMax3d& e1_t0,
+    const VectorMax3d& p_t1,
+    const VectorMax3d& e0_t1,
+    const VectorMax3d& e1_t1,
+    double& toi,
+    const double min_distance = 0.0,
+    const double tmax = 1.0,
+    const double conservative_rescaling = DEFAULT_CCD_CONSERVATIVE_RESCALING);
+
+/// @brief Computes the time of impact between a point and a triangle using continuous collision detection.
+/// @param p_t0 The initial position of the point.
+/// @param t0_t0 The initial position of the first vertex of the triangle.
+/// @param t1_t0 The initial position of the second vertex of the triangle.
+/// @param t2_t0 The initial position of the third vertex of the triangle.
+/// @param p_t1 The final position of the point.
+/// @param t0_t1 The final position of the first vertex of the triangle.
+/// @param t1_t1 The final position of the second vertex of the triangle.
+/// @param t2_t1 The final position of the third vertex of the triangle.
+/// @param[out] toi The time of impact between the point and the triangle.
+/// @param min_distance The minimum distance between two objects.
+/// @param tmax The maximum time to check for collisions.
+/// @param conservative_rescaling The conservative rescaling of the time of impact.
+/// @return True if a collision was detected, false otherwise.
+bool point_triangle_ccd(
+    const Eigen::Vector3d& p_t0,
+    const Eigen::Vector3d& t0_t0,
+    const Eigen::Vector3d& t1_t0,
+    const Eigen::Vector3d& t2_t0,
+    const Eigen::Vector3d& p_t1,
+    const Eigen::Vector3d& t0_t1,
+    const Eigen::Vector3d& t1_t1,
+    const Eigen::Vector3d& t2_t1,
+    double& toi,
+    const double min_distance = 0.0,
+    const double tmax = 1.0,
+    const double conservative_rescaling = DEFAULT_CCD_CONSERVATIVE_RESCALING);
+
+/// @brief Computes the time of impact between two edges using continuous collision detection.
+/// @param ea0_t0 The initial position of the first endpoint of the first edge.
+/// @param ea1_t0 The initial position of the second endpoint of the first edge.
+/// @param eb0_t0 The initial position of the first endpoint of the second edge.
+/// @param eb1_t0 The initial position of the second endpoint of the second edge.
+/// @param ea0_t1 The final position of the first endpoint of the first edge.
+/// @param ea1_t1 The final position of the second endpoint of the first edge.
+/// @param eb0_t1 The final position of the first endpoint of the second edge.
+/// @param eb1_t1 The final position of the second endpoint of the second edge.
+/// @param[out] toi The time of impact between the two edges.
+/// @param min_distance The minimum distance between two objects.
+/// @param tmax The maximum time to check for collisions.
+/// @param conservative_rescaling The conservative rescaling of the time of impact.
+/// @return True if a collision was detected, false otherwise.
+bool edge_edge_ccd(
+    const Eigen::Vector3d& ea0_t0,
+    const Eigen::Vector3d& ea1_t0,
+    const Eigen::Vector3d& eb0_t0,
+    const Eigen::Vector3d& eb1_t0,
+    const Eigen::Vector3d& ea0_t1,
+    const Eigen::Vector3d& ea1_t1,
+    const Eigen::Vector3d& eb0_t1,
+    const Eigen::Vector3d& eb1_t1,
+    double& toi,
+    const double min_distance = 0.0,
+    const double tmax = 1.0,
+    const double conservative_rescaling = DEFAULT_CCD_CONSERVATIVE_RESCALING);
+
+/// @brief Computes the time of impact between two objects using additive continuous collision detection.
+/// @param distance_squared A function that computes the squared distance between the two objects at a given time.
+/// @param[out] toi The time of impact between the two objects.
+/// @param min_distance The minimum distance between the objects.
+/// @param tmax The maximum time to check for collisions.
+/// @param conservative_rescaling The amount to rescale the objects by to ensure conservative advancement.
+/// @return True if a collision was detected, false otherwise.
+bool additive_ccd(
+    VectorMax12d x,
+    const VectorMax12d& dx,
+    const std::function<double()>& distance_squared,
+    const double max_disp_mag,
+    double& toi,
+    const double min_distance = 0.0,
+    const double tmax = 1.0,
+    const double conservative_rescaling = DEFAULT_CCD_CONSERVATIVE_RESCALING);
+
+} // namespace ipc::additive_ccd
\ No newline at end of file
diff --git a/src/ipc/ccd/ccd.hpp b/src/ipc/ccd/ccd.hpp
index d3082ee2a..edaae03ff 100644
--- a/src/ipc/ccd/ccd.hpp
+++ b/src/ipc/ccd/ccd.hpp
@@ -14,6 +14,20 @@ static constexpr double DEFAULT_CCD_CONSERVATIVE_RESCALING = 0.8;
 
 // 2D
 
+/// @brief Computes the time of impact between a point and an edge in 2D using continuous collision detection.
+/// @param p_t0 The initial position of the point.
+/// @param e0_t0 The initial position of the first endpoint of the edge.
+/// @param e1_t0 The initial position of the second endpoint of the edge.
+/// @param p_t1 The final position of the point.
+/// @param e0_t1 The final position of the first endpoint of the edge.
+/// @param e1_t1 The final position of the second endpoint of the edge.
+/// @param[out] toi The time of impact between the point and the edge.
+/// @param min_distance The minimum distance between the objects.
+/// @param tmax The maximum time to check for collisions.
+/// @param tolerance The error tolerance for the time of impact.
+/// @param max_iterations The maximum number of iterations to perform.
+/// @param conservative_rescaling The conservative rescaling of the time of impact.
+/// @return True if a collision was detected, false otherwise.
 bool point_edge_ccd_2D(
     const Eigen::Vector2d& p_t0,
     const Eigen::Vector2d& e0_t0,
@@ -30,6 +44,18 @@ bool point_edge_ccd_2D(
 
 // 3D
 
+/// @brief Computes the time of impact between two points in 3D using continuous collision detection.
+/// @param p0_t0 The initial position of the first point.
+/// @param p1_t0 The initial position of the second point.
+/// @param p0_t1 The final position of the first point.
+/// @param p1_t1 The final position of the second point.
+/// @param[out] toi The time of impact between the two points.
+/// @param min_distance The minimum distance between the objects.
+/// @param tmax The maximum time to check for collisions.
+/// @param tolerance The error tolerance for the time of impact.
+/// @param max_iterations The maximum number of iterations to perform.
+/// @param conservative_rescaling The conservative rescaling of the time of impact.
+/// @return True if a collision was detected, false otherwise.
 bool point_point_ccd(
     const Eigen::Vector3d& p0_t0,
     const Eigen::Vector3d& p1_t0,
@@ -42,6 +68,20 @@ bool point_point_ccd(
     const long max_iterations = DEFAULT_CCD_MAX_ITERATIONS,
     const double conservative_rescaling = DEFAULT_CCD_CONSERVATIVE_RESCALING);
 
+/// @brief Computes the time of impact between a point and an edge in 3D using continuous collision detection.
+/// @param p_t0 The initial position of the point.
+/// @param e0_t0 The initial position of the first endpoint of the edge.
+/// @param e1_t0 The initial position of the second endpoint of the edge.
+/// @param p_t1 The final position of the point.
+/// @param e0_t1 The final position of the first endpoint of the edge.
+/// @param e1_t1 The final position of the second endpoint of the edge.
+/// @param[out] toi The time of impact between the point and the edge.
+/// @param min_distance The minimum distance between the objects.
+/// @param tmax The maximum time to check for collisions.
+/// @param tolerance The error tolerance for the time of impact.
+/// @param max_iterations The maximum number of iterations to perform.
+/// @param conservative_rescaling The conservative rescaling of the time of impact.
+/// @return True if a collision was detected, false otherwise.
 bool point_edge_ccd_3D(
     const Eigen::Vector3d& p_t0,
     const Eigen::Vector3d& e0_t0,
@@ -56,6 +96,22 @@ bool point_edge_ccd_3D(
     const long max_iterations = DEFAULT_CCD_MAX_ITERATIONS,
     const double conservative_rescaling = DEFAULT_CCD_CONSERVATIVE_RESCALING);
 
+/// @brief Computes the time of impact between a point and a triangle in 3D using continuous collision detection.
+/// @param p_t0 The initial position of the point.
+/// @param t0_t0 The initial position of the first vertex of the triangle.
+/// @param t1_t0 The initial position of the second vertex of the triangle.
+/// @param t2_t0 The initial position of the third vertex of the triangle.
+/// @param p_t1 The final position of the point.
+/// @param t0_t1 The final position of the first vertex of the triangle.
+/// @param t1_t1 The final position of the second vertex of the triangle.
+/// @param t2_t1 The final position of the third vertex of the triangle.
+/// @param[out] toi The time of impact between the point and the triangle.
+/// @param min_distance The minimum distance between the objects.
+/// @param tmax The maximum time to check for collisions.
+/// @param tolerance The error tolerance for the time of impact.
+/// @param max_iterations The maximum number of iterations to perform.
+/// @param conservative_rescaling The conservative rescaling of the time of impact.
+/// @return True if a collision was detected, false otherwise.
 bool point_triangle_ccd(
     const Eigen::Vector3d& p_t0,
     const Eigen::Vector3d& t0_t0,
@@ -72,6 +128,22 @@ bool point_triangle_ccd(
     const long max_iterations = DEFAULT_CCD_MAX_ITERATIONS,
     const double conservative_rescaling = DEFAULT_CCD_CONSERVATIVE_RESCALING);
 
+/// @brief Computes the time of impact between two edges in 3D using continuous collision detection.
+/// @param ea0_t0 The initial position of the first endpoint of the first edge.
+/// @param ea1_t0 The initial position of the second endpoint of the first edge.
+/// @param eb0_t0 The initial position of the first endpoint of the second edge.
+/// @param eb1_t0 The initial position of the second endpoint of the second edge.
+/// @param ea0_t1 The final position of the first endpoint of the first edge.
+/// @param ea1_t1 The final position of the second endpoint of the first edge.
+/// @param eb0_t1 The final position of the first endpoint of the second edge.
+/// @param eb1_t1 The final position of the second endpoint of the second edge.
+/// @param[out] toi The time of impact between the two edges.
+/// @param min_distance The minimum distance between the objects.
+/// @param tmax The maximum time to check for collisions.
+/// @param tolerance The error tolerance for the time of impact.
+/// @param max_iterations The maximum number of iterations to perform.
+/// @param conservative_rescaling The conservative rescaling of the time of impact.
+/// @return True if a collision was detected, false otherwise.
 bool edge_edge_ccd(
     const Eigen::Vector3d& ea0_t0,
     const Eigen::Vector3d& ea1_t0,
@@ -90,6 +162,20 @@ bool edge_edge_ccd(
 
 // 2D or 3D
 
+/// @brief Computes the time of impact between a point and an edge in 2D or 3D using continuous collision detection.
+/// @param p_t0 The initial position of the point.
+/// @param e0_t0 The initial position of the first endpoint of the edge.
+/// @param e1_t0 The initial position of the second endpoint of the edge.
+/// @param p_t1 The final position of the point.
+/// @param e0_t1 The final position of the first endpoint of the edge.
+/// @param e1_t1 The final position of the second endpoint of the edge.
+/// @param[out] toi The time of impact between the point and the edge.
+/// @param min_distance The minimum distance between the objects.
+/// @param tmax The maximum time to check for collisions.
+/// @param tolerance The error tolerance for the time of impact.
+/// @param max_iterations The maximum number of iterations to perform.
+/// @param conservative_rescaling The conservative rescaling of the time of impact.
+/// @return True if a collision was detected, false otherwise.
 bool point_edge_ccd(
     const VectorMax3d& p_t0,
     const VectorMax3d& e0_t0,
diff --git a/src/ipc/utils/CMakeLists.txt b/src/ipc/utils/CMakeLists.txt
index ee2e471b1..9c14a00ef 100644
--- a/src/ipc/utils/CMakeLists.txt
+++ b/src/ipc/utils/CMakeLists.txt
@@ -9,7 +9,6 @@ set(SOURCES
   logger.cpp
   logger.hpp
   merge_thread_local.hpp
-  rational.hpp
   save_obj.cpp
   save_obj.hpp
   unordered_map_and_set.cpp
diff --git a/src/ipc/utils/intersection.cpp b/src/ipc/utils/intersection.cpp
index 8b4f0fea3..f7fd194b1 100644
--- a/src/ipc/utils/intersection.cpp
+++ b/src/ipc/utils/intersection.cpp
@@ -2,14 +2,15 @@
 
 #include <ipc/utils/eigen_ext.hpp>
 #include <ipc/config.hpp>
-#ifdef IPC_TOOLKIT_WITH_RATIONAL_INTERSECTION
-#include <ipc/utils/rational.hpp>
-#endif
 
 #include <igl/predicates/predicates.h>
 
 #include <Eigen/Geometry>
 
+#ifdef IPC_TOOLKIT_WITH_RATIONAL_INTERSECTION
+#include <rational/rational.hpp>
+#endif
+
 namespace ipc {
 
 #ifdef IPC_TOOLKIT_WITH_RATIONAL_INTERSECTION
@@ -20,7 +21,12 @@ bool is_edge_intersecting_triangle_rational(
     const Eigen::Vector3d& t1_float,
     const Eigen::Vector3d& t2_float)
 {
-    Vector3<Rational> e0, e1, t0, t1, t2;
+    using namespace rational;
+
+    typedef Eigen::Matrix<Rational, 3, 1, Eigen::ColMajor | Eigen::DontAlign>
+        Vector3r;
+
+    Vector3r e0, e1, t0, t1, t2;
 
     for (int d = 0; d < 3; ++d) {
         e0[d] = e0_float[d];
@@ -44,7 +50,7 @@ bool is_edge_intersecting_triangle_rational(
         - e1[1] * t1[2] * t2[0] - e1[2] * t0[0] * t1[1] + e1[2] * t0[0] * t2[1]
         + e1[2] * t0[1] * t1[0] - e1[2] * t0[1] * t2[0] - e1[2] * t1[0] * t2[1]
         + e1[2] * t1[1] * t2[0];
-    if (d.get_sign() == 0) {
+    if (d.sign() == 0) {
         return true;
     }
 
diff --git a/src/ipc/utils/rational.hpp b/src/ipc/utils/rational.hpp
deleted file mode 100644
index 23c178500..000000000
--- a/src/ipc/utils/rational.hpp
+++ /dev/null
@@ -1,138 +0,0 @@
-#pragma once
-
-#include <ipc/config.hpp>
-
-#ifdef IPC_TOOLKIT_WITH_RATIONAL_INTERSECTION
-
-#include <gmp.h>
-#include <iostream>
-
-namespace ipc {
-
-// https://cs.nyu.edu/acsys/cvc3/releases/1.5/doc/rational-gmp_8cpp-source.html
-class Rational {
-public:
-    mpq_t value;
-    void canonicalize() { mpq_canonicalize(value); }
-    int get_sign() const { return mpq_sgn(value); }
-
-    Rational()
-    {
-        mpq_init(value);
-        mpq_set_d(value, 0);
-    }
-
-    Rational(double d)
-    {
-        mpq_init(value);
-        mpq_set_d(value, d);
-        // canonicalize();
-    }
-
-    Rational(const mpq_t& v_)
-    {
-        mpq_init(value);
-        mpq_set(value, v_);
-        // canonicalize();
-    }
-
-    Rational(const Rational& other)
-    {
-        mpq_init(value);
-        mpq_set(value, other.value);
-    }
-
-    ~Rational() { mpq_clear(value); }
-
-    friend Rational operator-(const Rational& v)
-    {
-        Rational r_out;
-        mpq_neg(r_out.value, v.value);
-        return r_out;
-    }
-
-    friend Rational operator+(const Rational& x, const Rational& y)
-    {
-        Rational r_out;
-        mpq_add(r_out.value, x.value, y.value);
-        return r_out;
-    }
-
-    friend Rational operator-(const Rational& x, const Rational& y)
-    {
-        Rational r_out;
-        mpq_sub(r_out.value, x.value, y.value);
-        return r_out;
-    }
-
-    friend Rational operator*(const Rational& x, const Rational& y)
-    {
-        Rational r_out;
-        mpq_mul(r_out.value, x.value, y.value);
-        return r_out;
-    }
-
-    friend Rational operator/(const Rational& x, const Rational& y)
-    {
-        Rational r_out;
-        mpq_div(r_out.value, x.value, y.value);
-        return r_out;
-    }
-
-    Rational& operator=(const Rational& x)
-    {
-        if (this == &x)
-            return *this;
-        mpq_set(value, x.value);
-        return *this;
-    }
-
-    Rational& operator=(const double x)
-    {
-        mpq_set_d(value, x);
-        // canonicalize();
-        return *this;
-    }
-
-    friend bool operator<(const Rational& r, const Rational& r1)
-    {
-        return mpq_cmp(r.value, r1.value) < 0;
-    }
-
-    friend bool operator>(const Rational& r, const Rational& r1)
-    {
-        return mpq_cmp(r.value, r1.value) > 0;
-    }
-
-    friend bool operator<=(const Rational& r, const Rational& r1)
-    {
-        return mpq_cmp(r.value, r1.value) <= 0;
-    }
-
-    friend bool operator>=(const Rational& r, const Rational& r1)
-    {
-        return mpq_cmp(r.value, r1.value) >= 0;
-    }
-
-    friend bool operator==(const Rational& r, const Rational& r1)
-    {
-        return mpq_equal(r.value, r1.value);
-    }
-
-    friend bool operator!=(const Rational& r, const Rational& r1)
-    {
-        return !mpq_equal(r.value, r1.value);
-    }
-
-    double to_double() { return mpq_get_d(value); }
-
-    friend std::ostream& operator<<(std::ostream& os, const Rational& r)
-    {
-        os << mpq_get_d(r.value);
-        return os;
-    }
-};
-
-} // namespace ipc
-
-#endif


Test patch:
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 71fac9fdd..49f658eab 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -17,6 +17,11 @@ add_executable(ipc_toolkit_tests
   ccd/benchmark_ccd.cpp
   ccd/collision_generator.cpp
   ccd/test_ccd.cpp
+  ccd/test_ccd_benchmark.cpp
+  ccd/test_point_edge_ccd.cpp
+  ccd/test_point_triangle_ccd.cpp
+  ccd/test_edge_edge_ccd.cpp
+  ccd/test_ccd.cpp
   test_cfl.cpp
 
   # Test distances and dertivatives
@@ -68,6 +73,20 @@ target_link_libraries(ipc_toolkit_tests PUBLIC finitediff::finitediff)
 include(json)
 target_link_libraries(ipc_toolkit_tests PUBLIC nlohmann_json::nlohmann_json)
 
+if (IPC_TOOLKIT_TEST_CCD_BENCHMARK)
+  include(ccd_query_io)
+  target_link_libraries(ipc_toolkit_tests PUBLIC ccd_io::ccd_io)
+  target_compile_definitions(ipc_toolkit_tests PRIVATE IPC_TOOLKIT_TEST_CCD_BENCHMARK)
+  if(NOT (IPC_TOOLKIT_CCD_BENCHMARK_DIR STREQUAL ""))
+    message(STATUS "Using CCD benchmark directory: ${IPC_TOOLKIT_CCD_BENCHMARK_DIR}")
+    target_compile_definitions(ipc_toolkit_tests PRIVATE IPC_TOOLKIT_CCD_BENCHMARK_DIR="${IPC_TOOLKIT_CCD_BENCHMARK_DIR}")
+  endif()
+  if(NOT (IPC_TOOLKIT_CCD_NEW_BENCHMARK_DIR STREQUAL ""))
+    message(STATUS "Using new CCD benchmark directory: ${IPC_TOOLKIT_CCD_NEW_BENCHMARK_DIR}")
+    target_compile_definitions(ipc_toolkit_tests PRIVATE IPC_TOOLKIT_CCD_NEW_BENCHMARK_DIR="${IPC_TOOLKIT_CCD_NEW_BENCHMARK_DIR}")
+  endif()
+endif()
+
 if(IPC_TOOLKIT_WITH_CUDA)
   find_package(CUDAToolkit)
   target_link_libraries(ipc_toolkit_tests PRIVATE CUDA::cudart)
diff --git a/tests/barrier/test_barrier.cpp b/tests/barrier/test_barrier.cpp
index a82aea7ba..3e10045b9 100644
--- a/tests/barrier/test_barrier.cpp
+++ b/tests/barrier/test_barrier.cpp
@@ -45,7 +45,7 @@ double normalized_barrier_hessian(const double d, const double dhat)
 }
 } // namespace
 
-TEST_CASE("Test barrier derivatives", "[barrier]")
+TEST_CASE("Barrier derivatives", "[barrier]")
 {
     bool use_dist_sqr = GENERATE(false, true);
     double dhat = GENERATE_COPY(range(use_dist_sqr ? -2 : -5, 0));
@@ -115,7 +115,7 @@ TEST_CASE("Test barrier derivatives", "[barrier]")
     CHECK(fd::compare_gradient(fgrad, grad));
 }
 
-TEST_CASE("Test physical barrier", "[barrier]")
+TEST_CASE("Physical barrier", "[barrier]")
 {
     const bool use_dist_sqr = GENERATE(false, true);
     const double dhat =
diff --git a/tests/broad_phase/test_broad_phase.cpp b/tests/broad_phase/test_broad_phase.cpp
index 260ad6709..4765b5728 100644
--- a/tests/broad_phase/test_broad_phase.cpp
+++ b/tests/broad_phase/test_broad_phase.cpp
@@ -86,9 +86,9 @@ TEST_CASE("Vertex-Vertex Broad Phase", "[ccd][broad_phase][2D]")
 }
 
 #if defined(NDEBUG) || !(defined(WIN32) || defined(_WIN32) || defined(__WIN32))
-TEST_CASE("Entire 2D Mesh", "[ccd][broad_phase][2D]")
+TEST_CASE("Broad Phase: 2D Mesh", "[ccd][broad_phase][2D]")
 #else
-TEST_CASE("Entire 2D Mesh", "[ccd][broad_phase][2D][.]")
+TEST_CASE("Broad Phase: 2D Mesh", "[ccd][broad_phase][2D][.]")
 #endif
 {
     Eigen::MatrixXd tmp;
@@ -116,7 +116,7 @@ TEST_CASE("Entire 2D Mesh", "[ccd][broad_phase][2D][.]")
 }
 
 TEST_CASE(
-    "Test constraint_set.build() with codimensional points",
+    "Build constraint set with codimensional points",
     "[broad_phase][constraints]")
 {
     const double dhat = 1e-3;
diff --git a/tests/broad_phase/test_spatial_hash.cpp b/tests/broad_phase/test_spatial_hash.cpp
index e1ea32736..7ad8fabea 100644
--- a/tests/broad_phase/test_spatial_hash.cpp
+++ b/tests/broad_phase/test_spatial_hash.cpp
@@ -12,7 +12,7 @@
 
 using namespace ipc;
 
-TEST_CASE("Test build SpatialHash", "[spatial_hash][build]")
+TEST_CASE("Build SpatialHash", "[spatial_hash][build]")
 {
     Eigen::MatrixXd V0, V1;
     Eigen::MatrixXi E, F;
diff --git a/tests/ccd/test_ccd.cpp b/tests/ccd/test_ccd.cpp
index 78069fc51..0314f0b8a 100644
--- a/tests/ccd/test_ccd.cpp
+++ b/tests/ccd/test_ccd.cpp
@@ -2,273 +2,16 @@
 
 #include <ipc/ipc.hpp>
 #include <ipc/ccd/ccd.hpp>
+#include <ipc/ccd/additive_ccd.hpp>
 #include <ipc/ccd/point_static_plane.hpp>
 
 #include <test_utils.hpp>
 
-#include "collision_generator.hpp"
-
 #include <igl/read_triangle_mesh.h>
 #include <igl/edges.h>
 
 using namespace ipc;
 
-static const double EPSILON = std::numeric_limits<float>::epsilon();
-
-TEST_CASE("Point-edge 2D CCD", "[ccd]")
-{
-    Eigen::Vector2d p_t0, p_t1, e0_t0, e0_t1, e1_t0, e1_t1;
-    double toi_expected;
-    bool is_collision_expected;
-    SECTION("Edge becomes degenerate")
-    {
-        p_t0 << 0, 1;
-        e0_t0 << -1, 0;
-        e1_t0 << 1, 0;
-        SECTION("Edge degenerates before impact")
-        {
-            p_t1 << 0, -1;
-            e0_t1 << 3, 0;
-            e1_t1 << -3, 0;
-            // The edge will become degenerate at t=0.25
-        }
-        SECTION("Edge degenerates after impact")
-        {
-            p_t1 << 0, -1;
-            e0_t1 << 0.5, 0;
-            e1_t1 << -0.5, 0;
-            // The edge will become degenerate at t=2/3
-        }
-        // The point will collide with the edge at t=0.5
-        is_collision_expected = true;
-        toi_expected = 0.5;
-    }
-    SECTION("Edge moving right; point moving left")
-    {
-        p_t0 << -1, 0;
-        e0_t0 << 1, -1;
-        e1_t0 << 1, 1;
-
-        p_t1 << 1, 0;
-        e0_t1 << -1, -1;
-        e1_t1 << -1, 1;
-
-        is_collision_expected = true;
-        toi_expected = 0.5;
-    }
-    SECTION("Point on edge's line moving towards edge")
-    {
-        p_t0 << 0, 0;
-        e0_t0 << 0, 1;
-        e1_t0 << 0, 2;
-
-        p_t1 << 0, 2;
-        e0_t1 << 0, 1;
-        e1_t1 << 0, 2;
-
-        is_collision_expected = true;
-        toi_expected = 0.5;
-    }
-    SECTION("Point and edge moving parallel")
-    {
-        p_t0 << 0, 1;
-        e0_t0 << 1, 0;
-        e1_t0 << 1, 2;
-
-        p_t1 << 0, 2;
-        e0_t1 << 1, 1;
-        e1_t1 << 1, 3;
-
-        is_collision_expected = false;
-    }
-    SECTION("Point moving right; edge stretching vertically")
-    {
-        e1_t0 << 1, -1;
-        e1_t1 << 1, -2;
-        SECTION("Swap vertices order e_0 = [0, 2]")
-        {
-            p_t0 << 1, 1;
-            e0_t0 << 0, 0;
-
-            p_t1 << 1, 2;
-            e0_t1 << 1, 0;
-
-            std::swap(p_t0, e0_t0);
-            std::swap(p_t1, e0_t1);
-        }
-        SECTION("Swap vertices order e_0 = [1, 2]")
-        {
-            p_t0 << 0, 0;
-            e0_t0 << 1, 1;
-
-            p_t1 << 1, 0;
-            e0_t1 << 1, 2;
-        }
-        is_collision_expected = true;
-        toi_expected = 1.0;
-    }
-    SECTION("Point-point")
-    {
-        p_t0 << 1.11111, 0.5;
-        p_t1 << 0.888889, 0.5;
-        e0_t0 << 1, 0.5;
-        e1_t0 << 1, 0.75;
-        e0_t1 = e0_t0;
-        e1_t1 = e1_t0;
-
-        is_collision_expected = true;
-        toi_expected = 0.5;
-    }
-
-    double toi, alpha;
-    bool is_colliding = point_edge_ccd(
-        p_t0, e0_t0, e1_t0, p_t1, e0_t1, e1_t1, toi,
-        /*min_distance=*/0.0,
-        /*tmax=*/1.0,
-        /*tolerance=*/1e-6,
-        /*max_iterations=*/1e7,
-        /*conservative_rescaling=*/1.0);
-    REQUIRE(is_colliding == is_collision_expected);
-    if (is_collision_expected) {
-        CAPTURE(
-            p_t0.transpose(), e0_t0.transpose(), e1_t0.transpose(),
-            p_t1.transpose(), e0_t1.transpose(), e1_t1.transpose());
-        CHECK(toi <= toi_expected);
-    }
-}
-
-// ---------------------------------------------------
-// Helper functions
-// ---------------------------------------------------
-
-/// Compares the time of impact of different implementations
-/// against the expected time of impact
-void check_toi(
-    const Eigen::Vector2d& p_t0,
-    const Eigen::Vector2d& e0_t0,
-    const Eigen::Vector2d& e1_t0,
-    const Eigen::Vector2d& p_t1,
-    const Eigen::Vector2d& e0_t1,
-    const Eigen::Vector2d& e1_t1,
-    const double toi_expected)
-{
-    double toi_actual;
-    // check autodiff code
-    toi_actual = -1.0;
-    bool has_collision = point_edge_ccd(
-        p_t0, e0_t0, e1_t0, p_t1, e0_t1, e1_t1, toi_actual,
-        /*min_distance=*/0.0,
-        /*tmax=*/1.0,
-        /*tolerance=*/1e-6,
-        /*max_iterations=*/1e7,
-        /*conservative_rescaling=*/1.0);
-    CHECK(has_collision);
-    CHECK(toi_actual <= toi_expected);
-}
-
-// ---------------------------------------------------
-// Tests
-// ---------------------------------------------------
-#ifdef IPC_TOOLKIT_WITH_CORRECT_CCD
-TEST_CASE("Point-edge 2D ToI", "[ccd][toi]")
-#else
-TEST_CASE("Point-edge 2D ToI", "[ccd][toi][!mayfail]")
-#endif
-{
-    Eigen::Vector2d p_t0, e0_t0, e1_t0;
-    Eigen::Vector2d dp, de0, de1;
-
-    double expected_toi = -1;
-
-    SECTION("perpendicular") // (alpha=0.5)
-    {
-        p_t0 << 0.0, 1.0;
-        e0_t0 << -1.0, 0.0;
-        e1_t0 << 1.0, 0.0;
-
-        // touches, intersects, passes-trough
-        auto vel = GENERATE(1.0 + 1e-6, 2.0, 4.0);
-        // moving direction:
-        //  -> ->, -> -, -> <-, - <-, <- <-
-        auto j = GENERATE(0, 1, 2, 3, 4);
-        // extension, no-deform, compression,
-        auto dx = GENERATE(0.5, 0.0, -0.5);
-
-        dp << 0.0, -(3 - j) * vel / 2.0;
-        de0 << -dx, (j - 1.0) * vel / 2.0;
-        de1 << dx, (j - 1.0) * vel / 2.0;
-
-        expected_toi = 1.0 / vel;
-    }
-
-    SECTION("double impact") // (rotating edge)
-    {
-        e0_t0 << -1.0, 0.0;
-        e1_t0 << 1.0, 0.0;
-        p_t0 << 0.0, 0.5;
-
-        de0 << 1.6730970740318298, 0.8025388419628143;
-        de1 << -1.616142749786377, -0.6420311331748962;
-        dp << 0.0, -1.0;
-
-        expected_toi = 0.4482900963;
-    }
-
-    SECTION("random")
-    {
-        // This times out on Windows debug
-#if defined(NDEBUG) || !(defined(WIN32) || defined(_WIN32) || defined(__WIN32))
-        auto impact = GENERATE(random_impacts(100, /*rigid=*/true));
-#else
-        auto impact = GENERATE(random_impacts(10, /*rigid=*/true));
-#endif
-
-        p_t0 = impact.p_t0;
-        e0_t0 = impact.e0_t0;
-        e1_t0 = impact.e1_t0;
-
-        dp = impact.dp;
-        de0 = impact.de0;
-        de1 = impact.de1;
-
-        expected_toi = impact.toi;
-    }
-
-#ifdef NDEBUG           // This case takes forever to run
-    SECTION("parallel") //  (alpha=0 || alpha = 1)
-    {
-        p_t0 << 0.5, 0.0;
-        e0_t0 << -0.5, 0.0;
-        e1_t0 << -1.5, 0.0;
-
-        // touches, intersects, passes-trough
-        // double vel = GENERATE(1.0 + 1e-6, 2.0, 4.0);
-        double vel = 2.0;
-        // moving: both (same), ij, both (op), kl, both (same)
-        // double j = GENERATE(0, 1, 2, 3, 4);
-        double j = 0;
-        // extension, no-deform, compression,
-        // double dy = GENERATE(0.5, 0.0, -0.5);
-        double dy = 0.0;
-
-        dp << -(3 - j) * vel / 2.0, 0.0;
-        de0 << (j - 1.0) * vel / 2.0, 0.0;
-        // we only move one so we don't change the toi
-        de1 << (j - 1.0) * vel / 2.0, dy;
-
-        expected_toi = 1.0 / vel;
-    }
-#endif
-
-    Eigen::Vector2d p_t1 = p_t0 + dp;
-    Eigen::Vector2d e0_t1 = e0_t0 + de0;
-    Eigen::Vector2d e1_t1 = e1_t0 + de1;
-
-    check_toi(p_t0, e0_t0, e1_t0, p_t1, e0_t1, e1_t1, expected_toi);
-    // Flip the edges
-    check_toi(p_t0, e1_t0, e0_t0, p_t1, e1_t1, e0_t1, expected_toi);
-}
-
 TEST_CASE("Repeated CCD", "[ccd][repeat]")
 {
     const double FIRST_TOL = 1e-6, SECOND_TOL = 1e-7;
@@ -369,172 +112,6 @@ TEST_CASE("Repeated CCD", "[ccd][repeat]")
     } while (has_collisions_repeated && stepsize_repeated != 1.0);
 }
 
-// This test takes too long
-// TEST_CASE("Point-Triangle 3D CCD", "[ccd][point-triangle][timeout]")
-// {
-//     // point
-//     double v0z = GENERATE(0.0, -1.0);
-//     Eigen::Vector3d v0(0, 1, v0z);
-//     // triangle = (v1, v2, v3)
-//     Eigen::Vector3d v1(-1, 0, 1);
-//     Eigen::Vector3d v2(1, 0, 1);
-//     Eigen::Vector3d v3(0, 0, -1);
-//
-//     // displacements
-//     double u0y =
-//         -GENERATE(-1.0, 0.0, 0.5 - EPSILON, 0.5, 0.5 + EPSILON, 1.0, 2.0);
-//     double u0z = GENERATE(-EPSILON, 0.0, EPSILON);
-//     Eigen::Vector3d u0(0, u0y, u0z);
-//     double u1y =
-//         GENERATE(-1.0, 0.0, 0.5 - EPSILON, 0.5, 0.5 + EPSILON, 1.0, 2.0);
-//     Eigen::Vector3d u1(0, u1y, 0);
-//
-//     bool is_collision_expected = ((-u0y + u1y >= 1) && (v0z + u0z >=
-//     v3.z()));
-//
-//     double toi;
-//     bool is_colliding = point_triangle_ccd(
-//         v0, v1, v2, v3, v0 + u0, v1 + u1, v2 + u1, v3 + u1, toi);
-//
-//     CAPTURE(v0z, u0y, u1y, u0z, EPSILON);
-//     REQUIRE(is_colliding >= is_collision_expected);
-// }
-
-// Disable these tests in debug because they take a long time.
-#ifdef NDEBUG
-
-TEST_CASE("Edge-Edge 3D CCD", "[ccd][edge-edge]")
-{
-    // e0 = (v0, v1)
-    Eigen::Vector3d v0(-1, -1, 0);
-    Eigen::Vector3d v1(1, -1, 0);
-    // e2 = (v2, v3)
-    // double e1x = GENERATE(
-    //     -1 - EPSILON, -1, -1 + EPSILON, -0.5, 0, 0.5, 1 - EPSILON, 1,
-    //     1 + EPSILON);
-    double e1x = GENERATE(0);
-    Eigen::Vector3d v2(e1x, 1, -1);
-    Eigen::Vector3d v3(e1x, 1, 1);
-
-    // displacements
-    double y_displacement =
-        GENERATE(-1.0, 0.0, 1 - EPSILON, 1.0, 1 + EPSILON, 2.0);
-
-    Eigen::Vector3d u0, u1;
-    bool is_collision_expected;
-    SECTION("moving")
-    {
-        u0 << 0, y_displacement, 0;
-        u1 << 0, -y_displacement, 0;
-        is_collision_expected = y_displacement >= 1.0 && e1x >= -1 && e1x <= 1;
-    }
-    SECTION("fixed")
-    {
-        u0 << 0, 2 * y_displacement, 0;
-        u1.setZero();
-        is_collision_expected = y_displacement >= 2.0 && e1x >= -1 && e1x <= 1;
-    }
-
-    double toi;
-    bool is_colliding =
-        edge_edge_ccd(v0, v1, v2, v3, v0 + u0, v1 + u0, v2 + u1, v3 + u1, toi);
-
-    CAPTURE(y_displacement, e1x);
-    CHECK(is_colliding >= is_collision_expected);
-}
-
-// Only run this slow case on Linux and macOS
-#if !defined(WIN32)
-
-TEST_CASE("Zhongshi test case", "[ccd][point-triangle][zhongshi]")
-{
-    double qy = GENERATE(-EPSILON, 0, EPSILON);
-
-    Eigen::Vector3d q;
-    q << 0, qy, 0;
-
-    Eigen::Vector3d b0;
-    b0 << 0, 0, 0;
-    Eigen::Vector3d b1;
-    b1 << 0, 1, 0;
-    Eigen::Vector3d b2;
-    b2 << 1, 0, 0;
-
-    Eigen::Vector3d t0;
-    t0 << 0, 0, 1;
-    Eigen::Vector3d t1;
-    t1 << 0, 1, 1;
-    Eigen::Vector3d t2;
-    t2 << 1, 0, 1;
-
-    Eigen::Vector3d q1;
-    q1 << 0, qy, 0;
-
-    bool is_collision_expected = q.y() >= 0;
-
-    double toi;
-    bool is_colliding = point_triangle_ccd(q, b0, b1, b2, q1, t0, t1, t2, toi);
-
-    CAPTURE(qy);
-    CHECK(is_colliding >= is_collision_expected);
-}
-
-#endif
-
-TEST_CASE("Bolun test case", "[ccd][point-triangle][bolun]")
-{
-    Eigen::Vector3d x0(0.1, 0.1, 0.1), x1(0, 0, 1), x2(1, 0, 1), x3(0, 1, 1),
-        x0b(0.1, 0.1, 0.1), x1b(0, 0, 0), x2b(0, 1, 0), x3b(1, 0, 0);
-
-    bool is_collision_expected = true;
-
-    double toi;
-    bool is_colliding =
-        point_triangle_ccd(x0, x1, x2, x3, x0b, x1b, x2b, x3b, toi);
-
-    CHECK(is_colliding >= is_collision_expected);
-}
-
-TEST_CASE("Dobule root test case", "[ccd][edge-edge][double-root]")
-{
-    const Eigen::Vector3d a0s(-3.0022200, 0.2362580, 0.0165247);
-    const Eigen::Vector3d a1s(-3.2347850, 0.8312380, -0.1151003);
-    const Eigen::Vector3d a0e(-2.8995600, 0.0345838, 0.0638580);
-    const Eigen::Vector3d a1e(-3.1716930, 0.6104858, -0.0713340);
-    const Eigen::Vector3d b0(-3.0319900, 0.3148750, 0.0000000);
-    const Eigen::Vector3d b1(-2.8548800, 0.0900349, 0.0000000);
-
-    bool is_collision_expected = true;
-
-    double toi;
-    bool is_colliding = edge_edge_ccd(a0s, a1s, b0, b1, a0e, a1e, b0, b1, toi);
-
-    CHECK(is_colliding >= is_collision_expected);
-}
-
-TEST_CASE("Double root test case 2", "[ccd][edge-edge][double-root]")
-{
-    const Eigen::Vector3d a0s(0, 0, 1);
-    const Eigen::Vector3d a1s(0, 1, 1);
-    Eigen::Vector3d a0e(1, 1, 0);
-    Eigen::Vector3d a1e(0, 0, 0);
-    const Eigen::Vector3d b0(0.1, 0.2, 2);
-    const Eigen::Vector3d b1(0.1, 0.2, -1);
-
-    double t = GENERATE(0.5, 0.8, 0.88, 0.9, 1.0);
-    a0e = (a0e - a0s) * t + a0s;
-    a1e = (a1e - a1s) * t + a1s;
-
-    bool is_collision_expected = true;
-
-    double toi;
-    bool is_colliding = edge_edge_ccd(a0s, a1s, b0, b1, a0e, a1e, b0, b1, toi);
-
-    CHECK(is_colliding >= is_collision_expected);
-}
-
-#endif
-
 TEST_CASE("Point-Plane CCD", "[ccd][point-plane]")
 {
     Eigen::Vector3d p_t0 = Eigen::Vector3d::Random();
@@ -564,91 +141,6 @@ TEST_CASE("Point-Plane CCD", "[ccd][point-plane]")
     CHECK(toi <= t);
 }
 
-TEST_CASE("No Zero ToI CCD", "[ccd][no-zero-toi]")
-{
-    Eigen::Vector3d p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1;
-    p_t0 << 0.0133653, 0.100651, -0.0215935;
-    t0_t0 << 0.0100485, 0.0950896, -0.0171013;
-    t1_t0 << 0.0130388, 0.100666, -0.0218112;
-    t2_t0 << 0.015413, 0.100554, -0.0202265;
-    p_t1 << 0.0133652999767858, 0.099670000268615, -0.0215934999996444;
-    t0_t1 << 0.0100484999799995, 0.0941086002577558, -0.0171012999972189;
-    t1_t1 << 0.0130387999724314, 0.0996850002629403, -0.0218111999936902;
-    t2_t1 << 0.0154129999740718, 0.0995730002646605, -0.020226499996014;
-
-    double toi;
-    bool is_impacting = point_triangle_ccd(
-        p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1, toi);
-
-    CAPTURE(toi);
-    CHECK(!is_impacting);
-}
-
-#ifdef IPC_TOOLKIT_WITH_CORRECT_CCD
-TEST_CASE("Slow EE CCD", "[ccd][edge-edge][slow]")
-#else
-TEST_CASE("Slow EE CCD", "[ccd][edge-edge][slow][!mayfail]")
-#endif
-{
-    Eigen::Vector3d ea0_t0, ea1_t0, eb0_t0, eb1_t0, ea0_t1, ea1_t1, eb0_t1,
-        eb1_t1;
-    ea0_t0 << 1, 0.50803125, 2.10835646075301e-18;
-    ea1_t0 << -2.38233935445388e-18, 0.50803125, 1;
-    eb0_t0 << -4.99999999958867e-07, 0.5, 0;
-    eb1_t0 << -4.99999999958867e-07, 0.5, 1;
-    ea0_t1 << 1, 0.47124375, 4.11078309465837e-18;
-    ea1_t1 << -2.8526707189104e-18, 0.47124375, 1;
-    eb0_t1 << -4.99999999958867e-07, 0.5, 0;
-    eb1_t1 << -4.99999999958867e-07, 0.5, 1;
-
-    // BENCHMARK("compute toi")
-    // {
-    //     double toi;
-    //     bool is_impacting = edge_edge_ccd(
-    //         ea0_t0, ea1_t0, eb0_t0, eb1_t0, ea0_t1, ea1_t1, eb0_t1, eb1_t1,
-    //         toi);
-    // };
-
-    double tol = DEFAULT_CCD_TOLERANCE;
-    long max_iter = DEFAULT_CCD_MAX_ITERATIONS;
-    while (tol < 1) {
-        double toi;
-        bool is_impacting = edge_edge_ccd(
-            ea0_t0, ea1_t0, eb0_t0, eb1_t0, ea0_t1, ea1_t1, eb0_t1, eb1_t1, toi,
-            0.0, 1.0, tol, max_iter);
-        tol *= 10;
-
-        CAPTURE(toi);
-        CHECK(is_impacting);
-    }
-}
-
-TEST_CASE("Slow EE CCD 2", "[ccd][edge-edge][slow][thisone]")
-{
-    Eigen::Vector3d ea0_t0, ea1_t0, eb0_t0, eb1_t0, ea0_t1, ea1_t1, eb0_t1,
-        eb1_t1;
-    ea0_t0 << 1.00002232466453, 0.500004786049044, -2.06727783590977e-05;
-    ea1_t0 << 1.64687846177844e-05, 0.499996645067319, 1.63939999009028e-05;
-    eb0_t0 << 1, 0.5, 0;
-    eb1_t0 << 0, 0.5, 0;
-    ea0_t1 << 1.00294282700155, 0.498652627047143, 0.003626320742036;
-    ea1_t1 << -0.00219276550735626, 0.500871179186644, -0.00315828804921928;
-    eb0_t1 << 1, 0.5, 0;
-    eb1_t1 << 0, 0.5, 0;
-
-    double tol = 1e-6;
-    long max_iter = 1e6;
-    double tmax = 2.8076171875e-03;
-
-    double toi;
-    bool is_impacting = edge_edge_ccd(
-        ea0_t0, ea1_t0, eb0_t0, eb1_t0, ea0_t1, ea1_t1, eb0_t1, eb1_t1, toi,
-        0.0, tmax, tol, max_iter);
-
-    CAPTURE(toi);
-    CHECK(is_impacting);
-}
-
 TEST_CASE("Squash Tet", "[ccd]")
 {
     const double dhat = 1e-3;
diff --git a/tests/ccd/test_ccd_benchmark.cpp b/tests/ccd/test_ccd_benchmark.cpp
new file mode 100644
index 000000000..1f02ed67a
--- /dev/null
+++ b/tests/ccd/test_ccd_benchmark.cpp
@@ -0,0 +1,263 @@
+#include <ipc/config.hpp>
+
+#ifdef IPC_TOOLKIT_TEST_CCD_BENCHMARK
+#include <catch2/catch_all.hpp>
+
+#include <ipc/ccd/ccd.hpp>
+#include <ipc/ccd/additive_ccd.hpp>
+#include <fmt/format.h>
+
+#include <igl/Timer.h>
+#include <ccd_io/read_ccd_queries.hpp>
+
+#include <iostream>
+
+using namespace ipc;
+
+void run_benchmark(
+    const std::function<bool(
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        double&)>& edge_edge_ccd,
+    const std::function<bool(
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        const Eigen::Vector3d&,
+        double&)>& point_triangle_ccd)
+{
+    logger().set_level(spdlog::level::err);
+
+    namespace fs = std::filesystem;
+    using Matrix8x3 = Eigen::Matrix<double, 8, 3, Eigen::RowMajor>;
+
+    std::vector<fs::path> csv_dirs;
+    bool is_data_split = true;
+    SECTION("Original Benchmark")
+    {
+        is_data_split = false;
+
+#ifdef IPC_TOOLKIT_CCD_BENCHMARK_DIR
+        const fs::path root_path(IPC_TOOLKIT_CCD_BENCHMARK_DIR);
+#else
+        const fs::path root_path(CCD_IO_SAMPLE_QUERIES_DIR);
+#endif
+
+        const std::vector<std::string> folders { {
+            "chain",
+            "cow-heads",
+            "golf-ball",
+            "mat-twist",
+            "erleben-sliding-spike",
+            "erleben-spike-wedge",
+            "erleben-sliding-wedge",
+            "erleben-wedge-crack",
+            "erleben-spike-crack",
+            "erleben-wedges",
+            "erleben-cube-cliff-edges",
+            "erleben-spike-hole",
+            "erleben-cube-internal-edges",
+            "erleben-spikes",
+            "unit-tests",
+        } };
+
+        const std::array<std::string, 2> subfolders = { {
+            "edge-edge",
+            "vertex-face",
+        } };
+
+        for (const std::string& folder : folders) {
+            for (const std::string& subfolder : subfolders) {
+                csv_dirs.push_back(root_path / folder / subfolder);
+            }
+        }
+    }
+#ifdef IPC_TOOLKIT_CCD_NEW_BENCHMARK_DIR
+    SECTION("New Benchmark")
+    {
+        bool is_data_split = true;
+
+        const fs::path root_path(IPC_TOOLKIT_CCD_NEW_BENCHMARK_DIR);
+
+        const std::vector<std::string> folders { {
+            "armadillo-rollers",
+            "cloth-ball",
+            "cloth-funnel",
+            "n-body-simulation",
+            "rod-twist",
+        } };
+
+        for (const std::string& folder : folders) {
+            csv_dirs.push_back(root_path / folder / "queries");
+        }
+    }
+#endif
+
+    int i = 0;
+    double total_time = 0;
+    int num_fp = 0, num_fn = 0;
+    igl::Timer timer;
+
+    for (const auto& csv_dir : csv_dirs) {
+        for (const auto& csv : fs::directory_iterator(csv_dir)) {
+            std::vector<ccd_io::CCDQuery> queries;
+            if (is_data_split) {
+                queries = ccd_io::read_ccd_queries(
+                    csv.path().string(),
+                    csv.path().parent_path().parent_path() / "mma_bool"
+                        / (csv.path().stem().string() + "_mma_bool.json"));
+            } else {
+                queries = ccd_io::read_ccd_queries(csv.path().string());
+            }
+
+            bool is_edge_edge = is_data_split
+                ? csv.path().stem().string().find("ee") != std::string::npos
+                : csv.path().parent_path().filename().string() == "edge-edge";
+
+            for (int qi = 0; qi < queries.size(); qi++) {
+                INFO(fmt::format("{}(Q{})", csv.path().string(), qi));
+                std::cout << "\r                      "
+                          << "\rTesting query: " << ++i << std::flush;
+
+                Eigen::Map<const Matrix8x3> V(&queries[qi].vertices[0][0]);
+                const bool expected_result = queries[qi].ground_truth;
+
+                bool result;
+                double toi;
+                timer.start();
+                if (is_edge_edge) {
+                    result = edge_edge_ccd(
+                        V.row(0), V.row(1), V.row(2), V.row(3), V.row(4),
+                        V.row(5), V.row(6), V.row(7), toi);
+                } else {
+                    result = point_triangle_ccd(
+                        V.row(0), V.row(1), V.row(2), V.row(3), V.row(4),
+                        V.row(5), V.row(6), V.row(7), toi);
+                }
+                timer.stop();
+                total_time += timer.getElapsedTimeInMicroSec();
+
+                if (result != expected_result) {
+                    if (result) {
+                        num_fp++;
+                    } else {
+                        num_fn++;
+                    }
+                }
+
+                if (result < expected_result) {
+                    fmt::print("\n");
+                }
+                CHECK(result >= expected_result); // false positive is ok
+            }
+        }
+    }
+    fmt::print("\n");
+    fmt::print("Total time: {}s\n", total_time / 1e6);
+    fmt::print("Average time: {}μs\n", total_time / i);
+    fmt::print("False positives: {}\n", num_fp);
+    fmt::print("False negatives: {}\n", num_fn);
+}
+
+#ifdef IPC_TOOLKIT_WITH_CORRECT_CCD
+TEST_CASE(
+    "Run CCD Benchmark on TI CCD",
+    "[ccd][benchmark][3D][point-triangle][edge-edge][.]")
+{
+    fmt::print("Tight Inclusion CCD:\n");
+#else
+TEST_CASE(
+    "Run CCD Benchmark on FP CCD",
+    "[ccd][benchmark][3D][point-triangle][edge-edge][!mayfail][.]")
+{
+    fmt::print("Floating-Point CCD:\n");
+#endif
+
+    run_benchmark(
+        [](const Eigen::Vector3d& ea0_t0, const Eigen::Vector3d& ea1_t0,
+           const Eigen::Vector3d& eb0_t0, const Eigen::Vector3d& eb1_t0,
+           const Eigen::Vector3d& ea0_t1, const Eigen::Vector3d& ea1_t1,
+           const Eigen::Vector3d& eb0_t1, const Eigen::Vector3d& eb1_t1,
+           double& toi) -> bool {
+            return edge_edge_ccd(
+                ea0_t0, ea1_t0, eb0_t0, eb1_t0, ea0_t1, ea1_t1, eb0_t1, eb1_t1,
+                toi);
+        },
+        [](const Eigen::Vector3d& p_t0, const Eigen::Vector3d& t0_t0,
+           const Eigen::Vector3d& t1_t0, const Eigen::Vector3d& t2_t0,
+           const Eigen::Vector3d& p_t1, const Eigen::Vector3d& t0_t1,
+           const Eigen::Vector3d& t1_t1, const Eigen::Vector3d& t2_t1,
+           double& toi) -> bool {
+            return point_triangle_ccd(
+                p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1, toi);
+        });
+}
+
+TEST_CASE(
+    "Run CCD Benchmark on ACCD",
+    "[ccd][benchmark][3D][point-triangle][edge-edge][.]")
+{
+    fmt::print("Additive CCD:\n");
+
+    run_benchmark(
+        [](const Eigen::Vector3d& ea0_t0, const Eigen::Vector3d& ea1_t0,
+           const Eigen::Vector3d& eb0_t0, const Eigen::Vector3d& eb1_t0,
+           const Eigen::Vector3d& ea0_t1, const Eigen::Vector3d& ea1_t1,
+           const Eigen::Vector3d& eb0_t1, const Eigen::Vector3d& eb1_t1,
+           double& toi) -> bool {
+            return additive_ccd::edge_edge_ccd(
+                ea0_t0, ea1_t0, eb0_t0, eb1_t0, ea0_t1, ea1_t1, eb0_t1, eb1_t1,
+                toi);
+        },
+
+        [](const Eigen::Vector3d& p_t0, const Eigen::Vector3d& t0_t0,
+           const Eigen::Vector3d& t1_t0, const Eigen::Vector3d& t2_t0,
+           const Eigen::Vector3d& p_t1, const Eigen::Vector3d& t0_t1,
+           const Eigen::Vector3d& t1_t1, const Eigen::Vector3d& t2_t1,
+           double& toi) -> bool {
+            return additive_ccd::point_triangle_ccd(
+                p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1, toi);
+        });
+}
+
+// TEST_CASE("Failing Benchmark Cases", "[ccd]")
+// {
+//     using Matrix8x3 = Eigen::Matrix<double, 8, 3, Eigen::RowMajor>;
+
+//     const static std::vector<std::string> paths;
+//     const static std::vector<std::vector<int>> qids;
+
+//     for (int i = 0; i < paths.size(); i++) {
+//         const std::vector<ccd_io::CCDQuery> queries =
+//         ccd_io::read_ccd_queries(
+//             std::string(CCD_IO_SAMPLE_QUERIES_DIR) + paths[i]);
+//         for (const auto& qi : qids[i]) {
+//             Eigen::Map<const Matrix8x3> V(&queries[qi].vertices[0][0]);
+//             const bool expected_result = queries[qi].ground_truth;
+
+//             bool result;
+//             double toi;
+//             // if (subfolder == "edge-edge") {
+//             //     result = edge_edge_ccd(
+//             //         V.row(0), V.row(1), V.row(2), V.row(3), V.row(4),
+//             //         V.row(5), V.row(6), V.row(7), toi);
+//             // } else {
+//             result = additive_ccd::point_triangle_ccd(
+//                 V.row(0), V.row(1), V.row(2), V.row(3), V.row(4), V.row(5),
+//                 V.row(6), V.row(7), toi);
+//             // }
+//             CHECK(result >= expected_result); // false positive is ok
+//         }
+//     }
+// }
+#endif
\ No newline at end of file
diff --git a/tests/ccd/test_edge_edge_ccd.cpp b/tests/ccd/test_edge_edge_ccd.cpp
new file mode 100644
index 000000000..16fe7c867
--- /dev/null
+++ b/tests/ccd/test_edge_edge_ccd.cpp
@@ -0,0 +1,162 @@
+#include <catch2/catch_all.hpp>
+
+#include <ipc/config.hpp>
+#include <ipc/ccd/ccd.hpp>
+#include <ipc/ccd/additive_ccd.hpp>
+
+using namespace ipc;
+
+static const double EPSILON = std::numeric_limits<float>::epsilon();
+
+#ifdef IPC_TOOLKIT_WITH_CORRECT_CCD
+TEST_CASE("Edge-Edge CCD", "[ccd][3D][edge-edge]")
+#else
+TEST_CASE("Edge-Edge CCD", "[ccd][3D][edge-edge][!mayfail]")
+#endif
+{
+    Eigen::Vector3d ea0_t0, ea1_t0, eb0_t0, eb1_t0, ea0_t1, ea1_t1, eb0_t1,
+        eb1_t1;
+    bool is_collision_expected;
+    bool conservative_check = true;
+
+    double tol = DEFAULT_CCD_TOLERANCE;
+    long max_iter = DEFAULT_CCD_MAX_ITERATIONS;
+    double tmax = 1;
+
+#if !defined(WIN32) || defined(NDEBUG)
+    SECTION("General")
+    {
+        double uy = GENERATE(-1.0, 0.0, 1 - EPSILON, 1.0, 1 + EPSILON, 2.0);
+        double e1x = GENERATE(
+            -1 - EPSILON, -1, -1 + EPSILON, -0.5, 0, 0.5, 1 - EPSILON, 1,
+            1 + EPSILON);
+
+        ea0_t0 << -1, -1, 0;
+        ea1_t0 << 1, -1, 0;
+        eb0_t0 << e1x, 1, -1;
+        eb1_t0 << e1x, 1, 1;
+
+        Eigen::Vector3d u0, u1;
+        SECTION("moving")
+        {
+            u0 << 0, uy, 0;
+            u1 << 0, -uy, 0;
+            is_collision_expected = uy >= 1.0 && e1x >= -1 && e1x <= 1;
+        }
+        SECTION("fixed")
+        {
+            u0 << 0, 2 * uy, 0;
+            u1.setZero();
+            is_collision_expected = uy >= 2.0 && e1x >= -1 && e1x <= 1;
+        }
+
+        ea0_t1 = ea0_t0 + u0;
+        ea1_t1 = ea1_t0 + u0;
+        eb0_t1 = eb0_t0 + u1;
+        eb1_t1 = eb1_t0 + u1;
+    }
+    SECTION("Double root test case 1")
+    {
+        ea0_t0 << -3.0022200, 0.2362580, 0.0165247;
+        ea1_t0 << -3.2347850, 0.8312380, -0.1151003;
+        eb0_t0 << -3.0319900, 0.3148750, 0.0000000;
+        eb1_t0 << -2.8548800, 0.0900349, 0.0000000;
+        ea0_t1 << -2.8995600, 0.0345838, 0.0638580;
+        ea1_t1 << -3.1716930, 0.6104858, -0.0713340;
+        eb0_t1 = eb0_t0;
+        eb1_t1 = eb1_t0;
+
+        is_collision_expected = true;
+    }
+    SECTION("Double root test case 2")
+    {
+        ea0_t0 << 0, 0, 1;
+        ea1_t0 << 0, 1, 1;
+        eb0_t0 << 0.1, 0.2, 2;
+        eb1_t0 << 0.1, 0.2, -1;
+        ea0_t1 << 1, 1, 0;
+        ea1_t1 << 0, 0, 0;
+        eb0_t1 = eb0_t0;
+        eb1_t1 = eb1_t0;
+
+        const double t = GENERATE(0.5, 0.8, 0.88, 0.9, 1.0);
+        ea0_t1 = (ea0_t1 - ea0_t0) * t + ea0_t0;
+        ea1_t1 = (ea1_t1 - ea1_t0) * t + ea1_t0;
+
+        bool is_collision_expected = true;
+    }
+#endif
+    SECTION("Slow Case 1")
+    {
+        ea0_t0 << 1, 0.50803125, 2.10835646075301e-18;
+        ea1_t0 << -2.38233935445388e-18, 0.50803125, 1;
+        eb0_t0 << -4.99999999958867e-07, 0.5, 0;
+        eb1_t0 << -4.99999999958867e-07, 0.5, 1;
+        ea0_t1 << 1, 0.47124375, 4.11078309465837e-18;
+        ea1_t1 << -2.8526707189104e-18, 0.47124375, 1;
+        eb0_t1 << -4.99999999958867e-07, 0.5, 0;
+        eb1_t1 << -4.99999999958867e-07, 0.5, 1;
+
+        tol = GENERATE(1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1);
+
+        is_collision_expected = true;
+        conservative_check = false;
+    }
+    SECTION("Slow Case 2")
+    {
+        ea0_t0 << 1.00002232466453, 0.500004786049044, -2.06727783590977e-05;
+        ea1_t0 << 1.64687846177844e-05, 0.499996645067319, 1.63939999009028e-05;
+        eb0_t0 << 1, 0.5, 0;
+        eb1_t0 << 0, 0.5, 0;
+        ea0_t1 << 1.00294282700155, 0.498652627047143, 0.003626320742036;
+        ea1_t1 << -0.00219276550735626, 0.500871179186644, -0.00315828804921928;
+        eb0_t1 << 1, 0.5, 0;
+        eb1_t1 << 0, 0.5, 0;
+
+        max_iter = 1e6;
+        tmax = 2.8076171875e-03;
+
+        is_collision_expected = true;
+        conservative_check = false;
+    }
+    SECTION("Adversarial ACCD Case")
+    {
+        const double d0 = GENERATE(1e-2, 1e-4, 1e-6, 1e-8);
+        ea0_t0 << -1.0, -d0 / 2, 0.0;
+        ea1_t0 << 1.0, -d0 / 2, 0.0;
+        eb0_t0 << 0.0, d0 / 2, -1.0;
+        eb1_t0 << 0.0, d0 / 2, 1.0;
+
+        const double dy = d0;
+        const double scale = 1e-3;
+
+        ea0_t1 = ea0_t0 + Eigen::Vector3d(-scale, dy, -scale);
+        ea1_t1 = ea1_t0 + Eigen::Vector3d(scale, dy, scale);
+        eb0_t1 = eb0_t0 + Eigen::Vector3d(scale, -dy, scale);
+        eb1_t1 = eb1_t0 + Eigen::Vector3d(-scale, -dy, -scale);
+
+        // this ternary operator is to force MSVC to use 1 or 0
+        is_collision_expected = dy >= d0 / 2 ? true : false;
+        conservative_check = false;
+    }
+    CAPTURE(int(is_collision_expected));
+
+    double toi;
+    bool is_colliding = edge_edge_ccd(
+        ea0_t0, ea1_t0, eb0_t0, eb1_t0, ea0_t1, ea1_t1, eb0_t1, eb1_t1, toi,
+        /*min_distance=*/0.0, tmax, tol, max_iter);
+    if (conservative_check) {
+        CHECK(is_colliding >= is_collision_expected);
+    } else {
+        CHECK(is_colliding == is_collision_expected);
+    }
+
+    is_colliding = additive_ccd::edge_edge_ccd(
+        ea0_t0, ea1_t0, eb0_t0, eb1_t0, ea0_t1, ea1_t1, eb0_t1, eb1_t1, toi,
+        /*min_distance=*/0.0, tmax);
+    if (conservative_check) {
+        CHECK(is_colliding >= is_collision_expected);
+    } else {
+        CHECK(is_colliding == is_collision_expected);
+    }
+}
diff --git a/tests/ccd/test_point_edge_ccd.cpp b/tests/ccd/test_point_edge_ccd.cpp
new file mode 100644
index 000000000..3588de6dd
--- /dev/null
+++ b/tests/ccd/test_point_edge_ccd.cpp
@@ -0,0 +1,266 @@
+#include <catch2/catch_all.hpp>
+
+#include <ipc/config.hpp>
+#include <ipc/ccd/ccd.hpp>
+#include <ipc/ccd/additive_ccd.hpp>
+
+#include "collision_generator.hpp"
+
+using namespace ipc;
+
+namespace {
+/// Compares the time of impact of different implementations
+/// against the expected time of impact
+void check_toi(
+    const Eigen::Vector2d& p_t0,
+    const Eigen::Vector2d& e0_t0,
+    const Eigen::Vector2d& e1_t0,
+    const Eigen::Vector2d& p_t1,
+    const Eigen::Vector2d& e0_t1,
+    const Eigen::Vector2d& e1_t1,
+    const double toi_expected)
+{
+    double toi;
+    bool is_colliding = point_edge_ccd(
+        p_t0, e0_t0, e1_t0, p_t1, e0_t1, e1_t1, toi,
+        /*min_distance=*/0.0, /*tmax=*/1.0, DEFAULT_CCD_TOLERANCE,
+        DEFAULT_CCD_MAX_ITERATIONS, /*conservative_rescaling=*/1.0);
+    CHECK(is_colliding);
+    CHECK(toi <= toi_expected);
+
+    is_colliding = additive_ccd::point_edge_ccd(
+        p_t0, e0_t0, e1_t0, p_t1, e0_t1, e1_t1, toi,
+        /*min_distance=*/0.0, /*tmax=*/1.0, /*conservative_rescaling=*/0.99);
+    CHECK(is_colliding);
+    CHECK(toi <= toi_expected);
+}
+} // namespace
+
+TEST_CASE("Point-edge 2D CCD", "[ccd][2D][point-edge]")
+{
+    Eigen::Vector2d p_t0, p_t1, e0_t0, e0_t1, e1_t0, e1_t1;
+    double toi_expected;
+    bool is_collision_expected;
+    SECTION("Edge becomes degenerate")
+    {
+        p_t0 << 0, 1;
+        e0_t0 << -1, 0;
+        e1_t0 << 1, 0;
+        SECTION("Edge degenerates before impact")
+        {
+            p_t1 << 0, -1;
+            e0_t1 << 3, 0;
+            e1_t1 << -3, 0;
+            // The edge will become degenerate at t=0.25
+        }
+        SECTION("Edge degenerates after impact")
+        {
+            p_t1 << 0, -1;
+            e0_t1 << 0.5, 0;
+            e1_t1 << -0.5, 0;
+            // The edge will become degenerate at t=2/3
+        }
+        // The point will collide with the edge at t=0.5
+        is_collision_expected = true;
+        toi_expected = 0.5;
+    }
+    SECTION("Edge moving right; point moving left")
+    {
+        p_t0 << -1, 0;
+        e0_t0 << 1, -1;
+        e1_t0 << 1, 1;
+
+        p_t1 << 1, 0;
+        e0_t1 << -1, -1;
+        e1_t1 << -1, 1;
+
+        is_collision_expected = true;
+        toi_expected = 0.5;
+    }
+    SECTION("Point on edge's line moving towards edge")
+    {
+        p_t0 << 0, 0;
+        e0_t0 << 0, 1;
+        e1_t0 << 0, 2;
+
+        p_t1 << 0, 2;
+        e0_t1 << 0, 1;
+        e1_t1 << 0, 2;
+
+        is_collision_expected = true;
+        toi_expected = 0.5;
+    }
+    SECTION("Point and edge moving parallel")
+    {
+        p_t0 << 0, 1;
+        e0_t0 << 1, 0;
+        e1_t0 << 1, 2;
+
+        p_t1 << 0, 2;
+        e0_t1 << 1, 1;
+        e1_t1 << 1, 3;
+
+        is_collision_expected = false;
+    }
+    SECTION("Point moving right; edge stretching vertically")
+    {
+        e1_t0 << 1, -1;
+        e1_t1 << 1, -2;
+        SECTION("Swap vertices order e_0 = [0, 2]")
+        {
+            p_t0 << 1, 1;
+            e0_t0 << 0, 0;
+
+            p_t1 << 1, 2;
+            e0_t1 << 1, 0;
+
+            std::swap(p_t0, e0_t0);
+            std::swap(p_t1, e0_t1);
+        }
+        SECTION("Swap vertices order e_0 = [1, 2]")
+        {
+            p_t0 << 0, 0;
+            e0_t0 << 1, 1;
+
+            p_t1 << 1, 0;
+            e0_t1 << 1, 2;
+        }
+        is_collision_expected = true;
+        toi_expected = 1.0;
+    }
+    SECTION("Point-point")
+    {
+        p_t0 << 1.11111, 0.5;
+        p_t1 << 0.888889, 0.5;
+        e0_t0 << 1, 0.5;
+        e1_t0 << 1, 0.75;
+        e0_t1 = e0_t0;
+        e1_t1 = e1_t0;
+
+        is_collision_expected = true;
+        toi_expected = 0.5;
+    }
+
+    CAPTURE(
+        p_t0.transpose(), e0_t0.transpose(), e1_t0.transpose(),
+        p_t1.transpose(), e0_t1.transpose(), e1_t1.transpose());
+
+    double toi;
+    bool is_colliding = point_edge_ccd(
+        p_t0, e0_t0, e1_t0, p_t1, e0_t1, e1_t1, toi,
+        /*min_distance=*/0.0,
+        /*tmax=*/1.0, DEFAULT_CCD_TOLERANCE, DEFAULT_CCD_MAX_ITERATIONS,
+        /*conservative_rescaling=*/1.0);
+    REQUIRE(is_colliding == is_collision_expected);
+    if (is_collision_expected) {
+        CHECK(toi <= toi_expected);
+    }
+
+    is_colliding = additive_ccd::point_edge_ccd(
+        p_t0, e0_t0, e1_t0, p_t1, e0_t1, e1_t1, toi,
+        /*min_distance=*/0.0, /*tmax=*/1.0, /*conservative_rescaling=*/0.99);
+    REQUIRE(is_colliding == is_collision_expected);
+    if (is_collision_expected) {
+        CHECK(toi <= toi_expected);
+    }
+}
+
+#ifdef IPC_TOOLKIT_WITH_CORRECT_CCD
+TEST_CASE("Point-edge 2D ToI", "[ccd][2D][point-edge][toi]")
+#else
+TEST_CASE("Point-edge 2D ToI", "[ccd][2D][point-edge][toi][!mayfail]")
+#endif
+{
+    Eigen::Vector2d p_t0, e0_t0, e1_t0;
+    Eigen::Vector2d dp, de0, de1;
+
+    double expected_toi = -1;
+
+    SECTION("perpendicular") // (alpha=0.5)
+    {
+        p_t0 << 0.0, 1.0;
+        e0_t0 << -1.0, 0.0;
+        e1_t0 << 1.0, 0.0;
+
+        // touches, intersects, passes-trough
+        auto vel = GENERATE(1.0 + 1e-6, 2.0, 4.0);
+        // moving direction:
+        //  -> ->, -> -, -> <-, - <-, <- <-
+        auto j = GENERATE(0, 1, 2, 3, 4);
+        // extension, no-deform, compression,
+        auto dx = GENERATE(0.5, 0.0, -0.5);
+
+        dp << 0.0, -(3 - j) * vel / 2.0;
+        de0 << -dx, (j - 1.0) * vel / 2.0;
+        de1 << dx, (j - 1.0) * vel / 2.0;
+
+        expected_toi = 1.0 / vel;
+    }
+
+    SECTION("double impact") // (rotating edge)
+    {
+        e0_t0 << -1.0, 0.0;
+        e1_t0 << 1.0, 0.0;
+        p_t0 << 0.0, 0.5;
+
+        de0 << 1.6730970740318298, 0.8025388419628143;
+        de1 << -1.616142749786377, -0.6420311331748962;
+        dp << 0.0, -1.0;
+
+        expected_toi = 0.4482900963;
+    }
+
+    SECTION("random")
+    {
+        // This times out on Windows debug
+#if defined(NDEBUG) || !(defined(WIN32) || defined(_WIN32) || defined(__WIN32))
+        auto impact = GENERATE(random_impacts(100, /*rigid=*/true));
+#else
+        auto impact = GENERATE(random_impacts(10, /*rigid=*/true));
+#endif
+
+        p_t0 = impact.p_t0;
+        e0_t0 = impact.e0_t0;
+        e1_t0 = impact.e1_t0;
+
+        dp = impact.dp;
+        de0 = impact.de0;
+        de1 = impact.de1;
+
+        expected_toi = impact.toi;
+    }
+
+#ifdef NDEBUG           // This case takes forever to run
+    SECTION("parallel") //  (alpha=0 || alpha = 1)
+    {
+        p_t0 << 0.5, 0.0;
+        e0_t0 << -0.5, 0.0;
+        e1_t0 << -1.5, 0.0;
+
+        // touches, intersects, passes-trough
+        // double vel = GENERATE(1.0 + 1e-6, 2.0, 4.0);
+        double vel = 2.0;
+        // moving: both (same), ij, both (op), kl, both (same)
+        // double j = GENERATE(0, 1, 2, 3, 4);
+        double j = 0;
+        // extension, no-deform, compression,
+        // double dy = GENERATE(0.5, 0.0, -0.5);
+        double dy = 0.0;
+
+        dp << -(3 - j) * vel / 2.0, 0.0;
+        de0 << (j - 1.0) * vel / 2.0, 0.0;
+        // we only move one so we don't change the toi
+        de1 << (j - 1.0) * vel / 2.0, dy;
+
+        expected_toi = 1.0 / vel;
+    }
+#endif
+
+    Eigen::Vector2d p_t1 = p_t0 + dp;
+    Eigen::Vector2d e0_t1 = e0_t0 + de0;
+    Eigen::Vector2d e1_t1 = e1_t0 + de1;
+
+    check_toi(p_t0, e0_t0, e1_t0, p_t1, e0_t1, e1_t1, expected_toi);
+    // Flip the edges
+    check_toi(p_t0, e1_t0, e0_t0, p_t1, e1_t1, e0_t1, expected_toi);
+}
\ No newline at end of file
diff --git a/tests/ccd/test_point_triangle_ccd.cpp b/tests/ccd/test_point_triangle_ccd.cpp
new file mode 100644
index 000000000..c7fee8d4a
--- /dev/null
+++ b/tests/ccd/test_point_triangle_ccd.cpp
@@ -0,0 +1,112 @@
+#include <catch2/catch_all.hpp>
+
+#include <ipc/config.hpp>
+#include <ipc/ccd/ccd.hpp>
+#include <ipc/ccd/additive_ccd.hpp>
+
+using namespace ipc;
+
+static const double EPSILON = std::numeric_limits<float>::epsilon();
+
+#ifdef IPC_TOOLKIT_WITH_CORRECT_CCD
+TEST_CASE("Point-Triangle CCD", "[ccd][3D][point-triangle]")
+#else
+TEST_CASE("Point-Triangle CCD", "[ccd][3D][point-triangle][!mayfail]")
+#endif
+{
+    Eigen::Vector3d p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1;
+    bool is_collision_expected;
+    bool conservative_check = true;
+
+    std::string name;
+    SECTION("General")
+    {
+        name = "General";
+        const double v0z = GENERATE(0.0, -1.0);
+        const double u0y =
+            -GENERATE(-1.0, 0.0, 0.5 - EPSILON, 0.5, 0.5 + EPSILON, 1.0, 2.0);
+        const double u0z = GENERATE(-EPSILON, 0.0, EPSILON);
+        const double u1y =
+            GENERATE(-1.0, 0.0, 0.5 - EPSILON, 0.5, 0.5 + EPSILON, 1.0, 2.0);
+        CAPTURE(v0z, u0y, u1y, u0z);
+
+        p_t0 << 0, 1, v0z;
+        t0_t0 << -1, 0, 1;
+        t1_t0 << 1, 0, 1;
+        t2_t0 << 0, 0, -1;
+
+        p_t1 = p_t0 + Eigen::Vector3d(0, u0y, u0z);
+        t0_t1 = t0_t0 + Eigen::Vector3d(0, u1y, 0);
+        t1_t1 = t1_t0 + Eigen::Vector3d(0, u1y, 0);
+        t2_t1 = t2_t0 + Eigen::Vector3d(0, u1y, 0);
+
+        is_collision_expected = ((-u0y + u1y >= 1) && (v0z + u0z >= t2_t0.z()));
+    }
+#ifdef NDEBUG
+#ifndef WIN32
+    SECTION("Zhongshi's test case")
+    {
+        name = "Zhongshi's test case";
+        const double qy = GENERATE(-EPSILON, 0, EPSILON);
+        CAPTURE(qy);
+
+        p_t0 << 0, qy, 0;
+        t0_t0 << 0, 0, 0;
+        t1_t0 << 0, 1, 0;
+        t2_t0 << 1, 0, 0;
+        p_t1 << 0, qy, 0;
+        t0_t1 << 0, 0, 1;
+        t1_t1 << 0, 1, 1;
+        t2_t1 << 1, 0, 1;
+
+        is_collision_expected = qy >= 0;
+    }
+#endif
+    SECTION("Bolun's test case")
+    {
+        name = "Bolun's test case";
+        p_t0 << 0.1, 0.1, 0.1;
+        t0_t0 << 0, 0, 1;
+        t1_t0 << 1, 0, 1;
+        t2_t0 << 0, 1, 1;
+        p_t1 << 0.1, 0.1, 0.1;
+        t0_t1 << 0, 0, 0;
+        t1_t1 << 0, 1, 0;
+        t2_t1 << 1, 0, 0;
+        is_collision_expected = true;
+    }
+#endif
+    SECTION("No Zero ToI CCD")
+    {
+        name = "No Zero ToI CCD";
+        p_t0 << 0.0133653, 0.100651, -0.0215935;
+        t0_t0 << 0.0100485, 0.0950896, -0.0171013;
+        t1_t0 << 0.0130388, 0.100666, -0.0218112;
+        t2_t0 << 0.015413, 0.100554, -0.0202265;
+        p_t1 << 0.0133652999767858, 0.099670000268615, -0.0215934999996444;
+        t0_t1 << 0.0100484999799995, 0.0941086002577558, -0.0171012999972189;
+        t1_t1 << 0.0130387999724314, 0.0996850002629403, -0.0218111999936902;
+        t2_t1 << 0.0154129999740718, 0.0995730002646605, -0.020226499996014;
+        is_collision_expected = false;
+        conservative_check = false;
+    }
+    INFO(name);
+
+    double toi;
+    bool is_colliding = point_triangle_ccd(
+        p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1, toi);
+
+    if (conservative_check) {
+        CHECK(is_colliding >= is_collision_expected);
+    } else {
+        CHECK(is_colliding == is_collision_expected);
+    }
+
+    is_colliding = additive_ccd::point_triangle_ccd(
+        p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1, toi);
+    if (conservative_check) {
+        CHECK(is_colliding >= is_collision_expected);
+    } else {
+        CHECK(is_colliding == is_collision_expected);
+    }
+}
diff --git a/tests/friction/test_force_jacobian.cpp b/tests/friction/test_force_jacobian.cpp
index af60e9c07..fc45d2f02 100644
--- a/tests/friction/test_force_jacobian.cpp
+++ b/tests/friction/test_force_jacobian.cpp
@@ -206,7 +206,7 @@ void check_friction_force_jacobian(
     CHECK(fd::compare_jacobian(-jac_force, hess_D));
 }
 
-TEST_CASE("Test friction force jacobian", "[friction][force-jacobian][thisone]")
+TEST_CASE("Friction force jacobian", "[friction][force-jacobian]")
 {
     const int x_case = GENERATE(0, 1);
     FrictionData data = friction_data_generator();
@@ -236,7 +236,7 @@ TEST_CASE("Test friction force jacobian", "[friction][force-jacobian][thisone]")
 }
 
 TEST_CASE(
-    "Test friction force jacobian on real data",
+    "Friction force jacobian on real data",
     "[friction][force-jacobian][real-data]")
 {
     bool use_convergent_formulation = GENERATE(true, false);
diff --git a/tests/friction/test_friction.cpp b/tests/friction/test_friction.cpp
index e99d65663..fb875f1ab 100644
--- a/tests/friction/test_friction.cpp
+++ b/tests/friction/test_friction.cpp
@@ -18,7 +18,7 @@
 
 using namespace ipc;
 
-TEST_CASE("Test friction gradient and hessian", "[friction][gradient][hessian]")
+TEST_CASE("Friction gradient and hessian", "[friction][gradient][hessian]")
 {
     FrictionData data = friction_data_generator();
     const auto& [V0, V1, E, F, collision_constraints, mu, epsv_times_h, dhat, barrier_stiffness] =
diff --git a/tests/test_cfl.cpp b/tests/test_cfl.cpp
index 7674c8aac..d3c070483 100644
--- a/tests/test_cfl.cpp
+++ b/tests/test_cfl.cpp
@@ -8,7 +8,7 @@
 
 using namespace ipc;
 
-TEST_CASE("Test compute_cfl_stepsize()", "[ccd][cfl]")
+TEST_CASE("Compute CFL stepsize", "[ccd][cfl]")
 {
     const double dhat = GENERATE(1e-6, 1e-3, 1e-1);
 
diff --git a/tests/test_has_intersections.cpp b/tests/test_has_intersections.cpp
index 815d0917c..fc9ea0372 100644
--- a/tests/test_has_intersections.cpp
+++ b/tests/test_has_intersections.cpp
@@ -80,7 +80,7 @@ bool combine_meshes(
     return true;
 }
 
-// TEST_CASE("Test HashGrid EF vs brute-force", "[intersection][brute_force]")
+// TEST_CASE("HashGrid EF vs brute-force", "[intersection][brute_force]")
 // {
 // #ifdef NDEBUG
 //     std::string mesh1_name = GENERATE("cube.obj", "bunny.obj");
@@ -177,7 +177,7 @@ bool combine_meshes(
 //     CAPTURE(mesh1_name, mesh2_name, R1, R2);
 // }
 
-TEST_CASE("Test has_intersections()", "[intersection][thisone]")
+TEST_CASE("Has intersections", "[intersection]")
 {
     std::string mesh1_name = GENERATE("cube.obj", "bunny.obj");
     std::string mesh2_name = GENERATE("cube.obj", "bunny.obj");
diff --git a/tests/test_ipc.cpp b/tests/test_ipc.cpp
index 1828c565a..37bf24887 100644
--- a/tests/test_ipc.cpp
+++ b/tests/test_ipc.cpp
@@ -66,7 +66,7 @@ TEST_CASE("Dummy test for IPC compilation", "[!benchmark][ipc]")
     };
 }
 
-TEST_CASE("Test IPC full gradient", "[ipc][gradient]")
+TEST_CASE("IPC full gradient", "[ipc][gradient]")
 {
     const BroadPhaseMethod method = GENERATE_BROAD_PHASE_METHODS();
     const bool use_convergent_formulation = GENERATE(true, false);
@@ -134,7 +134,7 @@ TEST_CASE("Test IPC full gradient", "[ipc][gradient]")
     CHECK(fd::compare_gradient(grad_b, fgrad_b));
 }
 
-TEST_CASE("Test IPC full hessian", "[ipc][hessian]")
+TEST_CASE("IPC full hessian", "[ipc][hessian]")
 {
     const BroadPhaseMethod method = GENERATE_BROAD_PHASE_METHODS();
     const bool use_convergent_formulation = GENERATE(true, false);
@@ -202,7 +202,7 @@ TEST_CASE("Test IPC full hessian", "[ipc][hessian]")
     CHECK(fd::compare_hessian(hess_b, fhess_b, 1e-3));
 }
 
-TEST_CASE("Test IPC shape derivative", "[ipc][shape_opt]")
+TEST_CASE("IPC shape derivative", "[ipc][shape_opt]")
 {
     nlohmann::json data;
     {
@@ -312,7 +312,7 @@ TEST_CASE("Benchmark IPC shape derivative", "[ipc][shape_opt][!benchmark]")
     };
 }
 
-TEST_CASE("Test convergent formulation", "[ipc][convergent]")
+TEST_CASE("Convergent formulation", "[ipc][convergent]")
 {
     const bool use_convergent_formulation = GENERATE(false, true);
     const double dhat = 1e-3;