增量势能接触法(IPC)有限元模拟实现杂谈

Posted on 2024-09-19

曾经心中的大山终于被撼动了一角，那些遥远的光芒现在仿佛触手可及。

$\text{0. 闲话：为什么要写这篇文章}$

增量势能接触法，缩写为IPC，是由CMU李旻辰教授（博士时期）开发的一种全新的碰撞处理算法（当然也有几年了），相比于以往所有的碰撞处理算法，
IPC能实现任意时间步下的稳定模拟，而且模型也简单优雅。秉持着“学习一个算法就必须亲自动手“的原则，我在今年五月末从UCB回国之后，开始整理资料，
准备在自己的小框架上实现一个迷你的IPC玩具，体验整个流程。

在2023年暑假的中科大计算机图形学暑期课程中，有幸听了蒋陈凡夫教授的“Variational Contact”讲座，第一次认识了IPC算法背后的优化原理，但是依然似懂非懂。
那个时候FEM也没写对，CCD也不会写，空间加速结构也没写过，连续介质力学学得一知半解。而在将近一年的积累之后，继续整理资料，逐渐意识到自己有了可以开始尝试实现IPC的基础。

在这个过程中，我的主要参考就是2020年原始的IPC论文。然而首先还有一些前置的知识和算法需要学习，除此以外，在和其他前沿大佬的交流之中，我意识到2020年论文中的某些做法已经有了替代，对于最原始的接触模型也有了更符合物理意义的解读。
同时，原始代码写得较为庞杂混乱~~就是一坨大屎山~~，私以为是不那么易读的。在实现自己的版本的过程中，所遇到的疑问，踩到的坑，方法的选择和自己的理解都记录在此。

$\text{1. 形变，Stable Neohookean，和求导}$

这部分主要参照的内容是 $2022$ 年的 $\text{SIGGRAPH Course Dynamic Deformables}$ ，其中提出了一种比较万金油的求导方法：首先使用形变梯度 $F$ 的主不变量表示应变能密度，然后分别求导。

与传统力学使用的柯西主不变量不同，这里采用的三个主不变量是：

I_1 = \text{tr}(S), I_2 = \text{tr}(F^TF), I_3 = \text{det}(F), F = RS

可以证明，我们新采用的这一套主不变量可以表示原本柯西主不变量的所有信息。而且这三个主不变量要求导也是很容易的。

为了稳定起见，选用 $\text{Stable Neohookean}$ 应变能密度，这个没啥好说的，稳定的原因是它只使用了 $I_1$ 和 $I_3$ 这两个对翻转敏感的主不变量，并且通过添加额外的项来保证它满足赛斯-希尔应变度量的条件。

只要把 $\text{Course}$ 里的求导抄一抄，差不多就行了，但是一定要考虑张量展开为矩阵所面对的主序问题，我的 $\frac{\partial F}{\partial x}$ 的计算一开始直接抄，是错的，实际上应该是下面这样：

if constexpr (Dim == 3) {
      Matrix<T, Dim * Dim, Dim * (Dim + 1)> result;
      result.setZero();
      const auto &r0 = m_Dm_inverse.row(0);
      const auto &r1 = m_Dm_inverse.row(1);
      const auto &r2 = m_Dm_inverse.row(2);
      auto s0 = m_Dm_inverse.col(0).sum();
      auto s1 = m_Dm_inverse.col(1).sum();
      auto s2 = m_Dm_inverse.col(2).sum();
      result.col(0).reshaped(3, 3).row(0) = Vector<T, Dim>(-s0, -s1, -s2);
      result.col(1).reshaped(3, 3).row(1) = Vector<T, Dim>(-s0, -s1, -s2);
      result.col(2).reshaped(3, 3).row(2) = Vector<T, Dim>(-s0, -s1, -s2);
      result.col(3).reshaped(3, 3).row(0) = r0;
      result.col(6).reshaped(3, 3).row(0) = r1;
      result.col(9).reshaped(3, 3).row(0) = r2;
      result.col(4).reshaped(3, 3).row(1) = r0;
      result.col(7).reshaped(3, 3).row(1) = r1;
      result.col(10).reshaped(3, 3).row(1) = r2;
      result.col(5).reshaped(3, 3).row(2) = r0;
      result.col(8).reshaped(3, 3).row(2) = r1;
      result.col(11).reshaped(3, 3).row(2) = r2;
      return result;
    }

$\text{2. 连续碰撞检测(CCD)方法的选择与实现}$

这是我第一次写 $\text{CCD}$ ，我一开始考虑的两种方法都来自于 $\text{2022}$ 年的 $\text{SIGGRAPH Course Friction and Contact}$ 。这里先解释一下这两种方法，然后我们分别解释它们为什么不行。

$\text{2.1. 一种很传统的方法}$

第一种方法是首先求解四点共面的时间，这一步需要求解一个一元三次方程。然后对于所有可能的共面时间，再求解图元是否有相交。因为线段和三角形都可以参数化表示，我们只需要求一个很小的线性方程组解出参数，然后检查参数的范围是否符合预期。

要稳健地求三次方程本身就够难写了，我参考了 $\text{SIGGRAPH 2022}$ 上 $\text{Cem Yuksel}$ 在 $\text{A Fast \& Robust Solution for Cubic \& Higher-Order Polynomials}$ 里提出的方法，对于三次方程简单来说，就是从二阶导零点那里分开，对两部分分别牛顿迭代。

static Real robustCubicNewton(const CubicPolynomial &poly, Real x, Real l,
                              Real r, Real tolerance) {
  const auto &[a, b, c, d] = poly;
  while (true) {
    Real f = evalCubic(poly, x);
    Real df = evalQuadratic({.a = 3 * a, .b = 2 * b, .c = c}, x);
    if (df == 0.0) [[unlikely]]
      return x;
    Real dx = -f / df;
    x = core::clamp(x + dx, l, r);
    if (std::abs(dx) < tolerance)
      return x;
  }
}

// returns all the real root within the interval [l, r]
CubicEquationRoots clampedCubicSolve(const CubicPolynomial &poly, Real l, Real r,
                                     Real tolerance) {
  const auto &[a, b, c, d] = poly;
  if (abs(a) <= 1e-10) {
    auto result = quadraticSolve({.a = b, .b = c, .c = d});
    if (result.infiniteSolutions)
      return infiniteCubicRoots();
    else if (result.numRoots == 0)
      return {};
    else if (result.numRoots == 1) {
      if (result.roots[0] >= l && result.roots[0] <= r)
        return CubicEquationRoots(result.roots[0]);
      return {};
    }
    Real x1 = result.roots[0];
    Real x2 = result.roots[1];
    CubicEquationRoots roots;
    if (x1 > x2)
      std::swap(x1, x2);
    if (l <= x1 && x1 <= r)
      roots.addRoot(x1);
    if (l <= x2 && x2 <= r)
      roots.addRoot(x2);
    return roots;
  }
  Real x1, x2;
  Real fl = evalCubic(poly, l);
  Real fr = evalCubic(poly, r);
  Real xc = -b / (3.0 * a);
  Real fxc = evalCubic(poly, xc);
  auto derivativeRoots = quadraticSolve({.a = 3 * a, .b = 2 * b, .c = c});
  assert(!derivativeRoots.infiniteSolutions);
  if (!derivativeRoots.numRoots) {
    if (fl * fr > 0.0)
      return {};
    if (xc > l && xc < r) {
      if (fl * fxc <= 0.0)
        return CubicEquationRoots(robustCubicNewton(poly, l, l, xc, tolerance));
      return CubicEquationRoots(robustCubicNewton(poly, r, xc, r, tolerance));
    }
    return CubicEquationRoots(robustCubicNewton(poly, (l + r) * 0.5, l, r, tolerance));
  }
  x1 = derivativeRoots.roots[0];
  x2 = derivativeRoots.roots[1];
  if (x1 > x2)
    std::swap(x1, x2);
  Real fx1 = evalCubic(poly, x1);
  CubicEquationRoots roots{};
  if (l < x1) {
    Real bound = std::min(r, x1);
    Real fbound = bound == r ? fr : fx1;
    if (fl * fbound <= 0.0)
      roots.addRoot(robustCubicNewton(poly, l, l, bound, tolerance));
  }
  Real fx2 = evalCubic(poly, x2);
  Real lbound = std::max(l, x1);
  Real flbound = lbound == l ? fl : fx1;
  Real rbound = std::min(r, x2);
  Real frbound = rbound == r ? fr : fx2;
  if (lbound < rbound && flbound * frbound <= 0.0) {
    Real root;
    if (xc > lbound && xc < rbound) {
      if (flbound * fxc <= 0.0)
        root = robustCubicNewton(poly, (lbound + xc) * 0.5, lbound, xc, tolerance);
      else root = robustCubicNewton(poly, (rbound + xc) * 0.5, xc, rbound, tolerance);
    } else
      root = robustCubicNewton(poly, (lbound + rbound) * 0.5, lbound,
                               rbound, tolerance);
    roots.addRoot(root);
  }
  if (x2 < r) {
    Real bound = std::max(l, x2);
    Real fbound = bound == l ? fl : fx2;
    if (fbound * fr <= 0.0)
      roots.addRoot(robustCubicNewton(poly, r, bound, r, tolerance));
  }
  return roots;
}

应该说这个方法相当不错，基本上除了无穷多解，各种条件下都能稳健地求解，有无穷多解的时候特判一下就行。但是问题就出在无穷多解这里。

如果三次方程有无穷多解（也就是多项式恒为0），那么对应到实际物理情景，就应该是四点不仅在初始状态共面，而且它们的速度也共面。这个时候就必须同时求解关于 $t$ 和参数的方程组，来看图元在什么时候会相交。问题就在于这个方程组不是线性的，而是包括 $t$ 和参数的乘积项，这就不好弄了。

所以，在 $\text{Friction and Contact}$ 中给出了另外一个方法，同时求解关于 $t$ 和参数的多元非线性方程组的 $\text{Tight-Inclusion Based CCD}$ 。

$\text{2.2. TCCD真的那么好吗？}$

应该说我并没有完全搞明白 $\text{TCCD}$ 的原理，但是一个简单的问题摆在面前：这个方法必须维护一个区间的优先队列。这就很不可忍受了。对于 $\text{CPU}$ 可能还算好，但是对于 $\text{GPU}$ 这就太过分了。 $\text{GPU}$ 不喜欢花里胡哨的东西，后面 $\text{IPC}$ 都上 $\text{GPU}$ ，是自然不会用这个方法的。

$\text{2.3. ACCD，简单的就是好的}$

最后在大佬的建议下，选择了实现简单的 $\text{Additive CCD}$ 方法，它的想法很简单，沿着当前速度方向，走一个保证无穿插的小步长，直到距离小于阈值。

这个方法是非常契合 $\text{IPC}$ 的需求的，毕竟 $\text{IPC}$ 并不是真的关心什么时候碰了，而只是关心“我能沿当前方向安全地走多远”，有点儿预留也是OK的，而且这个方法实现很简单， $\text{GPU}$ 很喜欢。

std::optional<Real> CollisionDetector::runACCD(CCDMode mode,
                                               std::array<Vector<Real, 3>, 4> &x,
                                               std::array<Vector<Real, 3>, 4> &p,
                                               Real toi) const {
  auto computeDistance = [&]() -> Real {
    if (mode == CCDMode::EE)
      return std::sqrt(distanceSqrEdgeEdge(x[0], x[1], x[2], x[3]));
    else
      return std::sqrt(distanceSqrPointTriangle(x[0], x[1], x[2], x[3]));
  };
  Real lp = 0.0;
  if (mode == CCDMode::EE)
    lp = std::max(p[0].norm(), p[1].norm()) + std::max(p[2].norm(), p[3].norm());
  else lp = p[0].norm() + std::max(p[1].norm(), std::max(p[2].norm(), p[3].norm()));
  if (lp == 0.0)
    return std::nullopt;
  Real dis = computeDistance();
  Real g = s * dis;
  Real t = 0.0;
  Real tl = (1.0 - s) * (dis / lp);
  while (true) {
    for (int i = 0; i < 4; i++)
      x[i] += p[i] * tl;
    dis = computeDistance();
    if (dis < g + 1e-10) {
      [[unlikely]] if (t == 0.0) t = tl;
      break;
    }
    t += tl;
    if (t > toi)
      return std::nullopt;
    tl = 0.9 * dis / lp;
  }
  return t;
}

$\text{2.4. 全局的CCD}$

全局 $\text{CCD}$ ，或者说 $\text{CCD}$ 的宽阶段，就比较好说了，考虑到整体占用空间不确定，我个人用空间哈希用得不是很多，我更倾向于使用 $\text{LBVH}$ 。

把所有三角形的路径包围盒建出来做成 $text{LBVH}$ ，然后对每个三角形跑一遍 $\text{LBVH}$ 就行，跑到底了就对两个三角形做 $\text{CCD}$ ，这个就是很常规的对两个三角形所有可能的顶点和边做 $\text{CCD}$ 。

一开始我把所有的窄阶段询问都放在一个容器里，然后再统一并行地跑窄阶段，但是发现这个规模实在太过于巨大了。所以最后才改成了跑到底就做一次 $\text{CCD}$ 。

std::optional<Real> CollisionDetector::detect(const System &system, const VecXd &p) {
  m_bvh->update(SystemTriangleTrajectoryAccessor{system, p});
  std::atomic<Real> toi = 1.0;
  std::atomic<bool> hasContact{false};
  tbb::parallel_for (0, system.numTriangles(), [&](int i) {
    auto trajectory_bbox = computeTrajectoryBBox({.system = system, .p = p, .toi = toi, .idx = i});
    m_bvh->runSpatialQuery(
        [&](int primitiveID) -> bool {
          if (i >= primitiveID) return false;
          if (system.checkTriangleAdjacent(i, primitiveID)) return false;
          auto t = trianglePairCCD({.system = system, .p = p, .ta = i, .tb = primitiveID, .toi = toi});
          if (!t) return false;
          if (*t > toi) return false;
          hasContact.store(true);
          toi.store(*t);
          return true;
        },
        [&](const BBox<Real, 3> &bbox) -> bool {
          return trajectory_bbox.overlap(bbox);
        });
  });
  if (hasContact) return toi;
  return std::nullopt;
}

当然，这里的剪枝其实还可以更严格一点，比如如果最大的idx也小于i那么整个子树也不用跑了。

在最原始的 $\text{IPC}$ 中，似乎还有一步是计算一个类似于 $\text{CFL}$ 数的剪枝步长估计，用于加速整个碰撞检测。
据大佬所言，现在碰撞检测已经不是性能瓶颈了，而且这个剪枝在 $\text{GPU}$ 上也不好做，因此也没人这么做了。

$\text{3. 约束集的更新}$

众所周知， $\text{IPC}$ 里主要有边-边约束和点-三角形约束，虽然所有可能的约束有很多，
但是考虑到障碍函数的作用范围 $\hat{d}$ 很小，实际上真正“活跃”的约束并不多，因此 $\text{IPC}$ 会专门维护一个活跃的约束集。
按照原论文的伪代码，这个约束集在每次系统的构型发生变化时都要重新计算，这是好理解的。

但是我和大佬交流之后，大佬指出这实际上是没有必要的。在线搜索之前，我们首先会计算出一个“绝对无穿插”的步长，此时我们可以预计算这个步长之内所有潜在的约束。在线搜索的时候，就无需再更新，只需直接对于所有潜在约束计算即可。
这样做不仅速度更快，而且也能保证约束集的一致性。

void IpcImplicitEuler::computeEdgeEdgeConstraints(const ConstraintSetPrecomputeRequest &config) {
  const auto &[p, toi, dHat] = config;
  const auto &x = system().currentConfig();
  EdgeTrajectoryAccessor edge_accessor(system(), p, toi);
  edgesBVH->update(edge_accessor);
  for (int edge_index = 0; edge_index < system().edges().cols(); edge_index++) {
    auto q_edge_pa = x.segment<3>(system().edges()(0, edge_index) * 3);
    auto q_edge_pb = x.segment<3>(system().edges()(1, edge_index) * 3);
    BBox<Real, 3> edge_trajectory_bbox = BBox<Real, 3>({q_edge_pa(0), q_edge_pa(1), q_edge_pa(2)})
        .expand({q_edge_pb(0), q_edge_pb(1), q_edge_pb(2)});
    edgesBVH->runSpatialQuery(
        [&](int edge_idx) -> bool {
          if (system().checkEdgeAdjacent(edge_index, edge_idx))
            return false;
          constraintSet.eeConstraints.push_back(
              {.system = system(),
                  .ia = edge_index,
                  .ib = edge_idx});
          return true;
        }, [&](const BBox<Real, 3> &bbox) -> bool {
          return edge_trajectory_bbox.overlap(bbox);
        });
  }
}

void IpcImplicitEuler::computeVertexTriangleConstraints(const ConstraintSetPrecomputeRequest &config) {
  const auto &[p, toi, dHat] = config;
  trianglesBVH->update(ipc::SystemTriangleTrajectoryAccessor(system(), p, toi));
  for (int i = 0; i < system().numVertices(); i++) {
    BBox<Real, 3> vertex_trajectory_bbox = computeVertexTrajectoryBBox({system(), p, toi, i}).dilate(dHat);
    trianglesBVH->runSpatialQuery(
        [&](int triangle_idx) -> bool {
          for (int j = 0; j < 3; j++)
            if (system().surfaces()(j, triangle_idx) == i)
              return false;
          constraintSet.vtConstraints.push_back(
              {.system = system(),
                  .iv = i,
                  .it = triangle_idx});
          return true;
        }, [&](const BBox<Real, 3> &bbox) -> bool {
          return vertex_trajectory_bbox.overlap(bbox);
        });
  }
}

void IpcImplicitEuler::precomputeConstraintSet(const ConstraintSetPrecomputeRequest &config) {
  constraintSet.vtConstraints.clear();
  constraintSet.eeConstraints.clear();
  computeVertexTriangleConstraints(config);
  computeEdgeEdgeConstraints(config);
}

$\text{4. 自适应接触刚度更新？接触本构模型！}$

在原论文中，做完一步迭代以后，会根据当前构型更新接触刚度，也就是论文里的 $\kappa$ 。

但是跟大佬讨论过后，这一步实际上是早期 $\text{IPC}$ 理解不足导致的。这样的自适应刚性虽然能完成模拟，但是实际上并没有物理意义。

$\kappa$ 是否有真的物理意义呢？

我们知道， $\text{IPC}$ 中，任意的两个表面实际上没有真的接触，“障碍能量”保证的是两个表面在足够接近的时候会受到与 $\kappa$ 成正比的巨大斥力。
这实际上就是一种本构模型：相当于我们把表面的薄薄的一层“空气”看成是一种特殊的超线性材料， $\text{kappa}$ 就是这种材料的“刚度”。

因此，我们实际上可以把 $\kappa$ 设计为一个固定值。这样我们就有了一个具有物理意义的接触本构模型。

然后，这里需要指出原论文里一个（在现在看来）有问题的地方，原论文对于障碍增广增量势能的表达式是：

E = ... + h^2\Psi, B = E + \kappa \sum\limits_{c\in C} b(c)

现在想想，障碍能量能和 $\Psi$ 分开吗？实际上是不能的，它们应该一起乘上 $h^2$ 。

所以应该是：

B = ... + h^2(\Psi + \kappa \sum\limits_{c\in C} b(c))

$\text{5. 边-边障碍能量的修饰子}$

当边平行的时候，如果我们还计算边边距离，并据此计算梯度，这是不可导的。我们希望能够在两条边接近平行时，将距离的类型退化为点-边距离。在这个过程中，我们希望边边距离引起的能量平滑地趋于零。

因为如此，我们才需要对边边距离乘上一个修饰子。

$\text{6. ipc-toolkit的距离函数}$

这里需要指出的是，虽然ipc-toolkit所有求距离的函数都写着distance，但它们实际上都求的是距离平方，而且原论文里也是提到了，为了尽可能避免开方，都使用距离平方计算。

$\text{7. 后话}$

到此算是初步走完了大部分的 $\text{IPC}$ 流程，也彻底了解了没能完整实现 $\text{FEM}$ 和 $\text{CCD}$ 的遗憾。接下来还有带耗散的 $\text{IPC}$ ， $\text{Codim-IPC}$ ， $\text{GPU-IPC}$ ，这些东西仿佛都不是那么遥远了。

0. 闲话：为什么要写这篇文章\text{0. 闲话：为什么要写这篇文章}0. 闲话：为什么要写这篇文章

1. 形变，Stable Neohookean，和求导\text{1. 形变，Stable Neohookean，和求导}1. 形变，Stable Neohookean，和求导

2. 连续碰撞检测(CCD)方法的选择与实现\text{2. 连续碰撞检测(CCD)方法的选择与实现}2. 连续碰撞检测(CCD)方法的选择与实现

2.1. 一种很传统的方法\text{2.1. 一种很传统的方法}2.1. 一种很传统的方法

2.2. TCCD真的那么好吗？\text{2.2. TCCD真的那么好吗？}2.2. TCCD真的那么好吗？

2.3. ACCD，简单的就是好的\text{2.3. ACCD，简单的就是好的}2.3. ACCD，简单的就是好的

2.4. 全局的CCD\text{2.4. 全局的CCD}2.4. 全局的CCD

3. 约束集的更新\text{3. 约束集的更新}3. 约束集的更新

4. 自适应接触刚度更新？接触本构模型！\text{4. 自适应接触刚度更新？接触本构模型！}4. 自适应接触刚度更新？接触本构模型！

5. 边-边障碍能量的修饰子\text{5. 边-边障碍能量的修饰子}5. 边-边障碍能量的修饰子

6. ipc-toolkit的距离函数\text{6. ipc-toolkit的距离函数}6. ipc-toolkit的距离函数

7. 后话\text{7. 后话}7. 后话