Offset Geometric Contact

.gitignore

@@ -670,4 +670,5 @@ python/update_bindings.py
 
 # External test data
 tests/data
+notebooks/*.png
 CMakeGraphVizOptions.cmake

CMakeLists.txt

@@ -192,13 +192,13 @@ endif()
 # Faster unordered map
 if(IPC_TOOLKIT_WITH_ROBIN_MAP)
   include(robin_map)
-  target_link_libraries(ipc_toolkit PUBLIC tsl::robin_map)
+  target_link_libraries(ipc_toolkit PRIVATE tsl::robin_map)
 endif()
 
 # Hashes
 if(IPC_TOOLKIT_WITH_ABSEIL)
   include(abseil)
-  target_link_libraries(ipc_toolkit PUBLIC absl::hash)
+  target_link_libraries(ipc_toolkit PRIVATE absl::hash)
 endif()
 
 # Intervals
@@ -297,7 +297,7 @@ endif()
 # Xcode
 ################################################################################
 
-if (CMAKE_GENERATOR STREQUAL "Xcode")
+if (IPC_TOOLKIT_TOPLEVEL_PROJECT AND CMAKE_GENERATOR STREQUAL "Xcode")
   set(CMAKE_XCODE_SCHEME_LAUNCH_CONFIGURATION "${CMAKE_BUILD_TYPE}")
   set_target_properties(ipc_toolkit PROPERTIES XCODE_GENERATE_SCHEME ON)
   if(IPC_TOOLKIT_BUILD_TESTS)

cmake/ipc_toolkit/ipc_toolkit_tests_data.cmake

@@ -1,4 +1,4 @@
-# ipc-toolkit-tests-data (https://github.com/ipc-sim/ipc-toolkit-tests-data)
+# ipc-toolkit-test-data (https://github.com/ipc-sim/ipc-toolkit-test-data)
 # License: MIT
 
 if(TARGET ipc_toolkit_test_data_download)
@@ -29,7 +29,7 @@ else()
         PREFIX "${FETCHCONTENT_BASE_DIR}/tests/data"
         SOURCE_DIR ${IPC_TOOLKIT_TESTS_DATA_DIR}
 
-        GIT_REPOSITORY https://github.com/ipc-sim/ipc-toolkit-tests-data.git
+        GIT_REPOSITORY https://github.com/ipc-sim/ipc-toolkit-test-data.git
         GIT_TAG 0c6ce752d020db32f62da37d7b8ca9f73c204b36
 
         CONFIGURE_COMMAND ""

docs/source/_static/graphviz/dependencies.dot

@@ -47,7 +47,7 @@ digraph "IPC Toolkit Dependencies" {
     "node5" -> "node3" [color = "#BE6562";];
     // ipc_toolkit -> igl_predicates
     "node6" [label = "robin_map\n(tsl::robin_map)";shape = box;style = "rounded,filled";fillcolor = "#FFE6CC";color = "#DAA52D";];
-    "node5" -> "node6" [color = "#8FB976";];
+    "node5" -> "node6" [color = "#BE6562";];
     // ipc_toolkit -> robin_map
     "node7" [label = "scalable_ccd\n(scalable_ccd::scalable_ccd)";shape = box;style = "rounded,filled";fillcolor = "#D5E8D4";color = "#8FB976";];
     "node7" -> "node0" [color = "#8FB976";];
@@ -77,5 +77,5 @@ digraph "IPC Toolkit Dependencies" {
     "node5" -> "node11" [color = "#BE6562";];
     // ipc_toolkit -> tight_inclusion
     "node12" [label = "absl_hash\n(absl::hash)";shape = box;style = "rounded,filled";fillcolor = "#D5E8D4";color = "#8FB976";];
-    "node5" -> "node12" [color = "#8FB976";];
+    "node5" -> "node12" [color = "#BE6562";];
 }

docs/source/cpp-api/barrier/TwoStageBarrier.rst

@@ -0,0 +1,5 @@
+Two-Stage Barrier
+=================
+
+.. doxygenclass:: ipc::TwoStageBarrier
+    :allow-dot-graphs:

docs/source/cpp-api/barrier/index.rst

@@ -10,6 +10,7 @@ Barrier
    NormalizedBarrier.rst
    ClampedLogSqBarrier.rst
    CubicBarrier.rst
+   TwoStageBarrier.rst
 
 Barrier functions and functionals.
 
@@ -21,16 +22,19 @@ Types
 
     * - Name
       - Description
-    * - :doc:`Barrier <Barrier>`
+    * - :cpp:class:`ipc::Barrier`
       - Base class for barrier functions.
-    * - :doc:`ClampedLogBarrier <ClampedLogBarrier>`
+    * - :cpp:class:`ipc::ClampedLogBarrier`
       - Smoothly clamped log barrier functions from [Li et al.].
-    * - :doc:`ClampedLogSqBarrier <ClampedLogSqBarrier>`
+    * - :cpp:class:`ipc::ClampedLogSqBarrier`
       - Clamped log barrier with a quadratic log term from [Huang et al.].
-    * - :doc:`CubicBarrier <CubicBarrier>`
+    * - :cpp:class:`ipc::CubicBarrier`
       - Cubic barrier function from [Ando 2024].
-    * - :doc:`NormalizedBarrier <NormalizedBarrier>`
+    * - :cpp:class:`ipc::NormalizedBarrier`
       - Normalized barrier function from [Li et al.].
+    * - :cpp:class:`ipc::TwoStageBarrier`
+      - Two-stage barrier function from [Chen et al. 2025].
+
 
 Functions
 ---------
@@ -40,21 +44,21 @@ Functions
 
     * - Name
       - Description
-    * - :func:`ipc::barrier`
+    * - :cpp:func:`ipc::barrier`
       - Evaluate the barrier function.
-    * - :func:`ipc::barrier_first_derivative`
+    * - :cpp:func:`ipc::barrier_first_derivative`
       - Derivative of the barrier function.
-    * - :func:`ipc::barrier_second_derivative`
+    * - :cpp:func:`ipc::barrier_second_derivative`
       - Second derivative of the barrier function.
-    * - :func:`ipc::barrier_force_magnitude`
+    * - :cpp:func:`ipc::barrier_force_magnitude`
       - Compute the barrier force magnitude.
-    * - :func:`ipc::barrier_force_magnitude_gradient`
+    * - :cpp:func:`ipc::barrier_force_magnitude_gradient`
       - Compute the gradient of the barrier force magnitude.
-    * - :func:`ipc::initial_barrier_stiffness`
+    * - :cpp:func:`ipc::initial_barrier_stiffness`
       - Compute the initial barrier stiffness.
-    * - :func:`ipc::update_barrier_stiffness`
+    * - :cpp:func:`ipc::update_barrier_stiffness`
       - Update the barrier stiffness based on the current state.
-    * - :func:`ipc::semi_implicit_stiffness`
+    * - :cpp:func:`ipc::semi_implicit_stiffness`
       - Compute the semi-implicit stiffness for all collisions.
 
 Function Details

notebooks/ogc.py

@@ -0,0 +1,194 @@
+import pathlib
+from find_ipctk import ipctk
+import polyscope as ps
+from polyscope import imgui
+import meshio
+import numpy as np
+
+import argparse
+
+parser = argparse.ArgumentParser()
+parser.add_argument(
+    "mesh",
+    type=pathlib.Path,
+    help="Path to the mesh file (e.g., .ply)"
+)
+args = parser.parse_args()
+
+mesh = meshio.read(args.mesh)
+
+# Rotate 90 degrees around the x-axis
+R = np.array([[1, 0, 0],
+              [0, 0, -1],
+              [0, 1, 0]])
+
+cmesh = ipctk.CollisionMesh(
+    # mesh.points @ R.T,
+    mesh.points,
+    ipctk.edges(mesh.cells_dict["triangle"]),
+    mesh.cells_dict["triangle"],
+)
+
+min_edge_length = np.min(
+    np.linalg.norm(cmesh.rest_positions[cmesh.edges[:, 0]] -
+                   cmesh.rest_positions[cmesh.edges[:, 1]], axis=1))
+
+max_dhat = 2e-1
+dhat = 2e-3  # 0.5 * min_edge_length
+use_ogc = True
+use_area_weighting = True
+
+# Prebuild the candidates since dhat is the max it will ever be
+candidates = ipctk.Candidates()
+candidates.build(cmesh, cmesh.rest_positions,
+                 max_dhat, broad_phase=ipctk.BVH())
+
+collisions = ipctk.NormalCollisions()
+
+
+def contact_force():
+    global collisions
+    collisions.clear()
+    collisions.collision_set_type = (
+        ipctk.NormalCollisions.CollisionSetType.OGC if use_ogc else ipctk.NormalCollisions.CollisionSetType.IPC)
+    collisions.use_area_weighting = use_area_weighting
+    collisions.build(candidates, cmesh, cmesh.rest_positions, dhat)
+
+    B = ipctk.BarrierPotential(dhat)
+    f = -B.gradient(collisions, cmesh,
+                    cmesh.rest_positions).reshape(-1, 3, order="C")
+    f *= 1e5
+
+    return f
+
+
+ps.init()
+ps.set_give_focus_on_show(True)
+ps.set_ground_plane_mode("shadow_only")
+
+ps_mesh = ps.register_surface_mesh(
+    args.mesh.stem,
+    cmesh.rest_positions,
+    cmesh.faces,
+    edge_width=1.0,
+)
+
+f = contact_force()
+ps_force = ps_mesh.add_scalar_quantity(
+    "force",
+    np.linalg.norm(f, axis=1),
+    defined_on="vertices",
+    cmap="reds",
+    vminmax=(0, 1e-1),
+    enabled=True,
+)
+
+ps_force_vector = ps_mesh.add_vector_quantity(
+    "force_vector",
+    f,
+    defined_on="vertices",
+    radius=0.01,
+    enabled=True,
+    # vectortype="ambient",
+)
+
+# xi = 0
+ev_i = 0
+ps_ev = ps.register_curve_network(
+    "edge-vertex",
+    cmesh.rest_positions[0:3],  # Temporarily use first 3 vertices
+    np.array([[1, 2]]),
+    enabled=False,
+)
+ee_i = 0
+ps_ee = ps.register_curve_network(
+    "edge-edge",
+    cmesh.rest_positions[0:4],  # Temporarily use first 4 vertices
+    np.array([[0, 1], [2, 3]]),
+    enabled=True,
+)
+ps_ee.set_radius(dhat/2, relative=False)
+ps_ee_closest_points = ps.register_point_cloud(
+    "edge-edge-closest-points",
+    cmesh.rest_positions[0:2],  # Temporarily use first 2 vertices
+    enabled=True,
+)
+
+
+def callback():
+    global dhat, use_ogc, ev_i, ee_i, max_dhat
+    changed0, dhat = imgui.SliderFloat("dhat", dhat, 1e-3, max_dhat, "%.5f")
+    if changed0 and dhat > max_dhat:
+        max_dhat = dhat
+        candidates.build(cmesh, cmesh.rest_positions, max_dhat,
+                         broad_phase=ipctk.BVH())
+    changed1, use_ogc = imgui.Checkbox("use_ogc", use_ogc)
+    if changed0 or changed1:
+        f = contact_force()
+        ps_mesh.add_scalar_quantity("force", np.linalg.norm(f, axis=1))
+        ps_mesh.add_vector_quantity("force_vector", f)
+    imgui.Text(
+        f"Number of collisions: {len(collisions.vv_collisions)} VV, {len(collisions.ev_collisions)} EV,"
+        f" {len(collisions.ee_collisions)} EE, {len(collisions.fv_collisions)} FV")
+
+    if len(collisions.ev_collisions) > 0:
+        ev_i = min(ev_i, len(collisions.ev_collisions) - 1)
+        _, ev_i = imgui.SliderInt("Edge-Vertex Index", ev_i,
+                                  0, len(collisions.ev_collisions) - 1)
+        ps_ev.update_node_positions(cmesh.rest_positions[
+            collisions.ev_collisions[ev_i].vertex_ids(cmesh.edges, cmesh.faces)[:3]])
+    else:
+        ps_ev.set_enabled(False)
+
+    if len(collisions.ee_collisions) > 0:
+        ee_i = min(ee_i, len(collisions.ee_collisions) - 1)
+        changed2, ee_i = imgui.SliderInt("Edge-Edge Index", ee_i,
+                                         0, len(collisions.ee_collisions) - 1)
+        coeffs = collisions.ee_collisions[ee_i].compute_coefficients(
+            cmesh.rest_positions, cmesh.edges, cmesh.faces)
+        imgui.Text("Coeffs: {}".format(coeffs))
+        if changed0 or changed1 or changed2:
+            V_ee = cmesh.rest_positions[
+                collisions.ee_collisions[ee_i].vertex_ids(cmesh.edges, cmesh.faces)]
+            ps_ee.update_node_positions(V_ee)
+            if changed0:
+                ps_ee.set_radius(dhat/2, relative=False)
+            ps_ee_closest_points.update_point_positions(
+                np.vstack([
+                    coeffs[0] * V_ee[0] + coeffs[1] * V_ee[1],
+                    -coeffs[2] * V_ee[2] - coeffs[3] * V_ee[3]]))
+    else:
+        ps_ee.set_enabled(False)
+        ps_ee_closest_points.set_enabled(False)
+
+    # if changed2:
+    #     adj_vi = cmesh.edge_vertex_adjacencies[ei]
+    #     if len(adj_vi) == 1:
+    #         adj_vi = [adj_vi[0], adj_vi[0]]
+    #     ps_edge.update_node_positions(
+    #         cmesh.rest_positions[[
+    #             *cmesh.edges[ei],
+    #             *adj_vi
+    #         ]])
+    #
+    # _, xi = imgui.SliderInt("Vertex Index##mine", xi, 0,
+    #                         len(cmesh.rest_positions) - 1)
+    # x = cmesh.rest_positions[xi]
+    # imgui.Text(f"Vertex position: {x}")
+    # x0, x1 = cmesh.rest_positions[cmesh.edges[ei]]
+    # passed = True
+    # for vi in cmesh.edge_vertex_adjacencies[ei]:
+    #     x2 = cmesh.rest_positions[vi]
+    #     alpha = ipctk.point_edge_closest_point(x2, x0, x1)
+    #     imgui.Text(f"alpha: {alpha:.5f}")
+    #     p = (x1 - x0) * alpha + x0
+    #     imgui.Text(f"p: {p}, x2: {x2}, x: {x}")
+    #     imgui.Text(f"b: {(x - p).dot(p - x2)}")
+    #     passed = passed and (x - p).dot(p - x2) > 0
+    # imgui.Text(f"check_edge_feasible_region: {passed}")
+
+
+ps.set_user_callback(callback)
+
+
+ps.show()

python/_find_ipctk.py

@@ -1,20 +1,32 @@
-try:
-    import ipctk  # Try to import the built module
-except ImportError:
-    import sys
-    import pathlib
+import sys
+import pathlib
+import importlib.util
+
+# First, try to find the system-installed version
+spec = importlib.util.find_spec("ipctk")
+if spec is not None:
+    ipctk = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(ipctk)
+    print("Using system-installed ipctk module")
+else:
     repo_root = pathlib.Path(__file__).parents[1]
     possible_paths = [
         pathlib.Path("python").resolve(),
         repo_root / "build" / "python",
         repo_root / "build" / "release" / "python",
         repo_root / "build" / "debug" / "python",
     ]
+
     for path in possible_paths:
-        if path.exists() and len(list(path.glob("ipctk.*"))) > 0:
+        if path.exists():
             sys.path.append(str(path))
-            break
+            spec = importlib.util.find_spec("ipctk")
+            if spec is not None:
+                print(f"Using found ipctk module at {path}")
+                ipctk = importlib.util.module_from_spec(spec)
+                spec.loader.exec_module(ipctk)
+                break
+            else:
+                sys.path.pop()  # Remove the path if ipctk is not found
     else:
         raise ImportError("Could not find the ipctk module")
-    print(f"Using found ipctk module at {path}")
-    import ipctk  # Try again