Position Based Dynamics

limingyin18
421 阅读1分钟

介绍

基于位置求解。

方法

约束(Constraint)

对于NN个粒子，MM个约束：

C1(x)=0...CM(x)=0C_1(\mathbf{x})=0 \\ ... \\ C_M(\mathbf{x})=0

其中，x=[x1T,...,xNT]T\mathbf{x}=[\mathbf{x}^\mathsf{T}_1,...,\mathbf{x}^\mathsf{T}_N]^\mathsf{T}。

非线性高斯—桑德尔求解(The Non-Linear Gauss-Seidel Solver)

单独地求解每个约束C(x)C(\mathbf{x})。在一次迭代中，顶点位置的修正立即应用到下一个约束的求解。 给定x\mathbf{x}，求解Δx\Delta \mathbf{x}，使得C(x+Δx)=0C(\mathbf{x}+\Delta \mathbf{x})=0。 作线性近似：

C(x+Δx)≈C(x)+∇C(x)TΔx=0C(\mathbf{x}+\Delta \mathbf{x})\approx C(\mathbf{x}) + \nabla C(\mathbf{x})^{\mathsf{T}} \Delta \mathbf{x} = 0

由动量守恒定律可知Δx\Delta \mathbf{x}的方向为∇C\nabla C的方向，引入拉格朗日乘数子λ\lambda：

Δx=λM−1∇C(x)\Delta \mathbf{x} = \lambda \mathbf{M}^{-1}\nabla C(\mathbf{x})

其中，M=diag(m1,m2,...,mN)\mathbf{M}=diag(m_1, m_2, ..., m_N)。

解得

λ=−C(x)ΔC(x)TM−1ΔC(x)=−C(x)∑imi−1∥∇xiC(x)∥2\lambda = -\frac{C(\mathbf{x})}{\Delta C(\mathbf{x})^\mathsf{T}\mathbf{M}^{-1}\Delta C(\mathbf{x})} = -\frac{C(\mathbf{x})}{\sum_im_i^{-1}\|\nabla_{\mathbf{x}_i} C(\mathbf{x})\|^2}

对于粒子ii：

Δxi=λmi−1∇xiC(x)\Delta \mathbf{x}_i=\lambda m^{-1}_{i} \nabla_{\mathbf{x}_i}C(\mathbf{x})

指定约束

Stretching

距离约束：

C(x1,x2)=∥x1−x2∥−dC(\mathbf{x}_1,\mathbf{x}_2)=\|\mathbf{x}_1-\mathbf{x}_2\| - d

梯度：

∇C(x)=[∇x1C(x)∇x2C(x)]=[x1−x2∥x1−x2∥−x1−x2∥x1−x2∥]\nabla C(\mathbf{x}) = \left[\begin{matrix} \nabla_{\mathbf{x}_1}C(\mathbf{x}) \\ \nabla_{\mathbf{x}_2}C(\mathbf{x}) \end{matrix}\right] = \left[\begin{matrix}\frac{\mathbf{x}_1-\mathbf{x}_2}{\|\mathbf{x}_1-\mathbf{x}_2\|} \\ -\frac{\mathbf{x}_1-\mathbf{x}_2}{\|\mathbf{x}_1-\mathbf{x}_2\|}\end{matrix}\right]

拉格朗日乘数子：

λ=−∥x1−x2∥−dm1−1+m2−1\lambda = -\frac{\|\mathbf{x}_1-\mathbf{x}_2\| - d}{m_1^{-1}+m_2^{-1}}

位置修正

Δx1=−m1−1m1−1+m2−1(∥x1−x2∥−d)x1−x2∥x1−x2∥Δx2=m2−1m1−1+m2−1(∥x1−x2∥−d)x1−x2∥x1−x2∥\begin{aligned} \Delta \mathbf{x}_1 &= \frac{-m_1^{-1}}{m_1^{-1}+m_2^{-1}}(\|\mathbf{x}_1-\mathbf{x}_2\| - d)\frac{\mathbf{x}_1-\mathbf{x}_2}{\|\mathbf{x}_1-\mathbf{x}_2\|} \\ \Delta \mathbf{x}_2 &= \frac{m_2^{-1}}{m_1^{-1}+m_2^{-1}}(\|\mathbf{x}_1-\mathbf{x}_2\| - d)\frac{\mathbf{x}_1-\mathbf{x}_2}{\|\mathbf{x}_1-\mathbf{x}_2\|} \end{aligned}

算法

  • 第6行：为了连续碰撞检测，需要当前位置和预测位置。
  • 第15行：根据摩擦力系数和恢复力系数调整速度。

实现

void Stretching::solveConstraint()
{
    RigidBody &rb0 = mPBD.mRigiBodies[mBodyIndexes[0]];
    RigidBody &rb1 = mPBD.mRigiBodies[mBodyIndexes[1]];

    Vector3f x0x1 = rb0.mBaryCenter - rb1.mBaryCenter;
    float x0x1N = x0x1.norm();
    float lambda = -(x0x1N - distance) / (rb0.mInvMass + rb1.mInvMass);

    Vector3f deltaX0 = rb0.mInvMass * lambda / x0x1N * x0x1;
    Vector3f deltaX1 = -rb1.mInvMass * lambda / x0x1N * x0x1;

    rb0.mVelocity += deltaX0 / mPBD.dt;
    rb1.mVelocity += deltaX1 / mPBD.dt;
    rb0.mBaryCenter += deltaX0;
    rb1.mBaryCenter += deltaX1;
}

演示

参考

  1. M. Zwicker and C. Soler, Position-Based Simulation Methods in Computer Graphics, 2015, The Eurographics Association, 10.2312/egt.20151045
avatar
文章
阅读
粉丝