Bending Constraint
처음부터 모든 코드를 보면 헷갈리니까 조금씩 잘라서 보자.
/**
* Bending contraint class in which all bending constraints
* are derived from
*/
abstract class BendingConstraint extends Constraint {
bendingIds: Int32Array;
constructor(
positions: Float32Array,
invMass: Float32Array,
indices: Uint16Array,
neighbors: Float32Array,
compliance: number
) {
super(positions, invMass, indices, neighbors, compliance);
this.bendingIds = new Int32Array(this.getTriPairIds());
}
...bending constraint는 이름에서 알 수 있듯 옷감 등의 모델이 구부러지는 현상에 대한 constraint이다.
주석에서 derived라고 말한 것으로 보아 구부러짐에 대해서는 미분을 이용한 수식을 이용하다보다.
constraint에서 필요한 요소들을 초기화하고, 특별히 bendingid를 삼각형의 한 pair를 이루는 vertex의 id들을 설정해준다.
그림에서 id0 id1을 잇는 점선 부분이 구부러진다고 생각하자
abstract class BendingConstraint extends Constraint {
...
private getTriPairIds(): number[] {
const numTris = this.indices.length / 3; // Every 3 vertices is a triangle
const triPairIds = [];
for (let i = 0; i < numTris; i++) {
// triangles
for (let j = 0; j < 3; j++) {
// edges
// This is one edge of a triangle id0 ------- id1
const id0 = this.indices[3 * i + j];
const id1 = this.indices[3 * i + ((j + 1) % 3)];
// Check to see if there is a neighbor triangle
// See findTriNeighbors for details
const n = this.neighbors[3 * i + j];
// Neighbor found!
if (n >= 0) {
// Need to find opposite particle ids that are on opposite sides of shared edge
// Find the last vertice in this triangle
// this is the vertice of the triangle not on the shared edge.
const id2 = this.indices[3 * i + ((j + 2) % 3)];
// Neighbor triangle (using n, since that's the shared edge of the neighbor triangle)
const ni = Math.floor(n / 3); // The neighbot triangle
const nj = n % 3; // LOCAL edge, of the neighbor triangle. (so either 0, 1, 2)
// Similar to above, find the non-shared vertice
const id3 = this.indices[3 * ni + ((nj + 2) % 3)];
triPairIds.push(id0);
triPairIds.push(id1);
triPairIds.push(id2);
triPairIds.push(id3);
}
}
}
return triPairIds;
}
...Bending Constraint에는 solve가 없고, getTriPairIds는 한 변이 맞닿는 두 삼각형 pair vertex id 4개를 짝지어서 triPairIds에 정리해넣는 함수이다.
Bending Constraint는 추상클래스이기 때문에 이를 extends한 Performant Bending Constraint, Isometric Bending Constraint에 solve가 있고, 여기서 어떻게 시뮬레이션을 처리하는지 알아보도록 하자.
Performant Bending Constraint
initBendingLengths
export class PerformantBendingConstraint extends BendingConstraint {
...
constructor(
...
this.initBendingLengths();
}
...
private initBendingLengths() {
// Calculate and initialize rest lengths of bending constraints
for (let i = 0; i < this.bendingLengths.length; i++) {
// we know id2 and id3 in bendingIds are the vertices that we want
// to add distance constraints
// see getTriPairIds for details
const id0 = this.bendingIds[4 * i + 2];
const id1 = this.bendingIds[4 * i + 3];
this.bendingLengths[i] = Math.sqrt(
vecDistSquared(this.positions, id0, this.positions, id1)
);
}
}performent bending constraint의 경우 생성자에서 initBendingLengths라는 함수를 실행한다.
함수 안에서 id0과 id1을 정의하는데, 들어가는 인덱스가 (i+2) (i+3)인 걸로 보아
private getTriPairIds(): number[] {
...
triPairIds.push(id0);
triPairIds.push(id1);
triPairIds.push(id2);
triPairIds.push(id3);triPairIds에서 id2와 id3을 의미할 것이다.
즉 구부러지는 선을 기준으로 양쪽에 있는 vertex이다.
private initBendingLengths() {
// Calculate and initialize rest lengths of bending constraints
for (let i = 0; i < this.bendingLengths.length; i++) {
// we know id2 and id3 in bendingIds are the vertices that we want
// to add distance constraints
// see getTriPairIds for details
const id0 = this.bendingIds[4 * i + 2];
const id1 = this.bendingIds[4 * i + 3];
this.bendingLengths[i] = Math.sqrt(
vecDistSquared(this.positions, id0, this.positions, id1)
);
}
}그리고 그 두 vertex id2, id3의 길이를 구한다.
즉 initBendingLength에서는 bending되기 전 길이를 초기화한다.
solve
export class PerformantBendingConstraint extends BendingConstraint {
...
solve(dt: number) {
const alpha = this.compliance / dt / dt;
for (let i = 0; i < this.bendingLengths.length; i++) {
const id2 = this.bendingIds[4 * i + 2];
const id3 = this.bendingIds[4 * i + 3];
const w0 = this.invMass[id2];
const w1 = this.invMass[id3];
const w = w0 + w1;
if (w == 0.0) continue;
vecSetDiff(this.grads, 0, this.positions, id2, this.positions, id3);
const len = Math.sqrt(vecLengthSquared(this.grads, 0));
if (len == 0.0) continue;
vecScale(this.grads, 0, 1.0 / len);
const restLen = this.bendingLengths[i];
const C = len - restLen;
const s = -C / (w + alpha);
vecAdd(this.positions, id2, this.grads, 0, s * w0);
vecAdd(this.positions, id3, this.grads, 0, -s * w1);
}
}performent bending constraint의 solve 함수에서는 어떤 물리적 처리를 할까
Distance Constraint와 다른 점은 이번에는 수직 거리가 아니다
bending이 이루어지는 두 다른 삼각형의 점 간의 거리를 이동시키는 solve이다
그래서 edge length 대신 bending length를 사용하고 나머지는 비슷하다
vecSetDiff(this.grads, 0, this.positions, id2, this.positions, id3);
const len = Math.sqrt(vecLengthSquared(this.grads, 0));
if (len == 0.0) continue;
vecScale(this.grads, 0, 1.0 / len);
const restLen = this.bendingLengths[i];
const C = len - restLen;
// normalizing factor를 이용하지 않음
const s = -C / (w + alpha);
vecAdd(this.positions, id2, this.grads, 0, s * w0);
vecAdd(this.positions, id3, this.grads, 0, -s * w1);distance constraint와 다른 점은
delta lagrangian multiplier를 구할 때 normalizing factor를 곱해 정규화를 하지 않는다는 점
그 외에는 비슷하다
간단히 보면,
정규화를 하지 않는 것은 곡률 측면에서 계산을 이루기 때문에
방향성을 정규화로 없앨 필요가 없는 것이다.
Isometric Bending Constraint
ismoetric bending constraint에서는 새롭게 Q matrix에 대한 개념이 추가된다.
Q matrix는 bending의 최소단위가 되는 네 개의 파티클이 두 삼각형을 이루는 구조에서 사용된다.
export class IsometricBendingConstraint extends BendingConstraint {
Q: Float32Array;
constructor(
positions: Float32Array,
invMass: Float32Array,
indices: Uint16Array,
neighbors: Float32Array,
compliance: number
) {
super(positions, invMass, indices, neighbors, compliance);
this.Q = new Float32Array((16 * this.bendingIds.length) / 4);
this.initQ();
}
...isometric bending constraint의 class 생성자를 보면 Q는 4*4 matrix를 즉 16개 요소가 한 그룹이므로 16을 bending id에 곱한다
4로 나누는 것은 4개의 particle 각각에 대한 처리라고 볼 수 있다.
initQ()
export class IsometricBendingConstraint extends BendingConstraint {
...
private initQ() {
for (let i = 0; i < this.bendingIds.length / 4; i++) {
const ids = [
this.bendingIds[4 * i],
this.bendingIds[4 * i + 1],
this.bendingIds[4 * i + 2],
this.bendingIds[4 * i + 3],
];
const particles = ids.map((id) => {
const copy = new Float32Array(4);
vecCopy(copy, 0, this.positions, id);
return copy;
});
const cotTheta = (a: DataArray, b: DataArray) => {
const cosTheta = vecDot(a, 0, b, 0);
const cross = new Float32Array(3);
vecSetCross(cross, 0, a, 0, b, 0);
const sinTheta = vecNorm(cross);
return cosTheta / sinTheta;
};
const [p0, p1, p2, p3] = particles;
...isometric bending constraint 에서는 생성자에서 Q를 초기화하는 함수를 호출한다.
한 그룹을 이루는 4개 파티클씩 ids에 할당하고
position을 particles 변수에 잠깐 넣었다가 p0 p1 p2 p3 4개 단위로 position을 할당한다
이 4개의 vertex position p0 p1 p2 p3을 가지고 edge들의 길이를 구한다
const cotTheta = (a: DataArray, b: DataArray) => {
const cosTheta = vecDot(a, 0, b, 0);
const cross = new Float32Array(3);
vecSetCross(cross, 0, a, 0, b, 0);
const sinTheta = vecNorm(cross);
return cosTheta / sinTheta;
};isometirc bending은 theta 각도를 이용한 방법이기 때문에 cotan theta를 계산하는 함수를 정의한다
export class IsometricBendingConstraint extends BendingConstraint {
...
private initQ() {
...
const e0 = new Float32Array(3),
e1 = new Float32Array(3),
e2 = new Float32Array(3),
e3 = new Float32Array(3),
e4 = new Float32Array(3),
ne1 = new Float32Array(3),
ne2 = new Float32Array(3),
e01Cross = new Float32Array(3),
e03Cross = new Float32Array(3);
vecSetDiff(e0, 0, p1, 0, p0, 0);
vecSetDiff(e1, 0, p2, 0, p1, 0);
vecSetDiff(e2, 0, p0, 0, p2, 0);
vecSetDiff(e3, 0, p3, 0, p0, 0);
vecSetDiff(e4, 0, p1, 0, p3, 0);
vecCopy(ne1, 0, e1, 0);
vecCopy(ne2, 0, e2, 0);
vecScale(ne1, 0, -1);
vecScale(ne2, 0, -1);
vecSetCross(e01Cross, 0, e0, 0, e1, 0);
vecSetCross(e03Cross, 0, e0, 0, e3, 0);
...먼저 p0~p3 들을 이용해서 edge를 구한다.
vecSetDiff()를 이용해서 구한 edge들은 아래 그림과 같은 형태이다.
vecCopy(ne1, 0, e1, 0);
vecCopy(ne2, 0, e2, 0);
vecScale(ne1, 0, -1);
vecScale(ne2, 0, -1);edge1과 edge2를 복사한 ne1 ne2가 있다.
ne는 -1로 scale을 잡기 때문에 ne라는 변수는 negative edge일 것이다.
vecSetCross(e01Cross, 0, e0, 0, e1, 0);
vecSetCross(e03Cross, 0, e0, 0, e3, 0);edge 0,1의 외적값과 edge 0,3의 외적값을 구해놓고, 나중에 norm 즉 크기를 구할 때 사용된다.
vecSetCross(e01Cross, 0, e0, 0, e1, 0);
vecSetCross(e03Cross, 0, e0, 0, e3, 0);
const cot01 = cotTheta(e0, ne1);
const cot02 = cotTheta(e0, ne2);
const cot03 = cotTheta(e0, e3);
const cot04 = cotTheta(e0, e4);edge 0과 각 다른 선분사이의 cotangent 값들을 구한다. 이때 edge 1,2 는 negative로 구한다. (정확한 이유는 모르겠다)
const K = new Float32Array([
cot01 + cot04,
cot02 + cot03,
-cot01 - cot02,
-cot03 - cot04,
]);
const Km: number[][] = multiply4dColumnVectorByTranspose(K);
const Q = mat4.fromValues(
Km[0][0],
Km[0][1],
Km[0][2],
Km[0][3],
Km[1][0],
Km[1][1],
Km[1][2],
Km[1][3],
Km[2][0],
Km[2][1],
Km[2][2],
Km[2][3],
Km[3][0],
Km[3][1],
Km[3][2],
Km[3][3]
);cotan들로 구해낸 K와 K의 전치행렬 K^(-1)을 가지고 4*4의 행렬 Q matrix를 만든다
const A0 = 0.5 * vecNorm(e01Cross);
const A1 = 0.5 * vecNorm(e03Cross);edge0 edge1의 외적의 절반 크기를 A0이라고 하고, edge0 edge3의 외적의 절반 크기를 A1라고 한다
둘을 각각 그림에서 왼쪽 오른쪽 삼각형에 해당할 것이다
mat4.multiplyScalar(Q, Q, 3.0 / (A0 + A1));
let j = 0;
const qPart = i * 16;
while (j < 16) {
this.Q[qPart + j] = Q[j];
j++;
}
matrix Q에 3/(A0+A1)을 곱함으로써 최종적으로 행렬 Q를 구한 셈이다
Q의 16개 요소들을 전역 Q array에 각각 담아서 Q 초기화를 완성한다.
이렇게 구한 Q 행렬이 담긴 배열을 이용해서 solve에서는
C, gradient C, delta lagrangian multiplier를 계산함으로써 position을 변경한다.
solve()
solve(dt: number) {
const alpha = this.compliance / dt / dt;
const memo = new Float32Array(16);
for (let i = 0; i < this.bendingIds.length / 4; i++) {
let idx = i * 4;
const ids = [
this.bendingIds[idx++],
this.bendingIds[idx++],
this.bendingIds[idx++],
this.bendingIds[idx],
];
const qIdx = i * 16;
memo[0] = vecDot(this.positions, ids[0], this.positions, ids[0]);
memo[1] = vecDot(this.positions, ids[0], this.positions, ids[1]);
memo[2] = vecDot(this.positions, ids[0], this.positions, ids[2]);
memo[3] = vecDot(this.positions, ids[0], this.positions, ids[3]);
memo[4] = memo[1];
memo[5] = vecDot(this.positions, ids[1], this.positions, ids[1]);
memo[6] = vecDot(this.positions, ids[1], this.positions, ids[2]);
memo[7] = vecDot(this.positions, ids[1], this.positions, ids[3]);
memo[8] = memo[2];
memo[9] = memo[6];
memo[10] = vecDot(this.positions, ids[2], this.positions, ids[2]);
memo[11] = vecDot(this.positions, ids[2], this.positions, ids[3]);
memo[12] = memo[3];
memo[13] = memo[7];
memo[14] = memo[11];
memo[15] = vecDot(this.positions, ids[3], this.positions, ids[3]);
...우선 vertex 각각을 ids에 4개씩 할당하는 것 까지는 다른 constraint의 solve와 동일하다
isometric bending constraint에서는 이제 16개의 matrix Q의 id들도 사용한다
export class IsometricBendingConstraint extends BendingConstraint {
...
solve(dt: number) {
...
let C = 0;
{
for (let j = 0; j < 16; j++) {
const Q = this.Q[qIdx + j];
if (Q === 0.0) continue;
C += Q * memo[j];
}
}
// If zero, let's move on.
if (C === 0.0) continue;C를 계산한다
class의 Q array에서 16개를 가져와서 임시로 Q matrix를 만들고
positoin에 대한 외적값 행렬 memo와 곱해서 누적 값으로 C를 구한다
memo는 대각을 기준으로 대칭인 행렬이다
정확한 계산법에 대해서는 Hessian 행렬이나 jacobian 행렬 방법에 대한 사전지식이 필요할 것 같아 복잡해지니 일단 넘긴다
export class IsometricBendingConstraint extends BendingConstraint {
...
solve(dt: number) {
...
// Calculate grad
{
for (let j = 0; j < 4; j++) {
vecSetZero(this.grads, j);
}
for (let j = 0; j < 16; j++) {
const Q = this.Q[qIdx + j];
if (Q === 0.0) continue;
vecAdd(this.grads, (j / 4) << 0, this.positions, ids[j % 4], Q);
}
}
let sum = 0;
for (let j = 0; j < 4; j++) {
if (this.invMass[ids[j]] === 0.0) continue;
sum += this.invMass[ids[j]] * vecDot(this.grads, j, this.grads, j);
}
...이번에는 gradient C를 구한다
export class IsometricBendingConstraint extends BendingConstraint {
...
solve(dt: number) {
...
const deltaLagrangianMultiplier = -(0.5 * C) / (sum + alpha);
for (let j = 0; j < 4; j++) {
if (this.invMass[ids[j]] === 0.0) continue;
vecAdd(
this.positions,
ids[j],
this.grads,
j,
this.invMass[ids[j]] * deltaLagrangianMultiplier
);
}
}
}
...export function vecAdd(
a: DataArray,
anr: number,
b: DataArray,
bnr: number,
scale = 1.0
) {
anr *= 3;
bnr *= 3;
a[anr++] += b[bnr++] * scale;
a[anr++] += b[bnr++] * scale;
a[anr] += b[bnr] * scale;
}마지막으로 delta lagrangian multiplier를 구해서 역질량과 gradient에 곱함으로
다음 프레임에 update될 position을 구한다
PBD 논문의 수식을 온전히 이해하고 적용하는 것은 어렵겠지만 어떠한 절차로 이루어져있는지에 대해서 알아보았다
각 constraint마다 solve에서 시뮬레이션을 위해 처리하는 부분들에 대해 파악을 하였으니
Collision.ts에서 충돌 solve 과정에 대해서 알아보도록 하자
