弹性体（布料）仿真方法总结

本文禁止转载

B站：Heskey0

Dynamic Deformables: Implementation and Production Practicalities

Dynamic Deformables

Chapter 1. Introduction

Stable Neo-Hookean Flesh Simulation, Breannan Smith, Fernando de Goes, Theodore Kim (Technical Paper presented at SIGGRAPH 2018)
Better Collisions and Faster Cloth for Pixar’s Coco, David Eberle (Talk presented at SIGGRAPH 2018)
Clean Cloth Inputs: Removing Character Self-Intersections with Volume Simulation, Audrey Wong, David Eberle, Theodore Kim (Talk presented at SIGGRAPH 2018)
Robust Skin Simulation in Incredibles 2, Ryan Kautzman, Gordon Cameron, Theodore Kim (Talk presented at SIGGRAPH 2018)
Analytic Eigensystems for Isotropic Distortion Energies, Breannan Smith, Fernando de Goes, Theodore Kim (Technical Paper presented at SIGGRAPH 2019)
Anisotropic Elasticity for Inversion-Safety and Element Rehabilitation, Theodore Kim, Fernando de Goes, Hayley Iben (Technical Paper presented at SIGGRAPH 2019)

Chapter 2. Deformation Fundamentals

\(\Psi\): energy functions (scoring function)
\(tr(A)=a_{11}+a_{22}+...+a_{nn}\): trace
\(\phi(x)=Fx+t\): deformation map (affine function)
\(F=D_SD_m^{-1}\): deformation gradient
- polar decomposition: \(F=RS\).
\(C=F^TF\): The right Cauchy-Green tensor
\(E=\frac12(F^TF-I)\): Green's strain

Chapter 3. Tensors

double-contraction: (2nd-order + 2nd-order)
- \[A:B= \pmatrix{a_0&a_2\\a_1&a_3} \pmatrix{b_0&b_2\\b_1&b_3} =a_0b_0+a_1b_1+a_2b_2+a_3b_3 \]
double-contraction: (3rd-order + 2nd-order)
- \[A:B= \pmatrix{ \pmatrix{a_0&a_2\\a_1&a_3}\\ \pmatrix{a_4&a_6\\a_5&a_7}\\ \pmatrix{a_8&a_{10}\\a_9&a_{11}} } \pmatrix{b_0&b_2\\b_1&b_3} = \pmatrix{ a_0b_0+a_1b_1+a_2b_2+a_3b_3\\ a_4b_0+a_5b_1+a_6b_2+a_7b_3\\ a_8b_0+a_9b_1+a_{10}b_2+a_{11}b_3 } \]
- \((3\times1\times2\times2):(2\times2)=(3\times1)\).
- \(\part F/\part x\) is a 3rd-order tensor, not a matrix
flattening/vectorization operator \(vec(\cdot)\): convert any matrix to a vector, and any higher-order tensor to a matrix.

\[vec(A)=vec(\pmatrix{a_0&a_2\\a_1&a_3})=\pmatrix{a_0\\a_1\\a_2\\a_3} \]

for 3rd-order tensors:

\[vec(A)=vec(\pmatrix{ \pmatrix{A}\\ \pmatrix{B}\\ \pmatrix{C} })=\pmatrix{vec(A)&vec(B)&vec(C)} \]

for 4th-order tensors:

\[vec(A)=vec(\pmatrix{ \pmatrix{A} \pmatrix{B}\\ \pmatrix{C} \pmatrix{D} })=\pmatrix{vec(A)&vec(B)&vec(C)&vec(D)} \]

note that:

\[(vec(A))^Tvec(B)=A:B \]

Chapter 4. Force Gradients

\[\frac{\part \Psi}{\part x}=vec(\frac{\part F}{\part x})^Tvec(\frac{\part \Psi}{\part F})\\ \frac{\part \Psi}{\part x}=vec(\frac{\part E}{\part x})^Tvec(\frac{\part \Psi}{\part E})\\ \frac{\part^2 \Psi}{\part x^2}=vec(\frac{\part F}{\part x})^Tvec(\frac{\part^2 \Psi}{\part F^2})vec(\frac{\part F}{\part x}) \]

first Piola-Kirchhoff stress tensor:(PK1)

\[P=\frac{\part \Psi}{\part F} \]
second Piola-Kirchhoff stress tensor(PK2):

\[\frac{\part\Psi}{\part E} \]
energy Hessian:

\[H=vec(\frac{\part^2 \Psi}{\part F^2}) \]

4.1 Energy

Dirichlet energy:

\[\Psi=||F||^2_F \]

St. Venant Kirchhoff(stvk): (stretching term only)

\[\Psi=||E||_F^2 \]

St. Venant Kirchhoff(stvk): (complete)

\[\Psi=\mu||E||_F^2+\frac\lambda2(tr(E))^2 \]

The Lame Parameters: Yong's Modulus: \(E\), Poisson's ratio: \(\nu\)

\[\mu=\frac{E}{2(1+\nu)}\\ \lambda=\frac{E\nu}{(1+\nu)(1-2\nu)} \]

As-Rigit-As-Possible(ARAP):

\[\Psi=\frac\mu2||F-R||_F^2 \]

Neo-Hookean:

\[\Psi=\frac\mu2(||F||_F^2-3)-\mu\log(J)+\frac\lambda2(\log(J))^2 \]

\(J=\det(F)\)

4.2 Gradient (PK1)

gradient of the Dirichlet energy:

\[\frac{\part\Psi}{\part F}=2F \]

gradient of the St. Venant Kirchhoff: (stretching term only)

\[\frac{\part\Psi}{\part F}=FE \]

gradient of the St. Venant Kirchhoff: (complete)

\[\frac{\part\Psi}{\part F}=\mu FE+\lambda(tr(E))F \]

gradient of the As-Rigid-As-Possible(ARAP):

\[\frac{\part\Psi}{\part F}=\mu(F-R) \]

gradient of the Neo-Hookean:

\[\frac{\part\Psi}{\part F}=\mu\left(F-\frac1J\frac{\part J}{\part F}\right)+\lambda\frac{\log J}{J}\frac{\part J}{\part F} \]

3D场景中，设：\(3\times3\) 的矩阵 \(F=[f_0|f_1|f_2]\)，则：(\(\times\): cross product)

\[\frac{\part J}{\part F}=[f_1\times f_2|f_2\times f_0|f_0\times f_1] \]

4.3 Hessian

Hessian of the Dirichlet energy:

\[vec(\frac{\part P}{\part F})=2I_{9\times9} \]

Hessian of the Neo-Hookean:

\[vec(\frac{\part P}{\part F})=\mu I_{9\times9}+\left[\frac{\mu+\lambda(1-\log J)}{J^2}\right]g_Jg_J^T+\left[\frac{\lambda\log J-\mu}{J}\right]H_J \]

in which:

\[g_J=vec(\frac{\part J}{\part F}) \\H_J=vec(\frac{\part^2 J}{\part F^2}) \]
Gradient of \(J\):

\[\frac{\part J}{\part F}=\pmatrix{ f_{11}f_{22}-f_{12}f_{21}&f_{20}f_{12}-f_{22}f_{10}&f_{10}f_{21}-f_{11}f_{20}\\ f_{21}f_{02}-f_{22}f_{01}&f_{00}f_{22}-f_{02}f_{20}&f_{20}f_{01}-f_{21}f_{00}\\ f_{01}f_{12}-f_{02}f_{11}&f_{10}f_{02}-f_{12}f_{00}&f_{00}f_{11}-f_{01}f_{10} }\\ g_J=vec(\pmatrix{f_1\times f_2|f_2\times f_0|f_0\times f_1})\\ G_J=\pmatrix{ \pmatrix{\frac{\part J}{\part f_{00}}\frac{\part J}{\part F}}&\pmatrix{\frac{\part J}{\part f_{01}}\frac{\part J}{\part F}}&\pmatrix{\frac{\part J}{\part f_{02}}\frac{\part J}{\part F}}\\ \pmatrix{\frac{\part J}{\part f_{10}}\frac{\part J}{\part F}}&\pmatrix{\frac{\part J}{\part f_{11}}\frac{\part J}{\part F}}&\pmatrix{\frac{\part J}{\part f_{12}}\frac{\part J}{\part F}}\\ \pmatrix{\frac{\part J}{\part f_{20}}\frac{\part J}{\part F}}&\pmatrix{\frac{\part J}{\part f_{21}}\frac{\part J}{\part F}}&\pmatrix{\frac{\part J}{\part f_{22}}\frac{\part J}{\part F}} }\\ g_Jg_J^T=vec(G_J) \]
Hessian of \(J\): 3D场景中，设：\(3\times3\) 的矩阵 \(F=[f_0|f_1|f_2]\)，则：

\[H_J=\pmatrix{ 0_{3\times3}&-\hat f_2&\hat f_1\\ \hat f_2&0_{3\times3}&-\hat f_0\\ -\hat f_1&\hat f_0&0_{3\times3} }\\ \hat f_i=\pmatrix{ 0&-f_i[2]&f_i[1]\\ f_i[2]&0&-f_i[0]\\ -f_i[1]&f_i[0]&0 } \]

Chapter 5. A Better Way for Isotropic Solids

5.1 Invariants

Invariants with respect to what?

Understanding these quantities intuitively is about to become important, so let’s look at two different ways to think about them.

5.1.1 Invariants as rotation removers

we looked at two different approaches to removing the rotation.

StVK: (stretch only)

\[\Psi=\frac12||F^TF-I||_F^2=\frac12||S^2-I||_F^2 \]

since \(S\) is symmetric,

\[F^TF=(RS)^TRS=S^TR^TRS=S^2 \]

得出 \(I_C=tr(F^TF)=tr(S^2)\).

ARAP:

\[\begin{aligned} \Psi&=||F-R||_F^2\\ &=||R^T(F-R)||_F^2\\ &=||S-I||_F^2 \end{aligned} \]

Using rotation invariance of \(||\cdot||_F\):

\[||F-R||_F^2=||R^T(F-R)||_F^2 \]

Then comparing the two ways:

The ultimate measure of how much stretching and squashing is really going on in \(F\) can be obtained using its SVD, which is \(F=U\Sigma V^T\). \(\Sigma\) shows exactly how much streching is going on along each of the 3D directions. We can interpret the invariants in terms of these singular values.

\[\Sigma=\pmatrix{\sigma_x&0&0\\0&\sigma_y&0\\0&0&\sigma_z} \]

由于 \(S=V^T\Sigma V\)，可以推出

StVK以及ARAP的能量分别为 \(\Psi=1/2||\Sigma^2-I||^2_F\) 和 \(\Psi=||\Sigma-I||_F^2\)，
\(I_C=tr(F^TF)=tr(S)=tr(\Sigma)\).

5.1.2 Invariants as Geometric Measurements

Looking at the \(tr(\Sigma^2)\) version of \(I_C\), this works out to:

\[I_C=\sigma_x^2+\sigma_y^2+\sigma_z^2 \]

the singular values are amount of stretching and squashing along each axis:

The matrix \(F\) describes the rotation and scaling of an infinitesimal cube of material. The singular values of \(F\), the \(σ_{\{x,y,z\}}\) labels above, describe the lengths of each side of that cube.

\(I_C\) sums up the squared lengths of each side of the cube: \(I_C = σ^2_x + σ^2_y + σ^2_z\) . I.e. it measures the overall edge lengths of the entire cube.

\(I_C,II_C,II_C^\star,III_C\) invariant can be written as (Cauchy-Green invariants):

\[I_C = σ^2_x + σ^2_y + σ^2_z\\ II_C=||F^TF||_F^2=tr(\Sigma^4)=\sigma_x^4+\sigma_y^4+\sigma_z^4\\ II_C^\star=\frac12(I_C^2-II_C)=(\sigma_x\sigma_y)^2+(\sigma_x\sigma_z)^2+(\sigma_y\sigma_z)^2\\ III_C=\det(F^TF)=(\sigma_x\sigma_y\sigma_z)^2 \]

上面几个式子全是平方，You can’t use a sum of squared values to express a sum of unsquared values，所以我们定义一组新的invariants.（section 5.3）

5.3 Building a Generic Hessian

SVD: \(F=U\Sigma V^T\).

\[\Sigma=\pmatrix{\sigma_x&0&0\\0&\sigma_y&0\\0&0&\sigma_z} \]

eighenpairs:

\[\begin{aligned} &\lambda_0=\frac2{\sigma_x+\sigma_y}&Q_0=\frac1{\sqrt{2}}U\pmatrix{0&-1&0\\1&0&0\\0&0&0}V^T\\ &\lambda_1=\frac2{\sigma_y+\sigma_z}&Q_1=\frac1{\sqrt{2}}U\pmatrix{0&0&0\\0&0&1\\0&-1&0}V^T\\ &\lambda_2=\frac2{\sigma_x+\sigma_z}&Q_2=\frac1{\sqrt{2}}U\pmatrix{0&0&1\\0&0&0\\-1&0&0}V^T \end{aligned} \]

invariants from Smith et al. (2019):

\[I_1=tr(S)\\ I_2=tr(S^2)\\ I_3=det(F)\\ g_1=vec(R)\\ g_2=vec(2F)\\ g_J=vec(\pmatrix{f_1\times f_2|f_2\times f_0|f_0\times f_1})\\ H_1=\sum^2_{i=0}\lambda_ivec(Q_i)vec(Q_i)^T\\ H_2=2I_{9\times9}\\ H_J=\pmatrix{ 0_{3\times3}&-\hat f_2&\hat f_1\\ \hat f_2&0_{3\times3}&-\hat f_0\\ -\hat f_1&\hat f_0&0_{3\times3} }\\ \]

energy Hessian:

\[vec(\frac{\part^2\Psi}{\part F^2})=\sum^3_{i=1}\left(\frac{\part^2\Psi}{\part I^2_i}g_ig_i^T+\frac{\part\Psi}{\part I_i}H_i\right)) \]

Finally, we follow the same three-step process:

Rewrite \(\Psi\) using \(I_1,I_2,I_3\).
Derive the scalar derivatives \(\frac{\part^2\Psi}{\part I_i},\frac{\part^2\Psi}{\part I_i^2}\).
compute \(vec(\frac{\part^2\Psi}{\part F^2})\).

5.3.1 Neo-Hookean

using \(I_3\) instead of \(J\):

\[vec(\frac{\part^2\Psi}{\part F^2})=\mu I_{9\times9}+\frac{\lambda(1-\log I_3)+\mu}{I_3^2}g_3g_3^T+\frac{\lambda\log I_3}{I_3}H_3 \]

5.3.2 ARAP

\[vec(\frac{\part^2\Psi}{\part F^2})=2I_{9\times9}-2H_1 \]

Matlab Code:

function [H] = ARAP_Hessian (F)
[U Sigma V] = svd_rv(F);

% get the twist modes
T0 = [0 -1 0; 1 0 0; 0 0 0];
T0 = (1 / sqrt(2)) * U * T0 * V’;

T1 = [0 0 0; 0 0 1; 0 -1 0];
T1 = (1 / sqrt(2)) * U * T1 * V’;

T2 = [0 0 1; 0 0 0; -1 0 0];
T2 = (1 / sqrt(2)) * U * T2 * V’;

% get the flattened versions
t0 = vec(T0);
t1 = vec(T1);
t2 = vec(T2);

% get the singular values in an order that is consistent
% with the numbering in [Smith et al. 2019]
s0 = Sigma (1 ,1);
s1 = Sigma (2 ,2);
s2 = Sigma (3 ,3);

H = 2 * eye(9 ,9);
H = H - (4 / (s0 + s1)) * (t0 * t0’);
H = H - (4 / (s1 + s2)) * (t1 * t1’);
H = H - (4 / (s0 + s2)) * (t2 * t2’);
end

5.3.3 Symmetric Dirichlet

\[\Psi=||F||_F^2+||F^{-1}||_F^2\\ vec(\frac{\part^2\Psi}{\part F^2})=\frac1{2I_3}g_1g_1^TT+\frac1{I_3}(\frac{I_1}2-2)H_1+(1-\frac1{4I_3})H_2+\\\frac1{(I_3)^3}(-4I_1+\frac{I_1^2-I2}2)g_3g_3^T+\frac1{(I_3)^3}(2I_1-\frac{I_1^2-I2}4)H_3 \]

Chapter 6. A Better Way for Anisotropic Solids

loose retelling of Kim et al. (2019)

6.1 What's Anisotropy?

isotropic solids: If you squish them down, they bulge out the exactly the same amount in all directions.

Most real-world solids don’t do this exactly; they can be stiffer in some directions and softer in others:

A length of hemp rope sags less, because it contains stiff fibers that all run in a preferred direction

Muscles are an anisotropic material. When fibers contract, they pull on the bones on either side of the muscle.

6.2 Invariants

we have some anisotropy direction \(a\in R^3\).

\[I_4=a^TSa\\ I_5=a^TS^TSa \]

We write a sign-extracting signum function \(S\) and its derivative \(\delta\):

\[S(x)=\left\{ \begin{aligned} &-1&if\,x<0\\ &0&if\,x=0\\ &1&if\,x>0 \end{aligned} \right.\\ \delta(x)=\left\{ \begin{aligned} &0&if\,x\ne0\\ &\infin&if\,x=0 \end{aligned} \right. \]

Hessian和PK1的公式：github的HOBAKv1仓库中src/Hyperelastic/Volume/ANISOTROPIC_XXX.cpp

Chapter 7. Thin Shell Forces

To model cloth/thin shell behavior we need to describe forces models that resist changes to the stretching, shearing and bending modes of deformation.

7.1 Handy derivatives of a few vector operators

Notation:

\(p_i\): column vector (position of \(i\)th vertex)
\(x_{ij}=p_j-p_i\): column vector
\(\hat x_{ij}={x_{ij}}/{||x_{ij}||}\): unit vector

Energy Functions, Forces and their Jacobians：

BW(1998) define an energy function \(E\) for each mode of deformation with a stiffness factor \(k_s\) to control how strong the resulting force is that resists the deformation:

\[E=\frac{k_s}{2}C(p)^TC(p) \]

in which:

\(C(p)\): the constraint function (is zero when there is no deformation in the given mode)

the restoring force (negative derivative of \(E\)):

\[f_i=-k_s\frac{\part C(p)^T}{\part p_i}C(p) \]

the second positional derivative of \(E\):

\[\frac{\part f_i}{\part p_j}=-k_s\left[\frac{\part C(p)}{\part p_i}\frac{\part C(p)^T}{\part p_j} +\frac{\part^2 C(p)}{\part p_i\part p_j}C(p)\right] \]

7.2 Strain, Energy and Energy hessian

Notation:

given a triangle with vertices \(p_i,p_j,p_k\),
- define the column vector \(\Delta p_1=p_j-p_i\) and \(\Delta p_2=p_k-p_i\).
- define the scalar \(\Delta u_1=u_j-u_i\) and \(\Delta u_2=u_k-u_i\).
- define the scalar \(\Delta v_1=v_j-v_i\) and \(\Delta v_2=v_k-v_i\).
- \(u,v\) 是点 \(p\) 的在二维板片中的坐标（MD中的衣服板片） (cloth is cut from bolts of fabric with inherent warp and weft directions，布料是横平竖直的，将布料裁剪成块之后才组装成了衣服)

We then define a linear mapping function from material space to world space \([w_u,w_v]\) in equation:

\[[w_u\quad w_v]=[\Delta p_1\quad\Delta p_2]\pmatrix{\Delta u_1&\Delta u_2\\\Delta v_1&\Delta v_2}^{-1} \]

Constraint function包含三部分:

stretch
shear
bending

三个部分都需要计算f关于位置的导数，公式请见原文

Chapter 8. Implicit Integration Methods

8.1 Backward Differentiation Methods in Fizt(Pixar)

The general \(s\)-order BDF formula(BDF: backward differentiation formula):

\[\sum_{k=0}^sa_ky_{n+k}=h\beta f(t_{n+s},y_{n+s}) \]

where:

\(s\): 使用 \(s+1\) 个时间点来求解

\(s=1\): 使用 \(x_{n+1},x_n\).

\(s=2\): 使用 \(x_{n+1},x_n,x_{n-1}\).

\(a_k,\beta\): coefficient

\(h\): time step

\(y_n=(x_n,v_n)\): state (position and velocity)

The first-order BDF form: (Backward Euler)

\[\pmatrix{x_{n+1}-x_n\\v_{n+1}-v_n}=h\pmatrix{v_{n+1}\\M^{-1}f(x_{n+1},v_{n+1})} \]

对 \(f(x_{n+1},v_{n+1})\) 使用泰勒展开：

\[f(x_{n+1},v_{n+1})\approx f(x_n,v_n)+\frac{\part f}{\part x}\Delta x+\frac{\part f}{\part v}\Delta v \]

可以获得两个线性系统：

\[(M-h\frac{\part f}{\part v}-h^2\frac{\part f}{\part x})\Delta v=h(f+h\frac{\part f}{\part x}v_n)\\ (M-h\frac{\part f}{\part v}-h^2\frac{\part f}{\part x})\Delta x=h^2 f_n+h(M-h\frac{\part f}{\part v})v_n \]

The second-order BDF form:

\[\pmatrix{\frac32x_{n+1}-2x_n+\frac12x_{n-1}\\\frac32v_{n+1}-2v_n+\frac12v_{n-1}}=h\pmatrix{v_{n+1}\\M^{-1}f(x_{n+1},v_{n+1})} \]

可以获得两个线性系统：

\[(M-\frac23h\frac{\part f}{\part v}-\frac49h^2\frac{\part f}{\part x})\Delta v=\frac13M(v_n-v_{n-1})+\frac23h\left(f(x_n,v_n)+\frac13\frac{\part f}{\part x}(x_n-x_{n-1}+2hv_n)\right)\\ (M-\frac23h\frac{\part f}{\part v}-\frac49h^2\frac{\part f}{\part x})\Delta x=\frac M3\left((x_n-x_{n-1})+\frac13h(8v_n-2v_{n-1})\right)+\frac29h\left(2hf_n-\frac{\part f}{\part v}(x_n-x_{n-1}+2hv_n)\right) \]

Chapter 9. Constrained Backward-Euler

using backward euler and solve sparse linear system \(Ax=b\).

9.1 Constraint Filtering in Baraf-Witkin Cloth

To precisely control the motion of a particle in any of its Degrees of freedom(DOFs), a mass modification scheme was developed.

Each particle has a filter \(S_i\):

\[S_i=\left\{ \begin{aligned} &I&\textrm{if zero DOFs are constrained}\\ &I-\hat p_i\hat p_i^T&\textrm{if one DOF (}\hat p_i\textrm{) is constrained}\\ &I-\hat p_i\hat p_i^T-\hat q_i\hat q_i^T&\textrm{if one DOFs (}\hat p_i,\hat q_i\textrm{) is constrained}\\ &0&\textrm{if all three DOFs are constrained} \end{aligned} \right. \]

These filters are applied within the inner loop of of the original Fizt solver. At each iteration the solution is projected onto a state which satisfies the constraints

9.1.1 Pre-filtered Preconditioned Conjugate Gradient(PPCG)

Tamstorf et al. (2015 TOG) Smoothed aggregation multigrid for cloth simulation

The formulation is:

\[\begin{aligned} (SAS^T+I-S)y&=Sc\\ y&=x-z\\ c&=b-Az \end{aligned} \]

The key equation is the first one:

The \(SAS^T\) term is a straightforward projection of the original matrix \(A\) into the subspace spanned by the filters \(S\).

The \(I-S\) term then serves to improve the conditioning of this subspace

9.2 Reverse Cuthill-McKee

In the following sections, we will discuss some techniques that we have used to improve the performance of system assembly and our PPCG solver.

Reverse Cuthill-McKee: 提高内存访问的连续性

R-C-McKee

Sparsity pattern: (a)Default mesh ordering (b) After Reverse Cuthill-McKee

Chapter 10. Collision Processing

10.1 Proximity Queries

We define a collision envelope around the surface by offset and inset distance parameters. To prevent against redundant contact data, we partition the space above and below each mesh face using this envelope. (说人话就是：在网格内外都加一层包皮，每个面附近的region称为一个cell)

2Ddetection

A 2D version of our collision cell strategy. On the left, we are looking at the red vertex, and seeing which line segments in the mesh it is in collision with. On the middle left, it is inside the blue cell, so the vertex is in collision with that segment. On the middle right, what if the mesh doubles back on itself? On the right, it is totally legal for the vertex to be inside cells of two segments, both the blue and green ones. The vertex is in collision with both segments.

3Ddetection

Particle-Kinematic Face Cell Check: (a) Cell region containing a particle in the envelope. (b) Particle is inside the cell, but outside envelope.（顶点在包皮内，则产生碰撞）

10.2 Collision Energies

10.2.1 Vertex-face energy

McAdams et al.(2011 TOG)Efficient elasticity for character skinning with contact and collisions.

In this paper, the collision energy is a spring energy defined in world-space as:

\[\Psi_{Mc}(x_p,x_s)=k(x_p-x_s)^TN(x_p-x_s) \]

where:

\(x_p\): the vertex in collision

\(x_s=\alpha x_1+\beta x_2+\gamma x_3\): the closest point on the colliding triangle

\((\alpha,\beta,\gamma)\): barycentric coordinate

\(N=nn^T\): this matrix is the outer product of the normal of the colliding triangle

On the left, the vertex-face collision configuration. In the middle, this is how we’ll label the vertices. The vertices on the face are ordered clockwise. On the right, we can use this ordering to compute a normal for the face.

An entirely equivalent energy is:

\[\Psi_{Mc}(x_p,x_s)=k(t^Tn)^2 \]

in which:

\(t=x_p-x_s=x_0-\alpha x_1-\beta x_2-\gamma x_3\)

derivative:

\[\frac{\part \Psi_{Mc}(x_p,x_s)}{\part x}=2k(t^Tn)\left[\frac{\part t^T}{\part x}n+\frac{\part n}{\part x}^Tt\right] \]

其中：

\[\frac{\part t}{\part x}=\left[I_{3\times3}\quad-\alpha I_{3\times3}\quad-\beta I_{3\times3}\quad-\gamma I_{3\times3} \right]\in R^{3\times12} \]

force:

\[f=-a\frac{\part \Psi_{Mc}(x_p,x_s)}{\part x} \]

where \(a\) denotes the area of the collision, which we set to the rest-state area of the triangle, and 1/3 the rest-state area of the all the triangles in the one-ring of the colliding vertex

Collision Hessian:

\[\frac{\part^2 \Psi_{Mc}(x_p,x_s)}{\part x^2}=2k\left[gg^T+t^Tn(t^T\frac{\part^2 n}{\part x^2}+\frac{\part t^T}{\part x}\frac{\part n}{\part x}+\frac{\part n^T}{\part x}\frac{\part t}{\part x})\right] \]

其中：

\[g=\frac{\part t^T}{\part x}n+\frac{\part n^T}{\part x}t\in R^{12} \]

10.2.2 Edge-Edge energy

energy:

\[\Psi_{e\times e}=k\cdot(\epsilon+(x_b-x_a)^Tn)^2 \]

where:

\(x_0,x_1\): endpoint of edge \(e_0=x_1-x_0\).

\(x_2,x_3\): endpoint of edge \(e_1=x_3-x_2\).

\(n=e_1\times e_0\).

\(x_a,x_b\): two closest point on \(e_0,e_1\).

\(x_a=a_0x_0+a_1x_1\).

\(x_b=b_0x_2+b_1x_3\).

gradient:

\[\frac{\part\Psi_{e\times e}}{\part x}=2k\left(\epsilon+(x_b-x_a)^Tn\right)\left[\frac{\part x_b-x_a}{\part x}n+(x_b-x_a)\frac{\part n}{\part x}\right] \]

其中：

\[\frac{\part x_b-x_a}{\part x}=[-a_0I_{3\times3}\quad-a_1I_{3\times3}\quad b_0I_{3\times3}\quad b_1I_{3\times3}] \]