Skip to content

Commit db256cc

Browse files
committed
Merge branch 'main' into ogc
2 parents 17084ce + 112712c commit db256cc

39 files changed

+500
-391
lines changed

docs/source/index.rst

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121
tutorials/getting_started.rst
2222
tutorials/advanced_friction.rst
2323
tutorials/convergent.rst
24+
tutorials/gcp.rst
2425
tutorials/ogc.rst
2526
tutorials/nonlinear_ccd.rst
2627
tutorials/adhesion.rst

docs/source/references.bib

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -118,3 +118,13 @@ @article{Chen2025Offset
118118
number = 4,
119119
note = {\url{https://ankachan.github.io/Projects/OGC/index.html}}
120120
}
121+
@article{Huang2025GCP,
122+
title = {Geometric Contact Potential},
123+
author = {Huang, Zizhou and Paik, Maxwell and Ferguson, Zachary and Panozzo, Daniele and Zorin, Denis},
124+
year = 2025,
125+
month = aug,
126+
journal = {ACM Transactions on Graphics},
127+
volume = 44,
128+
number = 4,
129+
note = {\url{https://doi.org/10.1145/3731142}}
130+
}

docs/source/tutorials/gcp.rst

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
.. _geometric-contact-potential-formulation:
2+
3+
Geometric Contact Potential
4+
======================
5+
6+
In addition to the original implementation of :cite:t:`Li2020IPC`, we also implement the Geometric Contact Potential (GCP) from :cite:t:`Huang2025GCP`.
7+
8+
GCP is implemented as separate collision and potential classes. A basic example of computing the GCP potential is as follows.
9+
10+
.. code-block:: c++
11+
12+
ipc::SmoothCollisions collisions;
13+
if (use_adaptive_dhat)
14+
collisions.compute_adaptive_dhat(collision_mesh, vertices, params);
15+
collisions.build(collision_mesh, vertices, params, use_adaptive_dhat);
16+
17+
ipc::SmoothContactPotential barrier_potential;
18+
double b = barrier_potential(collisions, mesh, vertices);
19+
20+
.. important::
21+
If `use_adaptive_dhat` is true, make sure to call `SmoothCollisions::compute_adaptive_dhat()` before `SmoothCollisions::build()`.
22+
23+
Technical Details
24+
-----------------
25+
26+
Similar to IPC, GCP potential has a distance barrier that converges to infinity as the distance converges to zero.
27+
28+
.. math::
29+
B(d; \hat{d}, p) := \frac{3}{2} B^3(\frac{2d}{\hat{d}}) d^{-p}
30+
31+
where :math:`\hat{d}` is the support size of the barrier function, :math:`p=d-1`, where :math:`d` is the dimension of the physics scene (2 or 3).
32+
33+
The main motivation of developing GCP is to eliminate the spurious collision forces in IPC when :math:`\hat{d}` is relatively large. The two important components to achieve this are *local minimum constraint* and *exterior direction constraint*. Please check out :cite:t:`Huang2025GCP` for more details, here we only give a high level idea.
34+
35+
The local minimum constraint is by noticing that the barrier potential is only necessary when the two primitives in a collision pair are the "distance local minima" on each surface they live. To achieve this, we consider the gradient of distance function along the tangent directions on the surface, and smoothly clamp the barrier to zero if the gradient of distance is positive along the tangent directions, since it suggests that there already exists a closer point on the surface that can trigger the barrier.
36+
37+
The exterior direction constraint is by noticing that, for a thin volumetric shell, the barrier potential is not needed for points on the opposite sides of the shell. To achieve this, we smoothly clamp the barrier based on the normal directions of two primitives in a collision pair. If the normal direction points to the opposite direction of the distance direction, we can skip this collision pair.
38+
39+
Both constraints not only reduces the number of collision pairs, but also increases the accuracy of collision handling by eliminating spurious forces that are non-physical, allowing wider range of :math:`\hat{d}` to be used in the simulation.
40+
41+
Parameter choices
42+
-----------------
43+
44+
Besides :math:`\hat{d}`, GCP also has other parameters to control the smoothness of local minimum and exterior direction constraints. For each constraint, there are :math:`\alpha` and :math:`\beta` to control the support region of the smooth Heaviside function (:math:`\alpha_t,\ \beta_t` for local minimum constraint, and :math:`\alpha_n,\ \beta_n` for exterior direction constraint). The smooth Heaviside function :math:`H(z;\alpha,\beta)\in C^1(\mathbb{R})` satisfies that
45+
46+
.. math::
47+
H(z;\alpha,\beta) = 0,\ \forall z < -\alpha \\
48+
H(z;\alpha,\beta) = 1,\ \forall z > \beta \\
49+
H'(z;\alpha,\beta) \geq 0,\ \forall z
50+
51+
Since :math:`H(z;\alpha,\beta)` takes the dot and cross products of unit vectors as inputs, we have :math:`-1\leq z\leq 1`. Thus, the basic requirement for parameters is:
52+
53+
.. math::
54+
0 \leq -\beta < \alpha \leq 1.
55+
56+
As :math:`\alpha + \beta` decreases, the support size of the Heaviside function reduces and the function becomes less smooth, thus, the overall potential becomes harder to optimize.
57+
58+
.. tip::
59+
For simplicity, we recommend to just use :math:`\beta_n = \beta_t = 0`. To make sure the potential landscape is smooth, the recommended parameter choices are :math:`\alpha_t\in [0.2, 0.9]`, :math:`\alpha_n = 0.1`.
60+
61+
Friction
62+
--------
63+
64+
We implement the original friction formulation of :cite:t:`Li2020IPC` following the same style, details are covered in the paper :cite:t:`Huang2025GCP`.

src/ipc/collisions/tangential/tangential_collisions.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -116,7 +116,7 @@ void TangentialCollisions::build(
116116

117117
void TangentialCollisions::build(
118118
const CollisionMesh& mesh,
119-
const Eigen::MatrixXd& vertices,
119+
Eigen::ConstRef<Eigen::MatrixXd> vertices,
120120
const SmoothCollisions& collisions,
121121
const SmoothContactParameters& params,
122122
const double normal_stiffness,

src/ipc/collisions/tangential/tangential_collisions.hpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -97,7 +97,7 @@ class TangentialCollisions {
9797
/// @param blend_mu Function to blend vertex-based coefficients of friction. Defaults to average.
9898
void build(
9999
const CollisionMesh& mesh,
100-
const Eigen::MatrixXd& vertices,
100+
Eigen::ConstRef<Eigen::MatrixXd> vertices,
101101
const SmoothCollisions& collisions,
102102
const SmoothContactParameters& params,
103103
const double normal_stiffness,

src/ipc/math/math.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -49,7 +49,7 @@ double opposite_direction_penalty(
4949
return Math<double>::smooth_heaviside(d.dot(t) / t.norm(), alpha, beta);
5050
}
5151

52-
GradType<6> opposite_direction_penalty_grad(
52+
GradientType<6> opposite_direction_penalty_grad(
5353
Eigen::ConstRef<Eigen::Vector3d> t,
5454
Eigen::ConstRef<Eigen::Vector3d> d,
5555
const double alpha,
@@ -109,7 +109,7 @@ double negative_orientation_penalty(
109109
return opposite_direction_penalty(n, d, alpha, beta);
110110
}
111111

112-
GradType<9> negative_orientation_penalty_grad(
112+
GradientType<9> negative_orientation_penalty_grad(
113113
Eigen::ConstRef<Eigen::Vector3d> t1,
114114
Eigen::ConstRef<Eigen::Vector3d> t2,
115115
Eigen::ConstRef<Eigen::Vector3d> d,

src/ipc/math/math.hpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -119,7 +119,7 @@ double opposite_direction_penalty(
119119
const double alpha,
120120
const double beta);
121121

122-
GradType<6> opposite_direction_penalty_grad(
122+
GradientType<6> opposite_direction_penalty_grad(
123123
Eigen::ConstRef<Eigen::Vector3d> t,
124124
Eigen::ConstRef<Eigen::Vector3d> d,
125125
const double alpha,
@@ -139,7 +139,7 @@ double negative_orientation_penalty(
139139
const double alpha,
140140
const double beta);
141141

142-
GradType<9> negative_orientation_penalty_grad(
142+
GradientType<9> negative_orientation_penalty_grad(
143143
Eigen::ConstRef<Eigen::Vector3d> t1,
144144
Eigen::ConstRef<Eigen::Vector3d> t2,
145145
Eigen::ConstRef<Eigen::Vector3d> d,

src/ipc/potentials/potential.hpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ template <class TCollisions> class Potential {
1212
using TCollision = typename TCollisions::value_type;
1313
/// @brief Maximum degrees of freedom per collision
1414
static constexpr int STENCIL_NDOF = 3 * TCollision::STENCIL_SIZE;
15-
using VectorMaxNd = Vector<double, Eigen::Dynamic, STENCIL_NDOF>;
15+
using VectorMaxNd = VectorMax<double, STENCIL_NDOF>;
1616
using MatrixMaxNd = MatrixMax<double, STENCIL_NDOF, STENCIL_NDOF>;
1717

1818
public:

src/ipc/potentials/tangential_potential.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -683,7 +683,7 @@ TangentialPotential::VectorMaxNd TangentialPotential::smooth_contact_force(
683683
assert(rest_positions.size() == velocities.size());
684684

685685
// const VectorMaxNd x = dof(rest_positions, edges, faces);
686-
// const Vector<double, -1, STENCIL_NDOF> u =
686+
// const VectorMax<double, STENCIL_NDOF> u =
687687
// dof(lagged_displacements, edges, faces);
688688
// const VectorMaxNd v = dof(velocities, edges, faces);
689689
const VectorMaxNd lagged_positions =

0 commit comments

Comments
 (0)