Extremum Sum Potential

CMakeLists.txt

@@ -100,6 +100,7 @@ option(IPC_TOOLKIT_WITH_ABSEIL                "Use Abseil's hash functions"
 option(IPC_TOOLKIT_WITH_FILIB                 "Use filib for interval arithmetic"              ON)
 option(IPC_TOOLKIT_WITH_INEXACT_CCD           "Use the original inexact CCD method of IPC"    OFF)
 option(IPC_TOOLKIT_WITH_PROFILER              "Enable performance profiler"                   OFF)
+option(IPC_TOOLKIT_WITH_GEOGRAM               "Use geogram for predicates / expansion types"  ON)
 option(IPC_TOOLKIT_WITH_TRACY                 "Enable Tracy frame profiler"                   OFF)
 option(IPC_TOOLKIT_WITH_MESHFEM_SPARSE        "Use MeshFEMSparse for block-accelerated assembly" ON)
 
@@ -302,6 +303,12 @@ if(IPC_TOOLKIT_WITH_FILIB)
   target_link_libraries(ipc_toolkit PUBLIC filib::filib)
 endif()
 
+# Geogram
+if(IPC_TOOLKIT_WITH_GEOGRAM)
+  include(geogram)
+  target_link_libraries(ipc_toolkit PRIVATE geogram::geogram)
+endif()
+
 # Block-accelerated Hessian assembly
 if(IPC_TOOLKIT_WITH_MESHFEM_SPARSE)
   include(meshfem_sparse)

cmake/find/FindSIMD.cmake

@@ -339,10 +339,10 @@ function(test_avx_availability)
 	elseif(CMAKE_CXX_COMPILER_ID MATCHES "Intel")
 		set(AVX_FLAGS "-xHost")
 	elseif(MSVC)
+		# /arch: is single-valued on MSVC: emitting both downgrades to AVX.
 		if(DETECTED_AVX_20)
 			SET(AVX_FLAGS "${AVX_FLAGS} /arch:AVX2")
-		endif()
-		if(DETECTED_AVX_10)
+		elseif(DETECTED_AVX_10)
 			SET(AVX_FLAGS "${AVX_FLAGS} /arch:AVX")
 		endif()
 	endif()

cmake/recipes/geogram.cmake

@@ -0,0 +1,24 @@
+# geogram (https://github.com/BrunoLevy/geogram)
+# License: BSD 3-Clause License
+if(TARGET geogram::geogram)
+    return()
+endif()
+
+message(STATUS "Third-party: creating target 'geogram::geogram'")
+
+include(CPM)
+CPMAddPackage(
+    URI "gh:BrunoLevy/[email protected]"
+    OPTIONS
+        "GEOGRAM_WITH_GRAPHICS OFF"
+        "GEOGRAM_WITH_LEGACY_NUMERICS OFF"
+        "GEOGRAM_WITH_HLBFGS OFF"
+        "GEOGRAM_WITH_TETGEN OFF"
+        "GEOGRAM_WITH_TRIANGLE OFF"
+        "GEOGRAM_WITH_LUA OFF"
+        "GEOGRAM_LIB_ONLY ON"
+)
+
+if(NOT TARGET geogram::geogram AND TARGET geogram)
+    add_library(geogram::geogram ALIAS geogram)
+endif()

docs/source/_static/graphviz/dependencies.dot

@@ -98,4 +98,7 @@ digraph "IPC Toolkit Dependencies" {
     // MeshFEMSparse -> Eigen3_Eigen
     "node17" -> "node9" [color = "#8FB976";];
     // MeshFEMSparse -> tbb
+    "node18" [label = "geogram\n(geogram::geogram)";shape = box;style = "rounded,filled";fillcolor = "#CCE7F8";color = "#6596B2";];
+    "node5" -> "node18" [color = "#BE6562"; style = "dashed";];
+    // ipc_toolkit -> geogram (IPC_TOOLKIT_WITH_GEOGRAM)
 }

docs/source/about/dependencies.rst

@@ -93,6 +93,12 @@ Additionally, IPC Toolkit may optionally use the following libraries:
       - `github.com/MeshFEM/MeshFEMSparse <https://github.com/MeshFEM/MeshFEMSparse>`_
       - |:white_check_mark:|
       - ``IPC_TOOLKIT_WITH_MESHFEM_SPARSE``
+    * - geogram
+      - Exact geometric predicates for the exact distance-type classification used by the ESP potential
+      - BSD-3-Clause
+      - `github.com/BrunoLevy/geogram <https://github.com/BrunoLevy/geogram>`_
+      - |:white_check_mark:|
+      - ``IPC_TOOLKIT_WITH_GEOGRAM``
     * - nlohmann/json
       - JSON parsing for profiler and tests
       - MIT
@@ -128,6 +134,9 @@ Some of these libraries are enabled by default, and some are not. You can enable
 .. note::
     ``MeshFEMSparse`` (and its transitive dependency ``MeshFEMCore``) is downloaded source-only and compiled into a minimal static library (matrix data structures and assembly routines; no sparse direct solvers). When enabled (the default), :cpp:func:`ipc::Potential::hessian` assembles through the block-CSC backend — several times faster than the triplet-based assembly, with identical results up to floating-point summation order — and a :cpp:class:`ipc::MeshFEMHessianAssembler` held across :cpp:func:`ipc::Potential::assemble_hessian` calls additionally reuses the sparsity pattern between assemblies. It requires ``IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT=RowMajor`` (the default; the option is automatically disabled otherwise).
 
+.. note::
+    ``geogram`` makes the distance-type classifiers in ``ipc/distance/distance_type_exact.hpp`` (used by the ESP potential) exact. With ``IPC_TOOLKIT_WITH_GEOGRAM`` set to ``OFF``, the same predicates are evaluated in floating point, still consistently across the classifiers, and a warning is logged the first time they are used.
+
 .. warning::
     ``filib`` is licensed under `LGPL-2.1 <https://github.com/zfergus/filib/blob/main/LICENSE>`_ and as such it is required to be dynamically linked. Doing so automatically is a challenge, so by default we use static linkage. Enabling dynamic linkage requires copying the ``.so``/``.dylib``/``.dll`` file to the binary directory or system path. To enable this, set the CMake option ``FILIB_BUILD_SHARED_LIBS`` to ``ON`` and add this CMake code to copy the shared library object to the binary directory:
 

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

@@ -0,0 +1,5 @@
+Inverse-Power Barrier
+=====================
+
+.. doxygenclass:: ipc::InversePowerBarrier
+    :allow-dot-graphs:

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

@@ -0,0 +1,5 @@
+Near-Far Barrier
+================
+
+.. doxygenclass:: ipc::NearFarBarrier
+    :allow-dot-graphs:

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

@@ -11,6 +11,8 @@ Barrier
    ClampedLogSqBarrier.rst
    CubicBarrier.rst
    TwoStageBarrier.rst
+   InversePowerBarrier.rst
+   NearFarBarrier.rst
 
 Barrier functions and functionals.
 
@@ -34,6 +36,10 @@ Types
       - Normalized barrier function from [Li et al.].
     * - :cpp:class:`ipc::TwoStageBarrier`
       - Two-stage barrier function from [Chen et al. 2025].
+    * - :cpp:class:`ipc::InversePowerBarrier`
+      - Inverse-power barrier :math:`h(d,\hat{d})/d^p` with a C² cubic B-spline window.
+    * - :cpp:class:`ipc::NearFarBarrier`
+      - Splits a base barrier into near and far parts (used by the ESP potential).
 
 
 Functions

docs/source/tutorials/gcp.rst

@@ -297,17 +297,17 @@ GCP is implemented as separate collision and potential classes. A basic example
             double beta_n = 0.0;    // exterior direction constraint offset
             int r = 2;              // barrier exponent (dimension - 1)
 
-            ipc::SmoothContactParameters params(dhat, alpha_t, beta_t, alpha_n, beta_n, r);
+            ipc::GCPParameters params(dhat, alpha_t, beta_t, alpha_n, beta_n, r);
 
             // Build collision set
             bool use_adaptive_dhat = true;
-            ipc::SmoothCollisions collisions;
+            ipc::GCPCollisions collisions;
             if (use_adaptive_dhat)
                 collisions.compute_adaptive_dhat(collision_mesh, vertices, params);
             collisions.build(collision_mesh, vertices, params, use_adaptive_dhat);
 
             // Compute potential
-            ipc::SmoothContactPotential barrier_potential(params);
+            ipc::GCPPotential barrier_potential(params);
             double b = barrier_potential(collisions, collision_mesh, vertices);
 
             // Compute gradient
@@ -331,17 +331,17 @@ GCP is implemented as separate collision and potential classes. A basic example
             beta_n = 0.0      # exterior direction constraint offset
             r = 2             # barrier exponent (dimension - 1)
 
-            params = ipctk.SmoothContactParameters(dhat, alpha_t, beta_t, alpha_n, beta_n, r)
+            params = ipctk.GCPParameters(dhat, alpha_t, beta_t, alpha_n, beta_n, r)
 
             # Build collision set
             use_adaptive_dhat = True
-            collisions = ipctk.SmoothCollisions()
+            collisions = ipctk.GCPCollisions()
             if use_adaptive_dhat:
                 collisions.compute_adaptive_dhat(collision_mesh, vertices, params)
             collisions.build(collision_mesh, vertices, params, use_adaptive_dhat)
 
             # Compute potential
-            barrier_potential = ipctk.SmoothContactPotential(params)
+            barrier_potential = ipctk.GCPPotential(params)
             b = barrier_potential(collisions, collision_mesh, vertices)
 
             # Compute gradient
@@ -351,15 +351,15 @@ GCP is implemented as separate collision and potential classes. A basic example
             hess = barrier_potential.hessian(collisions, collision_mesh, vertices)
 
 .. important::
-    If ``use_adaptive_dhat`` is true, make sure to call ``SmoothCollisions::compute_adaptive_dhat()`` **before** ``SmoothCollisions::build()``. Adaptive :math:`\hat{d}` computes per-element barrier extents based on the rest configuration to guarantee zero potential (and zero forces) in the undeformed state.
+    If ``use_adaptive_dhat`` is true, make sure to call ``GCPCollisions::compute_adaptive_dhat()`` **before** ``GCPCollisions::build()``. Adaptive :math:`\hat{d}` computes per-element barrier extents based on the rest configuration to guarantee zero potential (and zero forces) in the undeformed state.
 
 .. note::
-    Unlike ``NormalCollisions`` in IPC, ``SmoothCollisions`` must be rebuilt whenever vertex positions change, because the interaction set depends on the current geometry (normals, tangents, and distances).
+    Unlike ``NormalCollisions`` in IPC, ``GCPCollisions`` must be rebuilt whenever vertex positions change, because the interaction set depends on the current geometry (normals, tangents, and distances).
 
 Parameter Choices
 -----------------
 
-The ``SmoothContactParameters`` structure contains the following parameters:
+The ``GCPParameters`` structure contains the following parameters:
 
 .. list-table::
    :header-rows: 1
@@ -392,7 +392,7 @@ As :math:`\alpha + \beta` decreases, the support of the Heaviside function shrin
 
 Additional internal parameters that may affect behavior:
 
-- **Adaptive dhat ratio** (default ``0.5``): Controls the ratio :math:`\epsilon(x) / d_c(x, f_0)` in the adaptive barrier localization. Set it via ``SmoothContactParameters::set_adaptive_dhat_ratio()`` in C++ or the ``SmoothContactParameters.adaptive_dhat_ratio`` property in Python.
+- **Adaptive dhat ratio** (default ``0.5``): Controls the ratio :math:`\epsilon(x) / d_c(x, f_0)` in the adaptive barrier localization. Set it via ``GCPParameters::set_adaptive_dhat_ratio()`` in C++ or the ``GCPParameters.adaptive_dhat_ratio`` property in Python.
 - **Element measure** :math:`L`: For vertices, this is set to the average edge length around the vertex; for edges, to the edge length; for faces, :math:`L` is not needed. This determines the strength of the potential for low-dimensional contact (edge–edge, edge–vertex, vertex–vertex).
 
 Friction

python/src/bindings.cpp

@@ -98,6 +98,7 @@ PYBIND11_MODULE(ipctk, m)
     define_smooth_mu(m);
 
     define_smooth_potential(m);
+    define_esp_potential(m);
 
     // geometry
     define_angle(m);

python/src/candidates/candidates.cpp

@@ -12,7 +12,7 @@ void define_candidates(py::module_& m)
             "build",
             py::overload_cast<
                 const CollisionMesh&, Eigen::ConstRef<Eigen::MatrixXd>,
-                const double, BroadPhase*>(&Candidates::build),
+                const double, BroadPhase*, const bool>(&Candidates::build),
             R"ipc_Qu8mg5v7(
             Initialize the set of discrete collision detection candidates.
 
@@ -21,9 +21,10 @@ void define_candidates(py::module_& m)
                 vertices: Surface vertex positions (rowwise).
                 inflation_radius: Amount to inflate the bounding boxes.
                 broad_phase: Broad phase to use.
+                all_types: Also detect the candidate types only ESP needs.
             )ipc_Qu8mg5v7",
             "mesh"_a, "vertices"_a, "inflation_radius"_a = 0,
-            "broad_phase"_a = nullptr)
+            "broad_phase"_a = nullptr, "all_types"_a = false)
         .def(
             "build",
             py::overload_cast<