Skip to content

Commit 86af88a

Browse files
zferguszfergus
committed
Separate convergent formulation tutorial
1 parent 3db1d97 commit 86af88a

File tree

3 files changed

+112
-112
lines changed

3 files changed

+112
-112
lines changed

‎docs/source/index.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@
1111
:hidden:
1212

1313
tutorial/getting_started.rst
14+
tutorial/convergent.rst
1415
tutorial/nonlinear_ccd.rst
1516
tutorial/simulation.rst
1617
tutorial/misc.rst

‎docs/source/tutorial/convergent.rst‎

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
.. _convergent-contact-formulation:
2+
3+
Convergent Formulation
4+
======================
5+
6+
In addition to the original implementation of :cite:t:`Li2020IPC`, we also implement the convergent formulation of :cite:t:`Li2023Convergent`.
7+
8+
To enable the convergent formulation, we need to set ``use_convergent_formulation`` before building ``Constraints``:
9+
10+
.. md-tab-set::
11+
12+
.. md-tab-item:: C++
13+
14+
.. code-block:: c++
15+
16+
collision_constraints.set_use_convergent_formulation(true);
17+
collision_constraints.build(collision_mesh, vertices, dhat);
18+
19+
.. md-tab-item:: Python
20+
21+
.. code-block:: python
22+
23+
collision_constraints.use_convergent_formulation = True
24+
collision_constraints.build(collision_mesh, vertices, dhat)
25+
26+
.. important::
27+
The variable ``use_convergent_formulation`` should be set before calling ``CollisionConstraints::build`` for it to take effect. By default, it is ``false``.
28+
29+
Technical Details
30+
-----------------
31+
32+
*We briefly summarize the convergent formulation here for convenience.*
33+
34+
In order to derive a convergent formulation, we first define a continuous form of our barrier potential :math:`P`. For a surface :math:`\mathcal{S}` embedded in 3D space, we parameterize the surfaces by common (possibly discontinuous) coordinates :math:`u \in \tilde{M} \subset \mathbb{R}^2`, so that :math:`\mathcal{S}(u)` traverses the material points across all surfaces contiguously. The total contact potential is then
35+
36+
.. math::
37+
P(\mathcal{S})=\frac{1}{2} \int_{u \in \tilde{M}} \max _{v \in \tilde{M} \setminus{ }_r u} b(d(\mathcal{S}(u), \mathcal{S}(v)), \hat{d})~\mathrm{d} u,
38+
39+
where we define the operator :math:`\setminus_r: \mathcal{P}(\mathbb{R}^2) \times \mathbb{R} \times \mathbb{R}^2 \mapsto \mathcal{P}(\mathbb{R}^2)` to be
40+
41+
.. math::
42+
\tilde{M} \setminus_r u:=\left\{v \in \tilde{M} \mid\|u-v\|_2>r\right\}
43+
44+
with :math:`r \rightarrow 0`.
45+
46+
We then define our surface discretization with a triangulated boundary mesh geometry. As in the smooth case, we can parameterize the domain across all triangles with :math:`u \in \tilde{M}` so that :math:`p(u): \tilde{M} \mapsto \mathbb{R}^3` traverses all material points, across all triangles :math:`t ∈ T` in the triangle mesh contiguously. The corresponding surface contact potential is then
47+
48+
.. math::
49+
\frac{1}{2} \int_{u \in \tilde{M}} \max_{t \in T \backslash p(u)} b(d(p(u), t), \hat{d})~\mathrm{d} u,
50+
51+
where :math:`T \setminus p` is the set of boundary faces that do not contain the point :math:`p`.
52+
53+
Applying mesh vertices as nodes (and quadrature points), we numerically integrate the surface contact potential. For each nodal position :math:`x \in V` we then have a corresponding material space coordinate :math:`\bar{x} \in \bar{V}`. Piecewise linear integration of the surface barrier is then
54+
55+
.. math::
56+
\frac{1}{2} \sum_{\bar{x} \in \bar{V}} w_{\bar{x}} \max _{t \in T \backslash x(\bar{x})} b(d(x(\bar{x}), t), \hat{d}),
57+
58+
where :math:`w_{\bar{x}}` are the quadrature weights, each given by one-third of the sum of the areas (in material space) of the boundary triangles incident to :math:`\bar{x}`.
59+
60+
We next need to smoothly approximate the max operator in the contact potentials. However, common approaches such as an :math:`L^p`-norm or LogSumExp would decrease sparsity in subsequent numerical solves by increasing the stencil size per contact evaluation. We instead leverage the locality of our barrier function to approximate the max operator by removing duplicate distance pairs. Our resulting approximators for a triangulated surface is
61+
62+
.. math::
63+
\begin{aligned}
64+
\Psi_s(x) & =\sum_{t \in T \backslash x} b(d(x, t), \hat{d})-\sum_{e \in E_{\text {int }} \backslash x} b(d(x, e), \hat{d})+\sum_{x_2 \in V_{i n t} \backslash x} b\left(d\left(x, x_2\right), \hat{d}\right) \\
65+
& \approx \max _{t \in T \backslash x} b(d(x, t), \hat{d}),
66+
\end{aligned}
67+
68+
where :math:`V_{\text{int}} \subseteq V` is the subset of internal surface nodes and :math:`E_{\text{int}} \subseteq E` is the subset of internal surface edges (i.e., edges incident to two triangles). For locally convex regions this estimator is tight while remaining smooth. In turn, for nonconvex regions, it improves over direct summation.
69+
70+
The corresponding discrete barrier potential is then simply
71+
72+
.. math::
73+
P_s(V)= \frac{1}{2} \sum_{x \in V} w_x \Psi_s(x),
74+
75+
where we simplify with :math:`w_x = w_{\bar{x}}` defined appropriately, per domain, as covered above.
76+
77+
Please, see the `paper <https://arxiv.org/abs/2307.15908>`_ for more details (including the formulation for 2D curves and edge-edge contact) and evaluation.
78+
79+
The key difference between the original and the convergent formulations is that we (1) include area weights in the barrier potential and (2) include additional (negative) terms to cancel out the duplicate distance pairs. While this requires minor algorithmic changes, it produces considerably better results.
80+
81+
Physical Barrier
82+
----------------
83+
84+
We want our barrier potential to have the same units as our elastic potential (e.g., :math:`\text{J}`). Together with the area weighting (discussed above), this means the barrier should have units of pressure times distance (e.g., :math:`\text{Pa} \cdot \text{m}`). That is,
85+
86+
.. math::
87+
\text{Pa} \cdot \text{m} \cdot \text{m}^2 = \frac{\text{N}}{\text{m}^2} \cdot \text{m} \cdot \text{m}^2 = \text{N} \cdot \text{m} = \text{J}.
88+
89+
To achieve this, (when using the convergent formulation) we modify the barrier function to have units of distance:
90+
91+
.. math::
92+
b(d, \hat{d})=\left\{\begin{array}{lr}
93+
-\hat{d}\left(\frac{d}{\hat{d}}-1\right)^2 \ln \left(\frac{d}{\hat{d}}\right), & 0<d<\hat{d} \\
94+
0 & d \geq \hat{d}
95+
\end{array}\right.
96+
97+
.. note::
98+
This is equivalent to the original barrier function of :cite:p:`Li2020IPC` times :math:`1/\hat{d}^3` when using squared distances. Therefore, to simplify the implementation we only implement the original barrier function and multiply all constraints by :math:`1/\hat{d}^3`.
99+
100+
The barrier stiffness (:math:`\kappa`) then has units of pressure (e.g., :math:`\text{Pa}`), the same as Young's modulus (:math:`E`) in elasticity.
101+
This implies we can get good solver convergence even when using a fixed :math:`\kappa` by setting it relative to the material's Young's modulus (:math:`\kappa = 0.1 E` works well in many examples).
102+
The intention is to treat the barrier as a thin elastic region around the mesh, and having consistent units makes it easier to pick the stiffness for this "material".
103+
104+
.. _convergent-friction-formulation:
105+
106+
Friction
107+
--------
108+
109+
Just as with the :ref:`collision constraints <convergent-contact-formulation>`, we implement both the original friction formulation of :cite:t:`Li2020IPC` and the convergent formulation of :cite:t:`Li2023Convergent`.
110+
111+
The choice of formulation is dependent on how the fixed set of ``contact_constraints`` given to ``FrictionConstraints::build`` was built. If the ``contact_constraints`` were built using the convergent formulation, then the friction constraints will also use the convergent formulation. Otherwise, the original formulation will be used.

‎docs/source/tutorial/getting_started.rst‎

Lines changed: 0 additions & 112 deletions
Original file line numberDiff line numberDiff line change
@@ -327,109 +327,6 @@ To compute the adaptive barrier stiffness, we can use two functions: ``initial_b
327327
328328
# (next iteration)
329329
330-
.. _convergent-contact-formulation:
331-
332-
Convergent Formulation
333-
^^^^^^^^^^^^^^^^^^^^^^
334-
335-
In addition to the original implementation of :cite:t:`Li2020IPC`, we also implement the convergent formulation of :cite:t:`Li2023Convergent`.
336-
337-
To enable the convergent formulation, we need to set ``use_convergent_formulation``:
338-
339-
.. md-tab-set::
340-
341-
.. md-tab-item:: C++
342-
343-
.. code-block:: c++
344-
345-
collision_constraints.set_use_convergent_formulation(true);
346-
collision_constraints.build(collision_mesh, vertices, dhat);
347-
348-
.. md-tab-item:: Python
349-
350-
.. code-block:: python
351-
352-
collision_constraints.use_convergent_formulation = True
353-
collision_constraints.build(collision_mesh, vertices, dhat)
354-
355-
.. important::
356-
The variable ``use_convergent_formulation`` should be set before calling ``CollisionConstraints::build`` for it to take effect. By default, it is ``false``.
357-
358-
Technical Details
359-
'''''''''''''''''
360-
361-
*We briefly summarize the convergent formulation here for convenience.*
362-
363-
In order to derive a convergent formulation, we first define a continuous form of our barrier potential :math:`P`. For a surface :math:`\mathcal{S}` embedded in 3D space, we parameterize the surfaces by common (possibly discontinuous) coordinates :math:`u \in \tilde{M} \subset \mathbb{R}^2`, so that :math:`\mathcal{S}(u)` traverses the material points across all surfaces contiguously. The total contact potential is then
364-
365-
.. math::
366-
P(\mathcal{S})=\frac{1}{2} \int_{u \in \tilde{M}} \max _{v \in \tilde{M} \setminus{ }_r u} b(d(\mathcal{S}(u), \mathcal{S}(v)), \hat{d})~\mathrm{d} u,
367-
368-
where we define the operator :math:`\setminus_r: \mathcal{P}(\mathbb{R}^2) \times \mathbb{R} \times \mathbb{R}^2 \mapsto \mathcal{P}(\mathbb{R}^2)` to be
369-
370-
.. math::
371-
\tilde{M} \setminus_r u:=\left\{v \in \tilde{M} \mid\|u-v\|_2>r\right\}
372-
373-
with :math:`r \rightarrow 0`.
374-
375-
We then define our surface discretization with a triangulated boundary mesh geometry. As in the smooth case, we can parameterize the domain across all triangles with :math:`u \in \tilde{M}` so that :math:`p(u): \tilde{M} \mapsto \mathbb{R}^3` traverses all material points, across all triangles :math:`t ∈ T` in the triangle mesh contiguously. The corresponding surface contact potential is then
376-
377-
.. math::
378-
\frac{1}{2} \int_{u \in \tilde{M}} \max_{t \in T \backslash p(u)} b(d(p(u), t), \hat{d})~\mathrm{d} u,
379-
380-
where :math:`T \setminus p` is the set of boundary faces that do not contain the point :math:`p`.
381-
382-
Applying mesh vertices as nodes (and quadrature points), we numerically integrate the surface contact potential. For each nodal position :math:`x \in V` we then have a corresponding material space coordinate :math:`\bar{x} \in \bar{V}`. Piecewise linear integration of the surface barrier is then
383-
384-
.. math::
385-
\frac{1}{2} \sum_{\bar{x} \in \bar{V}} w_{\bar{x}} \max _{t \in T \backslash x(\bar{x})} b(d(x(\bar{x}), t), \hat{d}),
386-
387-
where :math:`w_{\bar{x}}` are the quadrature weights, each given by one-third of the sum of the areas (in material space) of the boundary triangles incident to :math:`\bar{x}`.
388-
389-
We next need to smoothly approximate the max operator in the contact potentials. However, common approaches such as an :math:`L^p`-norm or LogSumExp would decrease sparsity in subsequent numerical solves by increasing the stencil size per contact evaluation. We instead leverage the locality of our barrier function to approximate the max operator by removing duplicate distance pairs. Our resulting approximators for a triangulated surface is
390-
391-
.. math::
392-
\begin{aligned}
393-
\Psi_s(x) & =\sum_{t \in T \backslash x} b(d(x, t), \hat{d})-\sum_{e \in E_{\text {int }} \backslash x} b(d(x, e), \hat{d})+\sum_{x_2 \in V_{i n t} \backslash x} b\left(d\left(x, x_2\right), \hat{d}\right) \\
394-
& \approx \max _{t \in T \backslash x} b(d(x, t), \hat{d}),
395-
\end{aligned}
396-
397-
where :math:`V_{\text{int}} \subseteq V` is the subset of internal surface nodes and :math:`E_{\text{int}} \subseteq E` is the subset of internal surface edges (i.e., edges incident to two triangles). For locally convex regions this estimator is tight while remaining smooth. In turn, for nonconvex regions, it improves over direct summation.
398-
399-
The corresponding discrete barrier potential is then simply
400-
401-
.. math::
402-
P_s(V)= \frac{1}{2} \sum_{x \in V} w_x \Psi_s(x),
403-
404-
where we simplify with :math:`w_x = w_{\bar{x}}` defined appropriately, per domain, as covered above.
405-
406-
Please, see the `paper <https://arxiv.org/abs/2307.15908>`_ for more details (including the formulation for 2D curves and edge-edge contact) and evaluation.
407-
408-
The key difference between the original and the convergent formulations is that we (1) include area weights in the barrier potential and (2) include additional (negative) terms to cancel out the duplicate distance pairs. While this requires minor algorithmic changes, it produces considerably better results.
409-
410-
Physical Barrier
411-
''''''''''''''''
412-
413-
We want our barrier potential to have the same units as our elastic potential (e.g., :math:`\text{J}`). Together with the area weighting (discussed above), this means the barrier should have units of pressure times distance (e.g., :math:`\text{Pa} \cdot \text{m}`). That is,
414-
415-
.. math::
416-
\text{Pa} \cdot \text{m} \cdot \text{m}^2 = \frac{\text{N}}{\text{m}^2} \cdot \text{m} \cdot \text{m}^2 = \text{N} \cdot \text{m} = \text{J}.
417-
418-
To achieve this, (when using the convergent formulation) we modify the barrier function to have units of distance:
419-
420-
.. math::
421-
b(d, \hat{d})=\left\{\begin{array}{lr}
422-
-\hat{d}\left(\frac{d}{\hat{d}}-1\right)^2 \ln \left(\frac{d}{\hat{d}}\right), & 0<d<\hat{d} \\
423-
0 & d \geq \hat{d}
424-
\end{array}\right.
425-
426-
.. note::
427-
This is equivalent to the original barrier function of :cite:p:`Li2020IPC` times :math:`1/\hat{d}^3` when using squared distances. Therefore, to simplify the implementation we only implement the original barrier function and multiply all constraints by :math:`1/\hat{d}^3`.
428-
429-
The barrier stiffness (:math:`\kappa`) then has units of pressure (e.g., :math:`\text{Pa}`), the same as Young's modulus (:math:`E`) in elasticity.
430-
This implies we can get good solver convergence even when using a fixed :math:`\kappa` by setting it relative to the material's Young's modulus (:math:`\kappa = 0.1 E` works well in many examples).
431-
The intention is to treat the barrier as a thin elastic region around the mesh, and having consistent units makes it easier to pick the stiffness for this "material".
432-
433330
.. _modeling-thickness:
434331

435332
Modeling Thickness
@@ -562,15 +459,6 @@ We can also compute the first and second derivatives of the friction dissipative
562459
friction_potential_hess = friction_constraints.compute_potential_hessian(
563460
collision_mesh, velocity, epsv)
564461
565-
.. _convergent-friction-formulation:
566-
567-
Convergent Formulation
568-
^^^^^^^^^^^^^^^^^^^^^^
569-
570-
Just as with the :ref:`collision constraints <convergent-contact-formulation>`, we implement both the original friction formulation of :cite:t:`Li2020IPC` and the convergent formulation of :cite:t:`Li2023Convergent`.
571-
572-
The choice of formulation is dependent on how the fixed set of ``contact_constraints`` given to ``FrictionConstraints::build`` was built. If the ``contact_constraints`` were built using the convergent formulation, then the friction constraints will also use the convergent formulation. Otherwise, the original formulation will be used.
573-
574462
Continuous Collision Detection
575463
------------------------------
576464

0 commit comments

Comments
 (0)