(2014) Projective Dynamics论文摘要

本文禁止转载
B站：Heskey0

---

Projective Dynamics: Fusing Constraint Projections for Fast Simulation

Abstract

Our approach builds a bridge between nodal Finite Element methods and Position Based Dynamics. Inspired by continuum mechanics(连续介质力学), we derive a set of continuum-based potentials that can be efficiently incorporated within our solver.

Chapter 1. Introduction

​ The desirable qualities of PBD come at the cost of limited accuracy, because PBD is not rigorously derived from continuum mechanical principles.

​ The main advantage of our constraint-based potentials is that their structure enables an efficient local/global optimization (block coordinate descent)

• The local step consists of projecting every element onto the constraint manifold (solving a small nonlinear problem per element).
  • The local steps consists of small independent optimization problems, which can be all executed in parallel.
• The global step combines the results of individuals projections, finding a compromise between all of the individual constraints, while also taking into account global effects such as inertia and external forces

The local/global approach allows us to formulate an implicit integration solver

Chapter 2. Related Work

​ We also employ local/global alternation in our approach, but contrary to [Liu et al. 2013](fast mass spring system), which is limited to mass-spring systems and assumes only linear springs (Hooke’s law), we show how to generalize this concept employing projection onto constraint sets to simulate general nodal dynamical systems.

Chapter 3. Continuum Mechanics View

We start with the implicit time integration of FEM-discretized elastic models

3.1 Implicit Euler Solver

\[\min_{q_{n+1}}\frac1{2h^2}||M^{\frac12}(q_{n+1}-s_n)||_F^2+\sum_iW_i(q_{n+1}) \]

where

• \(f_{ext}\): The sum of the external forces.
• \(s_n=q_n+hv_n+h^2M^{-1}f_{ext}\).
• \(W_i(q)\): a scalar potential energy function
• \(||\cdot||_F\): the Frobenius norm

3.2 Nonlinear Elasticity

Numerous elastic potentials used in practice are formulated as a function of the strain using a (often nonlinear) material model \(\Psi(\cdot)\), resulting in elastic potentials \(W(q)=\Psi(E(q))\).

• From a geometric point of view, we can observe that \(E(q) = 0\) defines a constraint manifold of all possible undeformed configurations, while \(Ψ(E(q))\) measures how far the deformed configuration is from this manifold.

• the function \(\Psi(E(\cdot))\) defines both the constraint manifold \(E(·) = 0\) as its zero level set and the elastic potential given by its isolines. By introducing a projection variable \(p\) in the manifold, we can decouple the manifold definition from the elastic potential, modeled as the distance function \(d(q, p)\).

The simplified potentials:

\[W(q,p)=\frac w2||Aq-Bp||_F^2+\delta_C(p) \]

where:

• \(p\): the auxiliary variable.
• \(A,B\): constant matrices.
• \(w\): nonnegative weight.
• \(\delta_C(p)\): an indicator function that evaluates to \(0\) if \(C(p)=0\) and to \(+\infin\) otherwise.

3.3 Projective Implicit Euler Solver

minimize:

\[\frac1{2h^2}||M^{\frac12}(q-s_n)||_F^2+\sum_i\left[\frac{w_i}{2}||A_iS_iq-B_ip_i||_F^2+\delta_{C_i}(p_i)\right] \]

where:

• \(S_i\): a constant selection matrix that selects the vertices involved in the \(i\)th constraint.
• \(A_i=B_i=M_i^{1/2}\).

(1) Local solve

最小化每一个小约束：

\[\min_{p_i}W_i(p_i) \]

即：

\[\min_{p_i}\left[\frac{w_i}2||A_iS_iq-B_ip_i||_F^2+\delta_{C_i}(p_i)\right] \]

(2) Global solve

解线性系统：

\[\left(\frac{M}{h^2}+\sum_iw_iS^T_iA_i^TA_iS_i\right)q=\frac{M}{h^2}s_n+\sum_iw_iS^T_iA^T_iB_ip_i \]

选用\(A_i=B_i=M_i^{1/2}\)之后：

\[\left(\frac{M}{h^2}+\sum_iw_iS^T_iMS_i\right)q=\frac{M}{h^2}s_n+\sum_iw_iS^T_iMp_i \]

The system matrix is constant as long as the constraints are not changing and therefore can be prefactored at initialization, allowing for very efficient global solves

Chapter 4. Position Based Dynamics View

PBD和implicit Euler是相似的，区别在于PBD将优化目标中的势能项(elastic potential)换成了约束

4.1 Gauss-Seidel Solver

classical PBD performs three steps:

1. the positions are initialized by an explicit Euler step, ignoring internal forces
2. the positions are updated by projecting the current configuration consecutively on each constraint set respecting the mass weighting.
3. the velocities are updated as \(v_{n+1}=(q_{n+1}-q_{n})/h\).

the constraint resolution strategy of PBD actually implements a Gauss-Seidel type minimization on the energy:

\[\frac12\sum_i||M^{1/2}_i(S_iq-p_i)||^2_F+\delta_{C_i}(p_i) \]

where:

• \(M_i\): a lumped mass matrix \(M_i\) only involving the constraint's points.

the Gauss-Seidel approach minimizes this energy by optimizing each summand sequentially. (minimizing potentials of the form \(\frac12||M^{1/2}\Delta q||_F^2+\delta_C(q+\Delta q)\) where we introduce corrections \(\Delta q=p-q\) to simplify the derivation) Using Lagrange multipliers for the linearized constraint \(C(q)+tr(\nabla C(q)^T\Delta q)=0\), we define the Lagrangian

\[\frac12||M^{1/2}\Delta q||^2_F+\lambda(C(q)+tr(\nabla C(q)^T\Delta q)) \]

Using the critical point condition, we find the optimal direction \(\Delta q=-\lambda M^{-1}\nabla C(q)\). The Lagrange multiplier \(\lambda\) can then be found by requiring that the linearized constraint vanishes in this direction, i.e., \(C(q)-\lambda||M^{-1/2}\nabla C(q)||_F^2=0\), leading to the final update:

\[\Delta q=-M^{-1}\nabla C(q)\frac{C(q)}{||M^{-1/2}\nabla C(q)||_F^2} \]

Chapter 5. Continuum-Based Constraints

本节介绍一些continuous energies以及它们的discretizations \(W\)

5.1 Strain

discretized strain energy:

\[W(q,T)=\frac w2A||X_fX_g^{-1}-T||_F^2+\delta_M(T) \]

where:

• \(A\): triangle area
• \(T\): desired point-wise transformations
• \(X_f=[q_j-q_i,q_k-q_i]\in R^{2\times2}\)
  • \(X_f\) contains the triangle edges of the current configuration isometrically contains the triangle edges of the current configuration isometrically
  • \(X_g\) contains the triangle edges of the rest configuration

5.2 Area and Volume Preservation

Area and volume preservation is important for simulating incompressible materials.

\[\sigma_{\min}<\det(T)<\sigma_{\max} \]

If \(\sigma_{\min}<1\): the modeled material allows for compression.

If \(\sigma_{\max}>1\): the material allows for expansion.

5.3 Example-Based

Example-based simulation allows modeling artistic elastic material behavior by supplying a few deformation examples that the material should follow.

\[W(q,R,w)=\frac w2V||X_fX_g^{-3}-RX_h(w)X_g^{-1}||_F^2+\delta_{SO(3)}(R) \]

where:

• \(X_h(2)=X_g+\sum_iw_i(R_iX_{g_i}-X_g)\)
• \(w\): can either be defined locally per element or globally, resulting in local or global coupling of the deformation
• \(R_i\): precomputed rotation matrices
• \(g_i\): define the piecewise linear coordinate functions of the examples

5.4 Bending

\[W(q,R)=\frac w2A||X_fc-RX_gc||_2^2+\delta_{SO(3)}(R) \]

where:

• \(A\): the Voronoi area of the vertex
• \(X_f,X_g\): contain the one-ring edges of the vertex for the current configuration and for the rest configuration
• \(c\): stores the common cotangent weights divided by the Voronoi area.

Appendix: A Local Solver

Strain

When minimizing over \(T\) while keeping \(q\) fixed in the local step

\[\min_T||X_fX_g^{-1}-T||_F^2+\delta_M(T) \]

the optimization can be reformulated as

\[\min_{\Sigma^\star}||\Sigma-\Sigma^\star||_F^2\quad s.t.\quad\sigma_{\min}\le\Sigma_{ii}^\star\le\sigma_{\max} \]

where \(X_fX_g^{-1}=U\Sigma V^T\) and \(T=U\Sigma^\star V^T\). The optimal solution can be computed as \(\Sigma^\star\) being the singular values \(\Sigma\) clamped between \(\sigma_{\min}\) and \(\sigma_{\max}\). For tetrahedrons, if \(\det(X_fX_g^{-1})<0\), the last singular value is negated to avoid reflections.

Area and Volume

Similar to the strain constraint the local minimization of the volume constraint can be reformulated as:

\[\min_{\Sigma^\star}||\Sigma-\Sigma^\star||^2_F\quad s.t.\quad \sigma_{\min}\le\Pi_i\Sigma_{ii}^\star\le\sigma_{\max} \]

This problem can be further transformed in

\[\min_D||D||_2^2\quad s.t.\quad\Pi_i(\Sigma_{ii}+D_i)=\sigma \]

with \(\Sigma_{ii}^\star=\Sigma_{ii}+D_i\) and where \(\sigma=\sigma_{\min}\) when \(\Pi_i\Sigma_{ii}^\star<\sigma_{\min}\) and \(\sigma=\sigma_\max\) when \(\Pi_i\Sigma_{ii}^\star>\sigma_{\max}\). This constrained minimization can be solved by iteratively solving a quadratic programming problem by linearizing the constraint leading to a simple update rule:

\[D^{k+1}=\frac{\nabla C(D^k)^TD^k-C(D^k)}{||\nabla C(D^k)||_2^2}\nabla C(D^k) \]

where:

• \(C(D)=\Pi_i(\Sigma_{ii}+D_i)-\sigma\).

Example-Based

We solve the optimization

\[\min_{R,w}||X_fX_g^{-1}-RX_h(w)X_g^{-1}||^2_F+\delta_{SO(3)}(R) \]

using a local/global approach by minimizing over \(R\) and \(w\) iteratively. The minimization over R is solved using SVD following and solving over w corresponds to solve a simple linear system.

Bending

The local solve of the bending constraint can be formulated as

\[\min_R||v_f-Rv_g||_2^2+\delta_{SO(3)}(R) \]

where \(v_f=X_fc\) and \(v_g=X_gc\). This corresponds in finding a rotation \(R\) such that the rotated vector \(v_g\) matches best the vector \(v_f\) . While \(R\) could be found using SVD this problem has an easier closed form solution where \(Rv_g\) can be replaced by \(\frac{v_f||v_g||_2}{||v_f||_2}\).