이제 carmen씨의 코드를 바탕으로 cloth simulation을 동작시키는 test scene을 만들어 보려고 한다.
cloth model load하기
지금은 토끼 오브젝트를 load하고 있는데, cloth model로 변경하자.
cloth model의 경우 한 쪽만 렌더링하게 obj파일이 구성되어있는데, 이는 material을 double side로 설정해주면 된다.
clothObject.traverse((ele)=>{
if(ele instanceof Mesh){
ele.material = new MeshStandardMaterial({color: 'red', side: 2})
}
})cloth simulation 과정 분석을 시작하자
carmen씨의 webgpu기반 코드를 분석하여 cloth에 물리 시뮬레이션이 적용되도록 해볼 것이다.
import { ClothPhysicsObject } from "./PhysicsObject";
import { Mesh } from "./types";
/**
* Cloth that hangs (like from a clothesline)
* Uses position based dynamics constraints
*
* It is the user's responsibility to register constraints within the app.
*/
export default class Cloth extends ClothPhysicsObject {
constructor(mesh: Mesh, thickness: number) {
super(mesh, thickness);
this.init();
}
private init() {
// Set top of cloth to have a mass of 0 to hold still
// in order to get hanging from clothesline visual
{
// Variables to store top row
let minX = Number.MAX_VALUE;
let maxX = -Number.MAX_VALUE;
let maxY = -Number.MAX_VALUE;
for (let i = 0; i < this.numParticles; i++) {
minX = Math.min(minX, this.positions[3 * i]);
maxX = Math.max(maxX, this.positions[3 * i]);
maxY = Math.max(maxY, this.positions[3 * i + 1]);
}
// Thickness of the edge to zero out(?)
const eps = 0.000001;
for (let i = 0; i < this.numParticles; i++) {
const x = this.positions[3 * i];
const y = this.positions[3 * i + 1];
if (y > maxY - eps && (x < minX + eps || x > maxX - eps))
// if (y > maxY - eps)
this.invMass[i] = 0.0;
}
}
}
}중요한 것은 이 Cloth class를 기반으로 particle을 계산할 수 있는 스크립트를 잘 읽고 가져오는 것이다.
cloth 스크립트가 이루어지기 이전에 physics를 계산하는 스크립트를 준비해야한다. 이 물리엔진이 상당히 복잡하다.
PBD 기법을 바탕으로 만들기 때문에, collision과 constraint에 대한 코드들도 있고 인접 파티클 접근을 하는데 빠른 연산을 위한 hash기법도 들어가있어서 쉽지 않다.
어디가 어떻게 상호작용을 하는지 전반적인 흐름을 파악해보자.
읽어봐야할 스크립트는 Cloth.ts PhysicsObject.ts Hash.ts Constraint.ts Collision.ts 들이다.
PhysicsObject.ts
export default abstract class PhysicsObject {
numParticles: number;
positions: Float32Array;
prevPositions: Float32Array;
vels: Float32Array;
invMass: Float32Array;
normals: Float32Array;
indices: Uint16Array;
neighbors: Float32Array;
constraints: Constraint[];
constraintFactory: ConstraintFactory;
collisions: Collision[];
constructor(mesh: Mesh) {
this.numParticles = mesh.positions.length / 3;
this.positions = new Float32Array(mesh.positions);
this.normals = new Float32Array(mesh.normals);
this.prevPositions = new Float32Array(mesh.positions);
this.vels = new Float32Array(3 * this.numParticles);
this.invMass = new Float32Array(this.numParticles);
this.indices = new Uint16Array(mesh.indices);
this.constraints = [];
this.collisions = [];
this.invMass = this.initInvMass();
this.neighbors = this.findTriNeighbors();
this.constraintFactory = new ConstraintFactory(
this.positions,
this.invMass,
this.indices,
this.neighbors
);
}
...우선 physics object class를 추상화를 진행하고, 밑에서 cloth simulation의 형태에 맞추어서 추상화한 물리엔진을 재정의한다.
생성자에서는 각 필요한 요소들이 정의되는데, 지금은 얼핏 느낌만 오고 정확한 것은 알기 어렵다.
solve()
solve(dt: number) {
for (let i = 0; i < this.numParticles; i++) {
// Floor collision ( we currently don't have a need for it)
let y = this.positions[3 * i + 1];
const height = -0.7;
if (y < height) {
vecCopy(this.positions, i, this.prevPositions, i);
this.positions[3 * i + 1] = height;
}
}
for (const constraint of this.constraints) {
constraint.solve(dt);
}
for (const collision of this.collisions) {
collision.solve(dt);
}
}solve 함수는 단순히 해결하다는 의미는 아니고, PBD 논문의 알고리즘상 시뮬레이션의 각 순간순간마다 필요한 연산들을 진행하는 최소 단위를 solve라고 한다.
주석을 간단히 참고하면 파티클이 floor밑으로 떨어지지 않도록 y positoin을 체크하고 높이를 조정하는 단계이다.
그 이후에는 constraint에서 solve, collision에서 solve를 각각 실행해주는데,
export abstract class Constraint {
protected positions: Float32Array;
protected invMass: Float32Array;
protected indices: Uint16Array;
protected neighbors: Float32Array;
protected grads: Float32Array;
protected compliance: number;
constructor(
positions: Float32Array,
invMass: Float32Array,
indices: Uint16Array,
neighbors: Float32Array,
compliance: number
) {
this.positions = positions;
this.invMass = invMass;
this.indices = indices;
this.neighbors = neighbors;
this.compliance = compliance;
this.grads = new Float32Array(12);
}
/**
* Updates positions during an animation step
*/
abstract solve(dt: number): void;
}solve는 Constriant class에서도 마찬가지로 추상클래스로 정의되어 있고
보이는 각각의 constraint마다 다르게 정의된다.
DistanceConstraint solve()
export class DistanceConstraint extends Constraint {
...
solve(dt: number) {
const alpha = this.compliance / dt / dt;
for (let i = 0; i < this.edgeLengths.length; i++) {
const id0 = this.edgeIds[2 * i];
const id1 = this.edgeIds[2 * i + 1];
const w0 = this.invMass[id0];
const w1 = this.invMass[id1];
const w = w0 + w1;
if (w == 0.0) continue;
vecSetDiff(this.grads, 0, this.positions, id0, this.positions, id1);
const len = Math.sqrt(vecLengthSquared(this.grads, 0));
if (len == 0.0) continue;
const restLen = this.edgeLengths[i];
const C = len - restLen;
const normalizingFactor = 1.0 / len;
const s = (-C / (w + alpha)) * normalizingFactor;
vecAdd(this.positions, id0, this.grads, 0, s * w0);
vecAdd(this.positions, id1, this.grads, 0, -s * w1);
}
}DistanceConstraint class의 solve를 보자.
compliance는 물질이 변하기 쉬운 정도 즉 유연성을 의미하는 상수이다.
delta time 제곱으로 compliance를 나누어 alpha라고 하고, alpha는 s라는 값을 구할 때 사용된다.
반복문은 edge의 갯수만큼 순환하는데, distance constraint라는 이름으로 생각해보면 파티클끼리 떨어져있는 거리와 관련되어있는 constraint일 것이다. 예를 들면 edge의 길이보다 파티클이 더 멀리 떨어지게 되었을 때 원래 거리로 돌아오는 방향으로 속력을 만드는 constraint이라던가.
const id0 = this.edgeIds[2 * i];
const id1 = this.edgeIds[2 * i + 1];
const w0 = this.invMass[id0];
const w1 = this.invMass[id1];
const w = w0 + w1;
if (w == 0.0) continue;보통 2x, 2x+1, 2x+2 이렇게 세 개가 한 삼각형을 이룬다. edge를 계산할 때 모든 edge 쌍들에 대해서 각각 invMass라는 것을 구하고 이를 합산해서 0이면 나머지 과정을 스킵한다.
여기서 invMass는 inverse mass로 보통 물리식에서 많이 보이는 질량의 역수를 말한다.
코드를 보면 edge의 id를 구해놓고 같은 index로 파티클의 역질량을 가져오는 걸 볼 수 있는데
그렇다면 edge와 vertex간의 인덱스가 어느정도 일치한다는 것을 말한다. 즉 1번 vertex에는 1번 edge가 연결되어있고, vertex 2에는 edge 2가 연결되는 규칙이 유지된다는 뜻이다.
/**
* Subtract two 3-element vectors
*/
export function vecSetDiff(
dst: DataArray,
dnr: number,
a: DataArray,
anr: number,
b: DataArray,
bnr: number,
scale = 1.0
) {
dnr *= 3;
anr *= 3;
bnr *= 3;
dst[dnr++] = (a[anr++] - b[bnr++]) * scale;
dst[dnr++] = (a[anr++] - b[bnr++]) * scale;
dst[dnr] = (a[anr] - b[bnr]) * scale;
}vecSetDiff 함수는 vector a와 b의 차를 구하는 식이다. 요소 3개에 대해서 계산하므로 dnr, anr, bnr의 nr은 number로 생각된다.
a와 b 벡터의 차를 구하여 dst distance에 저장한다. 따라서 vecSetDiff의 dst에 해당하는 인자인 this.grads에 들어간다.
grads라는 것은 아무래도 gradients를 말하는 것 같다.
export class DistanceConstraint extends Constraint {
...
solve(dt: number) {
...
vecSetDiff(this.grads, 0, this.positions, id0, this.positions, id1);
const len = Math.sqrt(vecLengthSquared(this.grads, 0));
if (len == 0.0) continue;
const restLen = this.edgeLengths[i];
const C = len - restLen;
const normalizingFactor = 1.0 / len;
const s = (-C / (w + alpha)) * normalizingFactor;
vecAdd(this.positions, id0, this.grads, 0, s * w0);
vecAdd(this.positions, id1, this.grads, 0, -s * w1);
}
}구한 gradient를 이용해서 거리를 구한다.
제곱을 sqrt해서 length를 구하는 과정인데
/**
* Find the length of a 3-element vector within a DataArray
*/
export function vecLengthSquared(a: DataArray, anr: number): number {
anr *= 3;
let a0 = a[anr],
a1 = a[anr + 1],
a2 = a[anr + 2];
return a0 * a0 + a1 * a1 + a2 * a2;
}vecLengthSquared라는 함수에서는 마찬가지로 anr을 사용한다. 아까부터 anr에 3을 곱하던데, 왜일까.
norm? number? 어떤 변수인지 확실하진 않다.
// Get length of the grad vector before normalization
const length = Math.sqrt(grad[0] ** 2 + grad[1] ** 2 + grad[2] ** 2);그래서 carmen씨의 이전 버전의 코드를 뒤져봤다.
위 식이 vecLengthSquared에 해당하는 코드인데, 역시 anr에 해당하는 것은 grad array의 index였나보다.
예상대로 index number인 anr와 particle 및 gradient 배열 사이에 3배수의 대칭이 되는 것인듯 하다.
1번 인덱스에 vertex 0 1 2, 2번 인덱스에 vertex 3 4 5 이런 흐름이라면 3배수 대칭이 맞다.
돌아와서, 구한 length가 0이면 또 스킵한다. 역질량이 0(질량이 무한대)이거나 길이가 0이면 solve를 적용할 수 없다.
export class DistanceConstraint extends Constraint {
...
solve(dt: number) {
...
const restLen = this.edgeLengths[i];
const C = len - restLen;
const normalizingFactor = 1.0 / len;
const s = (-C / (w + alpha)) * normalizingFactor;
vecAdd(this.positions, id0, this.grads, 0, s * w0);
vecAdd(this.positions, id1, this.grads, 0, -s * w1);
}
}constraint에서는 edge의 길이(restLen)와 현재 변화된 파티클 사이의 길이(len)의 차이를 구해서 위치를 이동시키게 된다.
PBD의 distance constriant와 같이 기존 거리보다 떨어지면 떨어진 만큼 돌아오는 힘을 준다고 생각하면 된다.
그 다음에 len의 역수를 구해서 normalize하는데
위 식에서 del C가 normalizingFactor에 해당한다.
-C에서 구했던 역질량 w와 물질상수 alpha합을 나누고, normalizingFactor를 곱하여 s(lambda)를 구한다.
export class DistanceConstraint extends Constraint {
...
solve(dt: number) {
...
vecAdd(this.positions, id0, this.grads, 0, s * w0);
vecAdd(this.positions, id1, this.grads, 0, -s * w1);
}
}
...
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;
}구한 s값은 delta lagrangian multiplier라는 이름으로 carmen씨의 과거 코드에 명명되어있다.
이 값은 gradient 및 역질량과 곱함으로써 새로운 position을 정해준다.
다음에는 BendingConstraint에 대해서 알아보도록 하자.
