Page 309 - 《软件学报》2020年第10期

孙晓鹏 等:局部各向异性的薄壳收缩变形 3285

T
i 4 =tr(r fA),i 5 =tr(cA),
则各向异性 ARAP 能量密度场 Ψ ARAP 及 1st Piola-Kirchhoff 应力张量 ( ) R∈Pf 32× 为
() ⎞
2 ⎛ ∂Ψ Sign i
= μ /2 ( Ψ Sign i , () =Pf ARAP = μ 1− ⎜ 4 ⎟ fA ,
()) i −
ARAP 5 4 ∂f ⎜ ⎝ i 5 ⎟ ⎠
其中,Sign(·)为符号函数.进而,本文各向异性 ARAP 变形能量定义为
E = A ∫ Γ Ψ ARAP dX (4)
∂ E ∂ E ∂ E
)i ∈
记三角面片 f j 的 3 个顶点能量梯度 A ( 1,2,3 .为使模型向内凹陷,本文以新的梯度方向 A = A +k n j
A
∂x ji ∂ ji ∂x x ji
∂ E
代替 A , 其中,n j 为面 f j 的法向量,k A 为常量.
∂x ji
4 碰撞检测

针对薄壳收缩变形过程中发生的顶点与三角面片、边与边之间的自碰撞问题,本节采用连续碰撞检测,即
检测模型中任意点和面(点-面图元对),或任意两条边(边-边图元对),在某一特定的时间步长内有无自相交.2002
年 Bridson 等人 [22] 以求解 3 次多项式的方法检测点-面碰撞和边-边碰撞.考虑三维模型在某一时间步长内,任意
点和面、任意两条边都可能发生碰撞,在每一时间步长内求解 3 次多项式计算效率过低.

Fig.7 The vertex-triangle collision(left) and the edge-edge collision(right)
图 7 点-面碰撞(左)、与边-边碰撞(右)
为提高计算效率,本节首先采用轴向平行包围盒 [20] 与非渗透滤波器 [21] 等碰撞剔除算法作为预处理,剔除不
可能发生碰撞的图元对,降低误报率,然后采用 Bridson 方法检测碰撞.如图 7 图元对碰撞所示,首先定义阈值τ.
当点 i 与面 f j 的距离 d vf 小于τ时(如图 7 左图所示),点 i 将与面 f j 内的点 q j 发生碰撞,点 i 的位置为 x i ,点 q j
的位置为 x qj =α j1 x j1 +α j2 x j2 +α j3 x j3 ,其中,x j1 、x j2 、x j3 为面 f j 的 3 个顶点的位置,α j1 、α j2 、α j3 为面 f j 这 3 个顶点的
权值,且α j1 +α j2 +α j3 =1.
当边 e l 与边 e k 的距离 d ee 小于τ时(如图 7 右图所示),边 e l 上一点 q l 将与边 e k 上一点 q k 发生碰撞,点 q l 位置
为 x ql =α l1 x l1 +α l2 x l2 ,其中,x l1 和 x l2 分别为边 e l 的两个端点位置,α l1 和α l2 分别为边 e l 的两个端点权值,α l1 +α l2 =1;
点 q k 的位置为 x qk =α k1 x k1 +α k2 x k2 ,其中,x k1 和 x k2 分别为边 e k 的两个端点位置,α k1 和α k2 为边 e k 的两个端点权值,
α k1 +α k2 =1.
检测到发生碰撞的图元对后,分别定义点-面碰撞的约束函数和边-边碰撞的约束函数.若点 i 和面 f j 发生碰
撞,定义 C vf =|x i −(α j1 x j1 +α j2 x j2 +α j3 x j3 )|–δ;若边 e l 和边 e k 发生碰撞,则定义 C ee =|(α l1 x l1 +α l2 x l2 )–(α k1 x k1 +α k2 x k2 )|–δ,其
–4
中,δ=10 .
5 算法流程

本节基于第 3 节弹性变形能 E S 、弯曲变形能 E B 和局部各向异性变形能 E A 的定义,构造约束函数组,并给
出完整的局部各向异性薄壳收缩变形算法流程如下.
[5]
首先,指定各向异性变形能 E A 局部作用区域,只考虑模型所受外部压力,基于 Sympletic 欧拉法 计算每个
顶点的预测位置 p={p 0 ,p 1 ,…,p n }.