Babylon.js WebGPU Cloth Simulation (2)
前回アップした、XPBD(extended position-based dynamics) による Cloth シミュレーションの解説です。
１ GPU StorageBuffer
compute shader で扱うデータは、storage buffer を介して利用する。
mesh の各vertex point のデータは、Particle 構造体とし、そのメンバーとして pos(position), vel(velocity), prePos(previous position), invMass(inverse mass) を設定し、storage buffer の変数とする。
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
prePos : vec3<f32>,
@align(16) invMass: f32,
};
@group(0) @binding(0) var<storage, read_write> particles : array<Particle>;メンバー変数の alignment は vec3
alignment は、メンバー変数の最大値となるので、16 である。構造体のメモリー配置が 16 の倍数となるように、invMass 変数の前に @align(16) を入れる。これで正しくメモリー配置される。
CPU 側でデータを Storage Buffer にセットする場合、以下のようにする。
this.numParticles は Particle 数を表す変数。
this.particleData = new Float32Array(this.numParticles * 16);
for (let i = 0; i < this.numParticles; ++i) {
let p = this.points[i];
// pos
this.particleData[16 * i + 0] = p.pos.x;
this.particleData[16 * i + 1] = p.pos.y;
this.particleData[16 * i + 2] = p.pos.z;
this.particleData[16 * i + 3] = 1.0;
// vel
this.particleData[16 * i + 4] = p.vel.x;
this.particleData[16 * i + 5] = p.vel.y;
this.particleData[16 * i + 6] = p.vel.z;
this.particleData[16 * i + 7] = 1.0;
// prePos
this.particleData[16 * i + 8] = p.prevPos.x;
this.particleData[16 * i + 9] = p.prevPos.y;
this.particleData[16 * i + 10] = p.prevPos.z;
this.particleData[16 * i + 11] = 1.0;
// invMass
this.particleData[16 * i + 12] = p.invMass;
// dummy
this.particleData[16 * i + 13] = 0.0;
this.particleData[16 * i + 14] = 0.0;
this.particleData[16 * i + 15] = 0.0;
}
this.particleBuffers = new BABYLON.StorageBuffer(engine,
this.particleData.byteLength,
BABYLON.Constants.BUFFER_CREATIONFLAG_VERTEX |
BABYLON.Constants.BUFFER_CREATIONFLAG_READWRITE);
this.particleBuffers.update(this.particleData);GPU 側での alignment に合わせて、position データ等は 16 バイト（4 バイト X 4）にする。
// pos this.particleData[16 * i + 0] = p.pos.x; this.particleData[16 * i + 1] = p.pos.y; this.particleData[16 * i + 2] = p.pos.z; this.particleData[16 * i + 3] = 1.0;
invMass データも alignment 16 にしてあるので、もとの invMass データに加えて、12 バイト分のダミーのデータを入れる。
// invMass this.particleData[16 * i + 12] = p.invMass; // dummy this.particleData[16 * i + 13] = 0.0; this.particleData[16 * i + 14] = 0.0; this.particleData[16 * i + 15] = 0.0;
２ solve Constraints
GPU 計算で拘束を解く方法は、参照ページ “16 Simulation on the GPU” の “Graph Coloring” の方法で行なっている。
Stretching、Shearing、Bending Constraints で計算する２点は、以下のようにする。GPU 計算では、１回の dispatch（実行） が終了するまで buffer の値は更新されない。そのため、mesh point が重ならないように、passNr や iw、ih パラメータを使って計算する mesh point をコントロールする。
Stretching Constraints
Stretch 1
Stretch 2
分割数 が 6 の場合、Stretch 1 における passNr 0, 1 に対応した２点の配置は次のようになる。
Shearing Constraints
Shear 1
Shear 2
Bending Constraints
Bend 1
Bend 2
Stretch 1 を計算する compute shader
solveStretchConst1
struct Params {
passNr: u32,
};
@group(0) @binding(0) var<uniform> params : Params;
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
prePos : vec3<f32>,
@align(16) invMass: f32,
};
@group(0) @binding(1) var<storage, read_write> particles : array<Particle>;
@compute @workgroup_size(1)
// dispatch iw_max
fn main(@builtin(global_invocation_id) GlobalInvocationID : vec3<u32>) {
let numSubsteps: u32 = ${numSubsteps};
let dt: f32 = ${dt};
let sdt: f32 = dt / f32(numSubsteps);
let d0: f32 = ${this.defaultSpacingY};
let compliance: f32 = ${stretchCompliance};
let iw_max: u32 = ${this.iw_max};
let ih_max: u32 = ${this.ih_max};
let passNr: u32 = params.passNr;
var index : u32 = GlobalInvocationID.x;
for (var ih: u32 = passNr; ih <= (ih_max - 1u) - 2u + passNr; ih = ih + 2u) {
var id0: u32 = index + ih * iw_max;
var id1: u32 = index + (ih + 1u) * iw_max;
var w0: f32 = particles[id0].invMass;
var w1: f32 = particles[id1].invMass;
var w: f32 = w0 + w1;
if (w == 0.0) {return;}
var p0 : vec3<f32> = particles[id0].pos;
var p1 : vec3<f32> = particles[id1].pos;
var grad: vec3<f32> = p1 - p0;
var d: f32 = length(grad);
let alpha: f32 = compliance / sdt / sdt;
var lambda: f32 = (d - d0) / (w + alpha);
grad = normalize(grad);
var delta_p0: vec3<f32> = grad * (w0 * lambda);
particles[id0].pos = p0 + delta_p0;
var delta_p1: vec3<f32> = grad * (-w1 * lambda);
particles[id1].pos = p1 + delta_p1;
}
}拘束を解く式は、
J. Bender, M. Müller, M. Macklin
Position-Based Simulation Methods in Computer Graphics 2017
EUROGRAPHICS 2017 Tutorial (2017)
https://matthias-research.github.io/pages/publications/PBDTutorial2017-CourseNotes.pdf
Course notes の 4.2.6 XPBD と 5.1 Stretching の distance constraints の式を使う。
2点 ,
の位置補正の式は次のようになる。
: rest length ,
: compliance
,
: inverse mass of particles 1, 2
Shearing Bending にも同じ distance constraints の式を使う。
３ Uniform Buffer を用いた繰返し処理
拘束を解く際に用いる passNr や iw、ih の値は uniform buffer を通して shader に渡し、それらの値を変えて繰返し処理を行う。
uniform のある値に対し shader を計算し、次に新しい値に対して同じ shader を使って計算する場合は、別々の ComputeShader object として定義する。そして、それぞれの ComputeShader を使って、uniform の値を変えて計算する。
ComputeShader の定義
// cs solveStretchConst1
this.csSolveStretchConst1 = [];
for (let passNr = 0; passNr < 2; passNr++) {
const csSolve1 = new BABYLON.ComputeShader("csSolve1", engine,
{ computeSource: computeShader.solveStretchConst1 },
{ bindingsMapping: {
"params": { group: 0, binding: 0 },
"particle": { group: 0, binding: 1 },
}
});
csSolve1.setUniformBuffer("params", this.uBuffer);
csSolve1.setStorageBuffer("particle", this.particleBuffers);
this.csSolveStretchConst1.push(csSolve1);
}ComputeShader の実行
for (let passNr = 0; passNr < 2; passNr++) {
this.updateParams(passNr);
this.csSolveStretchConst1[passNr].dispatch(this.iw_max);
}UniformBuffer の設定とupdate メソッド
// Uniform buffer
this.uBuffer = new BABYLON.UniformBuffer(engine, undefined, undefined,
"uBuffer");
this.uBuffer.addUniform("passNr", 1);
updateParams(passNr) {
this.uBuffer.updateInt("passNr", passNr);
this.uBuffer.update();
}４ ワークグループ
compute shader の integrate と updateVel では GPU のワークグループサイズを 64 に設定している。
@workgroup_size(64)
integrate では、重力による速度の変化とそれによる位置の変化を計算している。updateVel では、拘束を解いた後の速度を計算している。この部分は、各 mesh point で同じ処理を行なっているので、並列計算が可能である。ワークグループサイズを 64 に設定し、同時処理する mesh point 毎の処理（ワークアイテム）の数を 64 にしている。
（ワークグループサイズの上限はより大きな値であるが、ここでは 64 にしている。）
shader を実行する際は、dispatch のワークグループ数を numParticles / 64 にする。
// Integrate this.csIntegrate.dispatch(Math.ceil(this.numParticles / 64));
Babylon.js WebGPU Cloth Simulation
Babylon.js WebGPU を用いて、XPBD(extended position-based dynamics) による Cloth シミュレーションを行なってみました。
XPBD の計算は、compute shader で GPU 計算を行なっています。
XPBD による Cloth シミュレーションは、下記を参照しました。
https://matthias-research.github.io/pages/
Ten Minute Physics
14 The secret of cloth simulation
16 Simulation on the GPU
14 では、XPBD は CPU 計算を行い、メッシュオブジェクトの検出に Three.js のレイキャスト を使っています。
16 では、Python 言語でプログラムを記述しており、Warp ライブラリーを用いて XPBD のGPU 計算を行なっています。メッシュオブジェクトの検出には、OpenGL と GLU、GLUT ライブラリーを用いています。
ここでは、この２つを組み合わせる形で、Babylon.js の WebGPU を用いて、compute shader で XPBD の GPU 計算を行い、レイキャストを用いてメッシュオブジェクトの検出を行なっている。
実行結果
プログラム
ClothSim-xpbd-WebGPU.html
<!DOCTYPE html>
<head>
<title>Babylon.js Cloth Sim</title>
<script src="https://preview.babylonjs.com/babylon.js"></script>
<script src="customMesh3D.js"></script>
</head>
<body>
<canvas id="renderCanvas" width="450" height="450"></canvas>
<script>
async function init() {
const canvas = document.getElementById("renderCanvas");
const engine = new BABYLON.WebGPUEngine(canvas);
await engine.initAsync();
const scene = new BABYLON.Scene(engine);
scene.clearColor = new BABYLON.Color3(0.8, 0.8, 0.8);
const camera = new BABYLON.ArcRotateCamera("camera",
-Math.PI / 2, Math.PI / 2, 8, new BABYLON.Vector3(0, 0.0, 0), scene);
camera.attachControl(canvas, true);
const light = new BABYLON.HemisphericLight("light1",
new BABYLON.Vector3(0, 1, 0), scene);
light.intensity = 0.8;
// Constants
let gravity = [0.0, -10.0, 0.0];
let dt = 1.0 / 60.0;
let numSubsteps = 5;
let cwidth = 2;
let cheight = 3;
let cdivisions = 50;
cdivisions = (cdivisions%2 == 0) ? cdivisions : cdivisions - 1;
let stretchCompliance = 0.0000001;
let bendingCompliance = 0.5;
class Cloth {
constructor(mesh, scene) {
let meshMat = new BABYLON.StandardMaterial("meshMat", scene);
meshMat.diffuseColor = BABYLON.Color3.Red();
meshMat.backFaceCulling = false;
this.customMesh = new BABYLON.Mesh("custom", scene);
this.customMesh.material = meshMat;
this.customMesh.userData = this; // for raycasting
let vertexData = new BABYLON.VertexData();
vertexData.positions = mesh.vertices;
vertexData.indices = mesh.indices;
vertexData.applyToMesh(this.customMesh, true);
this.numParticles = mesh.vertices.length / 3;
console.log("numParticles", this.numParticles);
this.iw_max = cdivisions + 1;
this.ih_max = cdivisions + 1;
let defaultMass = 1.0;
this.defaultSpacingX = cwidth / cdivisions;
this.defaultSpacingY = cheight / cdivisions;
// cpu data
this.points = [];
for (let i = 0; i < this.numParticles; ++i) {
let iw = i % this.iw_max;
let ih = Math.floor(i / this.iw_max);
let invMass = ((ih == 0) && (iw == 0) ||
(ih == 0) && (iw == this.iw_max-1)) ? 0 : 1 / defaultMass;
this.points.push({
pos: new BABYLON.Vector3(
mesh.vertices[3 * i + 0],
mesh.vertices[3 * i + 1],
mesh.vertices[3 * i + 2]),
vel: new BABYLON.Vector3(0, 0, 0),
prevPos: new BABYLON.Vector3(
mesh.vertices[3 * i + 0],
mesh.vertices[3 * i + 1],
mesh.vertices[3 * i + 2]),
invMass: invMass,
});
}
// buffer data
this.particleData = new Float32Array(this.numParticles * 16);
for (let i = 0; i < this.numParticles; ++i) {
let p = this.points[i];
// pos
this.particleData[16 * i + 0] = p.pos.x;
this.particleData[16 * i + 1] = p.pos.y;
this.particleData[16 * i + 2] = p.pos.z;
this.particleData[16 * i + 3] = 1.0;
// vel
this.particleData[16 * i + 4] = p.vel.x;
this.particleData[16 * i + 5] = p.vel.y;
this.particleData[16 * i + 6] = p.vel.z;
this.particleData[16 * i + 7] = 1.0;
// prePos
this.particleData[16 * i + 8] = p.prevPos.x;
this.particleData[16 * i + 9] = p.prevPos.y;
this.particleData[16 * i + 10] = p.prevPos.z;
this.particleData[16 * i + 11] = 1.0;
// invMass
this.particleData[16 * i + 12] = p.invMass;
// dummy
this.particleData[16 * i + 13] = 0.0;
this.particleData[16 * i + 14] = 0.0;
this.particleData[16 * i + 15] = 0.0;
}
console.log("particleData", Array.from(this.particleData));
this.particleBuffers = new BABYLON.StorageBuffer(engine,
this.particleData.byteLength,
BABYLON.Constants.BUFFER_CREATIONFLAG_VERTEX | BABYLON.Constants.BUFFER_CREATIONFLAG_READWRITE);
this.particleBuffers.update(this.particleData);
this.grabId = -1;
this.grabInvMass = 0.0;
this.grabIdData = new Int32Array([this.grabId]);
this.grabIdBuffer = new BABYLON.StorageBuffer(engine,
this.grabIdData.byteLength,
BABYLON.Constants.BUFFER_CREATIONFLAG_READWRITE);
this.grabIdBuffer.update(this.grabIdData);
this.grabInvMassData = new Float32Array([this.grabInvMass]);
this.grabInvMassBuffer = new BABYLON.StorageBuffer(engine,
this.grabInvMassData.byteLength,
BABYLON.Constants.BUFFER_CREATIONFLAG_READWRITE);
this.grabInvMassBuffer.update(this.grabInvMassData);
// Uniform buffer
this.uBuffer = new BABYLON.UniformBuffer(engine, undefined, undefined, "uBuffer");
this.uBuffer.addUniform("passNr", 1);
// uniform startGrab
this.grab1UBuffer = new BABYLON.UniformBuffer(engine, undefined, undefined, "grab1UBuffer");
this.grab1UBuffer.addUniform("pos", 3);
// unihorm moveGrabbed
this.grab12UBuffer = new BABYLON.UniformBuffer(engine, undefined, undefined, "grab12UBuffer");
this.grab12UBuffer.addUniform("pos", 3);
// uniform endGrab
this.grab2UBuffer = new BABYLON.UniformBuffer(engine, undefined, undefined, "grab2UBuffer");
this.grab2UBuffer.addUniform("vel", 3);
this.setComputeShader();
}
setComputeShader() {
const computeShader = {
integrate:`
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
prePos : vec3<f32>,
@align(16) invMass: f32,
};
@group(0) @binding(0) var<storage, read_write> particles : array<Particle>;
@compute @workgroup_size(64)
fn main(@builtin(global_invocation_id) GlobalInvocationID : vec3<u32>) {
let numSubsteps: u32 = ${numSubsteps};
let dt: f32 = ${dt};
let sdt: f32 = dt / f32(numSubsteps);
let gravity: f32 = ${gravity[1]};
var index : u32 = GlobalInvocationID.x;
let vInvMass: f32 = particles[index].invMass;
if (vInvMass == 0.0) { return; }
var vPos : vec3<f32> = particles[index].pos;
var vVel : vec3<f32> = particles[index].vel;
vVel.y += gravity * sdt;
particles[index].prePos = vPos;
vPos = vPos + vVel * sdt;
particles[index].pos = vPos;
particles[index].vel = vVel;
}`,
solveStretchConst1:`
struct Params {
passNr: u32,
};
@group(0) @binding(0) var<uniform> params : Params;
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
prePos : vec3<f32>,
@align(16) invMass: f32,
};
@group(0) @binding(1) var<storage, read_write> particles : array<Particle>;
@compute @workgroup_size(1)
// dispatch iw_max
fn main(@builtin(global_invocation_id) GlobalInvocationID : vec3<u32>) {
let numSubsteps: u32 = ${numSubsteps};
let dt: f32 = ${dt};
let sdt: f32 = dt / f32(numSubsteps);
let d0: f32 = ${this.defaultSpacingY};
let compliance: f32 = ${stretchCompliance};
let iw_max: u32 = ${this.iw_max};
let ih_max: u32 = ${this.ih_max};
let passNr: u32 = params.passNr;
var index : u32 = GlobalInvocationID.x;
for (var ih: u32 = passNr; ih <= (ih_max - 1u) - 2u + passNr; ih = ih + 2u) {
var id0: u32 = index + ih * iw_max;
var id1: u32 = index + (ih + 1u) * iw_max;
var w0: f32 = particles[id0].invMass;
var w1: f32 = particles[id1].invMass;
var w: f32 = w0 + w1;
if (w == 0.0) {return;}
var p0 : vec3<f32> = particles[id0].pos;
var p1 : vec3<f32> = particles[id1].pos;
var grad: vec3<f32> = p1 - p0;
var d: f32 = length(grad);
let alpha: f32 = compliance / sdt / sdt;
var lambda: f32 = (d - d0) / (w + alpha);
grad = normalize(grad);
var delta_p0: vec3<f32> = grad * (w0 * lambda);
particles[id0].pos = p0 + delta_p0;
var delta_p1: vec3<f32> = grad * (-w1 * lambda);
particles[id1].pos = p1 + delta_p1;
}
}`,
solveStretchConst2:`
struct Params {
passNr: u32,
};
@group(0) @binding(0) var<uniform> params : Params;
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
prePos : vec3<f32>,
@align(16) invMass: f32,
};
@group(0) @binding(1) var<storage, read_write> particles : array<Particle>;
@compute @workgroup_size(1)
// dispatch ih_max
fn main(@builtin(global_invocation_id) GlobalInvocationID : vec3<u32>) {
let numSubsteps: u32 = ${numSubsteps};
let dt: f32 = ${dt};
let sdt: f32 = dt / f32(numSubsteps);
let d0: f32 = ${this.defaultSpacingX};
let compliance: f32 = ${stretchCompliance};
let iw_max: u32 = ${this.iw_max};
let passNr: u32 = params.passNr;
var index : u32 = GlobalInvocationID.x;
for (var iw: u32 = passNr; iw <= (iw_max - 1u) - 2u + passNr; iw = iw + 2u) {
var id0: u32 = iw + index * iw_max;
var id1: u32 = (iw + 1u) + index * iw_max;
var w0: f32 = particles[id0].invMass;
var w1: f32 = particles[id1].invMass;
var w: f32 = w0 + w1;
if (w == 0.0) {return;}
var p0 : vec3<f32> = particles[id0].pos;
var p1 : vec3<f32> = particles[id1].pos;
var grad: vec3<f32> = p1 - p0;
var d: f32 = length(grad);
let alpha: f32 = compliance / sdt / sdt;
var lambda: f32 = (d - d0) / (w + alpha);
grad = normalize(grad);
var delta_p0: vec3<f32> = grad * (w0 * lambda);
particles[id0].pos = p0 + delta_p0;
var delta_p1: vec3<f32> = grad * (-w1 * lambda);
particles[id1].pos = p1 + delta_p1;
}
}`,
solveShearConst1:`
struct Params {
passNr: u32,
};
@group(0) @binding(0) var<uniform> params : Params;
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
prePos : vec3<f32>,
@align(16) invMass: f32,
};
@group(0) @binding(1) var<storage, read_write> particles : array<Particle>;
@compute @workgroup_size(1)
// dispatch iw_max-1
fn main(@builtin(global_invocation_id) GlobalInvocationID : vec3<u32>) {
let numSubsteps: u32 = ${numSubsteps};
let dt: f32 = ${dt};
let sdt: f32 = dt / f32(numSubsteps);
let d0X: vec2<f32> = vec2(${this.defaultSpacingX}, 0);
let d0Y: vec2<f32> = vec2(0, ${this.defaultSpacingY});
let d0: f32 = length(d0X - d0Y);
let compliance: f32 = ${stretchCompliance};
let iw_max: u32 = ${this.iw_max};
let ih_max: u32 = ${this.ih_max};
let passNr: u32 = params.passNr;
var index : u32 = GlobalInvocationID.x;
for (var ih: u32 = passNr; ih <= (ih_max - 1u) - 2u + passNr; ih = ih + 2u) {
var id0: u32 = index + ih * iw_max;
var id1: u32 = (index + 1u) + (ih + 1u) * iw_max;
var w0: f32 = particles[id0].invMass;
var w1: f32 = particles[id1].invMass;
var w: f32 = w0 + w1;
if (w == 0.0) {return;}
var p0 : vec3<f32> = particles[id0].pos;
var p1 : vec3<f32> = particles[id1].pos;
var grad: vec3<f32> = p1 - p0;
var d: f32 = length(grad);
let alpha: f32 = compliance / sdt / sdt;
var lambda: f32 = (d - d0) / (w + alpha);
grad = normalize(grad);
var delta_p0: vec3<f32> = grad * (w0 * lambda);
particles[id0].pos = p0 + delta_p0;
var delta_p1: vec3<f32> = grad * (-w1 * lambda);
particles[id1].pos = p1 + delta_p1;
}
}`,
solveShearConst2:`
struct Params {
passNr: u32,
};
@group(0) @binding(0) var<uniform> params : Params;
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
prePos : vec3<f32>,
@align(16) invMass: f32,
};
@group(0) @binding(1) var<storage, read_write> particles : array<Particle>;
@compute @workgroup_size(1)
// dispatch ih_max-1
fn main(@builtin(global_invocation_id) GlobalInvocationID : vec3<u32>) {
let numSubsteps: u32 = ${numSubsteps};
let dt: f32 = ${dt};
let sdt: f32 = dt / f32(numSubsteps);
let d0X: vec2<f32> = vec2(${this.defaultSpacingX}, 0);
let d0Y: vec2<f32> = vec2(0, ${this.defaultSpacingY});
let d0: f32 = length(d0X + d0Y);
let compliance: f32 = ${stretchCompliance};
let iw_max: u32 = ${this.iw_max};
let passNr: u32 = params.passNr;
var index : u32 = GlobalInvocationID.x;
for (var iw: u32 = passNr; iw <= (iw_max - 1u) - 2u + passNr; iw = iw + 2u) {
var id0: u32 = (iw + 1u) + index * iw_max;
var id1: u32 = iw + (index + 1u) * iw_max;
var w0: f32 = particles[id0].invMass;
var w1: f32 = particles[id1].invMass;
var w: f32 = w0 + w1;
if (w == 0.0) {return;}
var p0 : vec3<f32> = particles[id0].pos;
var p1 : vec3<f32> = particles[id1].pos;
var grad: vec3<f32> = p1 - p0;
var d: f32 = length(grad);
let alpha: f32 = compliance / sdt / sdt;
var lambda: f32 = (d - d0) / (w + alpha);
grad = normalize(grad);
var delta_p0: vec3<f32> = grad * (w0 * lambda);
particles[id0].pos = p0 + delta_p0;
var delta_p1: vec3<f32> = grad * (-w1 * lambda);
particles[id1].pos = p1 + delta_p1;
}
}`,
solveBendConst1:`
struct Params {
passNr: u32,
};
@group(0) @binding(0) var<uniform> params : Params;
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
prePos : vec3<f32>,
@align(16) invMass: f32,
};
@group(0) @binding(1) var<storage, read_write> particles : array<Particle>;
@compute @workgroup_size(1)
// dispatch iw_max
fn main(@builtin(global_invocation_id) GlobalInvocationID : vec3<u32>) {
let numSubsteps: u32 = ${numSubsteps};
let dt: f32 = ${dt};
let sdt: f32 = dt / f32(numSubsteps);
let spacingY: f32 = ${this.defaultSpacingY};
let d0: f32 = spacingY * 2;
let compliance: f32 = ${bendingCompliance};
let iw_max: u32 = ${this.iw_max};
let ih: u32 = params.passNr;
var index : u32 = GlobalInvocationID.x;
var id0: u32 = index + ih * iw_max;
var id1: u32 = index + (ih + 2u) * iw_max;
var w0: f32 = particles[id0].invMass;
var w1: f32 = particles[id1].invMass;
var w: f32 = w0 + w1;
if (w == 0.0) {return;}
var p0 : vec3<f32> = particles[id0].pos;
var p1 : vec3<f32> = particles[id1].pos;
var grad: vec3<f32> = p1 - p0;
var d: f32 = length(grad);
let alpha: f32 = compliance / sdt / sdt;
var lambda: f32 = (d - d0) / (w + alpha);
grad = normalize(grad);
var delta_p0: vec3<f32> = grad * (w0 * lambda);
particles[id0].pos = p0 + delta_p0;
var delta_p1: vec3<f32> = grad * (-w1 * lambda);
particles[id1].pos = p1 + delta_p1;
}`,
solveBendConst2:`
struct Params {
passNr: u32,
};
@group(0) @binding(0) var<uniform> params : Params;
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
prePos : vec3<f32>,
@align(16) invMass: f32,
};
@group(0) @binding(1) var<storage, read_write> particles : array<Particle>;
@compute @workgroup_size(1)
// dispatch ih_max
fn main(@builtin(global_invocation_id) GlobalInvocationID : vec3<u32>) {
let numSubsteps: u32 = ${numSubsteps};
let dt: f32 = ${dt};
let sdt: f32 = dt / f32(numSubsteps);
let spacingX: f32 = ${this.defaultSpacingX};
let d0: f32 = spacingX * 2;
let compliance: f32 = ${bendingCompliance};
let iw_max: u32 = ${this.iw_max};
let iw: u32 = params.passNr;
var index : u32 = GlobalInvocationID.x;
var id0: u32 = iw + index * iw_max;
var id1: u32 = (iw + 2u) + index * iw_max;
var w0: f32 = particles[id0].invMass;
var w1: f32 = particles[id1].invMass;
var w: f32 = w0 + w1;
if (w == 0.0) {return;}
var p0 : vec3<f32> = particles[id0].pos;
var p1 : vec3<f32> = particles[id1].pos;
var grad: vec3<f32> = p1 - p0;
var d: f32 = length(grad);
let alpha: f32 = compliance / sdt / sdt;
var lambda: f32 = (d - d0) / (w + alpha);
grad = normalize(grad);
var delta_p0: vec3<f32> = grad * (w0 * lambda);
particles[id0].pos = p0 + delta_p0;
var delta_p1: vec3<f32> = grad * (-w1 * lambda);
particles[id1].pos = p1 + delta_p1;
}`,
updateVel:`
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
prePos : vec3<f32>,
@align(16) invMass: f32,
};
@group(0) @binding(0) var<storage, read_write> particles : array<Particle>;
@compute @workgroup_size(64)
fn main(@builtin(global_invocation_id) GlobalInvocationID : vec3<u32>) {
let numSubsteps: u32 = ${numSubsteps};
let dt: f32 = ${dt};
let sdt: f32 = dt / f32(numSubsteps);
var index : u32 = GlobalInvocationID.x;
let vInvMass: f32 = particles[index].invMass;
if (vInvMass == 0.0) { return; }
var vPos : vec3<f32> = particles[index].pos;
var vPrePos : vec3<f32> = particles[index].prePos;
var vVel : vec3<f32> = (vPos - vPrePos) / sdt;
particles[index].vel = vVel;
}`,
startGrab:`
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
prePos : vec3<f32>,
@align(16) invMass: f32,
};
@group(0) @binding(0) var<storage, read_write> particles : array<Particle>;
@group(0) @binding(1) var<storage, read_write> grabId : i32;
@group(0) @binding(2) var<storage, read_write> grabInvMass : f32;
struct Params {
pos: vec3<f32>,
};
@group(0) @binding(3) var<uniform> params : Params;
@compute @workgroup_size(1)
// dispatch: 1
fn main(@builtin(global_invocation_id) GlobalInvocationID : vec3<u32>) {
var index: u32 = GlobalInvocationID.x;
let pos: vec3<f32> = params.pos;
let gId: i32 = grabId;
var gInvMass: f32 = 0.0;
if (gId >= 0) {
gInvMass = particles[gId].invMass;
particles[gId].invMass = 0.0;
particles[gId].pos = pos;
}
grabInvMass = gInvMass;
}`,
moveGrabbed:`
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
prePos : vec3<f32>,
@align(16) invMass: f32,
};
@group(0) @binding(0) var<storage, read_write> particles : array<Particle>;
@group(0) @binding(1) var<storage, read_write> grabId : i32;
struct Params {
pos: vec3<f32>,
};
@group(0) @binding(2) var<uniform> params : Params;
@compute @workgroup_size(1)
// dispatch: 1
fn main(@builtin(global_invocation_id) GlobalInvocationID : vec3<u32>) {
var index: u32 = GlobalInvocationID.x;
let pos: vec3<f32> = params.pos;
let gId: i32 = grabId;
if (gId >= 0) {
particles[gId].pos = pos;
}
}`,
endGrab:`
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
prePos : vec3<f32>,
@align(16) invMass: f32,
};
@group(0) @binding(0) var<storage, read_write> particles : array<Particle>;
@group(0) @binding(1) var<storage, read_write> grabId : i32;
@group(0) @binding(2) var<storage, read_write> grabInvMass : f32;
struct Params {
vel: vec3<f32>,
};
@group(0) @binding(3) var<uniform> params : Params;
@compute @workgroup_size(1)
// dispatch: 1
fn main(@builtin(global_invocation_id) GlobalInvocationID : vec3<u32>) {
var index: u32 = GlobalInvocationID.x;
let vel: vec3<f32> = params.vel;
let gId: i32 = grabId;
if (gId >= 0) {
particles[gId].invMass = grabInvMass;
particles[gId].vel = vel;
}
grabId = -1;
}`
};
// cs integrate
this.csIntegrate = new BABYLON.ComputeShader("csIntegrate", engine,
{ computeSource: computeShader.integrate },
{ bindingsMapping: {
"particle": { group: 0, binding: 0 },
}
});
this.csIntegrate.setStorageBuffer("particle", this.particleBuffers);
// cs solveStretchConst1
this.csSolveStretchConst1 = [];
for (let passNr = 0; passNr < 2; passNr++) {
const csSolve1 = new BABYLON.ComputeShader("csSolve1", engine,
{ computeSource: computeShader.solveStretchConst1 },
{ bindingsMapping: {
"params": { group: 0, binding: 0 },
"particle": { group: 0, binding: 1 },
}
});
csSolve1.setUniformBuffer("params", this.uBuffer);
csSolve1.setStorageBuffer("particle", this.particleBuffers);
this.csSolveStretchConst1.push(csSolve1);
}
// cs solveStretchConst2
this.csSolveStretchConst2 = [];
for (let passNr = 0; passNr < 2; passNr++) {
const csSolve2 = new BABYLON.ComputeShader("csSolve2", engine,
{ computeSource: computeShader.solveStretchConst2 },
{ bindingsMapping: {
"params": { group: 0, binding: 0 },
"particle": { group: 0, binding: 1 },
}
});
csSolve2.setUniformBuffer("params", this.uBuffer);
csSolve2.setStorageBuffer("particle", this.particleBuffers);
this.csSolveStretchConst2.push(csSolve2);
}
// cs solveShearConst1
this.csSolveShearConst1 = [];
for (let passNr = 0; passNr < 2; passNr++) {
const csShear1 = new BABYLON.ComputeShader("csShear1", engine,
{ computeSource: computeShader.solveShearConst1 },
{ bindingsMapping: {
"params": { group: 0, binding: 0 },
"particle": { group: 0, binding: 1 },
}
});
csShear1.setUniformBuffer("params", this.uBuffer);
csShear1.setStorageBuffer("particle", this.particleBuffers);
this.csSolveShearConst1.push(csShear1);
}
this.csSolveShearConst2 = [];
for (let passNr = 0; passNr < 2; passNr++) {
const csShear2 = new BABYLON.ComputeShader("csShear2", engine,
{ computeSource: computeShader.solveShearConst2 },
{ bindingsMapping: {
"params": { group: 0, binding: 0 },
"particle": { group: 0, binding: 1 },
}
});
csShear2.setUniformBuffer("params", this.uBuffer);
csShear2.setStorageBuffer("particle", this.particleBuffers);
this.csSolveShearConst2.push(csShear2);
}
// cs solveBendConst1
this.csSolveBendConst1 = [];
// passNr -> ih
for (let ih = 0; ih < (this.ih_max-1) - 1; ih++) {
const csBend1 = new BABYLON.ComputeShader("csBend1", engine,
{ computeSource: computeShader.solveBendConst1 },
{ bindingsMapping: {
"params": { group: 0, binding: 0 },
"particle": { group: 0, binding: 1 },
}
});
csBend1.setUniformBuffer("params", this.uBuffer);
csBend1.setStorageBuffer("particle", this.particleBuffers);
this.csSolveBendConst1.push(csBend1);
}
// cs solveBendConst2
this.csSolveBendConst2 = [];
// passNr -> iw
for (let iw = 0; iw < (this.iw_max-1) - 1; iw++) {
const csBend2 = new BABYLON.ComputeShader("csBend2", engine,
{ computeSource: computeShader.solveBendConst2 },
{ bindingsMapping: {
"params": { group: 0, binding: 0 },
"particle": { group: 0, binding: 1 },
}
});
csBend2.setUniformBuffer("params", this.uBuffer);
csBend2.setStorageBuffer("particle", this.particleBuffers);
this.csSolveBendConst2.push(csBend2);
}
// cs updateVel
this.csUpdateVel = new BABYLON.ComputeShader("csUpdateVel", engine,
{ computeSource: computeShader.updateVel },
{ bindingsMapping: {
"particle": { group: 0, binding: 0 },
}
});
this.csUpdateVel.setStorageBuffer("particle", this.particleBuffers);
// cs startGrab
this.csStartGrab = new BABYLON.ComputeShader("csStartGrab", engine,
{ computeSource: computeShader.startGrab },
{ bindingsMapping: {
"particle": { group: 0, binding: 0 },
"grabId": { group: 0, binding: 1 },
"grabInvMass": { group: 0, binding: 2 },
"params": { group: 0, binding: 3 },
}
});
this.csStartGrab.setStorageBuffer("particle", this.particleBuffers);
this.csStartGrab.setStorageBuffer("grabId", this.grabIdBuffer);
this.csStartGrab.setStorageBuffer("grabInvMass", this.grabInvMassBuffer);
this.csStartGrab.setUniformBuffer("params", this.grab1UBuffer);
// cs moveGrabbed
this.csMoveGrabbed = new BABYLON.ComputeShader("csMoveGrabbed", engine,
{ computeSource: computeShader.moveGrabbed },
{ bindingsMapping: {
"particle": { group: 0, binding: 0 },
"grabId": { group: 0, binding: 1 },
"params": { group: 0, binding: 2 },
}
});
this.csMoveGrabbed.setStorageBuffer("particle", this.particleBuffers);
this.csMoveGrabbed.setStorageBuffer("grabId", this.grabIdBuffer);
this.csMoveGrabbed.setUniformBuffer("params", this.grab12UBuffer);
// cs endGrab
this.csEndGrab = new BABYLON.ComputeShader("csEndGrab", engine,
{ computeSource: computeShader.endGrab },
{ bindingsMapping: {
"particle": { group: 0, binding: 0 },
"grabId": { group: 0, binding: 1 },
"grabInvMass": { group: 0, binding: 2 },
"params": { group: 0, binding: 3 },
}
});
this.csEndGrab.setStorageBuffer("particle", this.particleBuffers);
this.csEndGrab.setStorageBuffer("grabId", this.grabIdBuffer);
this.csEndGrab.setStorageBuffer("grabInvMass", this.grabInvMassBuffer);
this.csEndGrab.setUniformBuffer("params", this.grab2UBuffer);
}
updateParams(passNr) {
this.uBuffer.updateInt("passNr", passNr);
this.uBuffer.update();
}
simulate() {
// Integrate
this.csIntegrate.dispatch(Math.ceil(this.numParticles / 64));
// solve stretch constraints
for (let passNr = 0; passNr < 2; passNr++) {
this.updateParams(passNr);
this.csSolveStretchConst1[passNr].dispatch(this.iw_max);
}
for (let passNr = 0; passNr < 2; passNr++) {
this.updateParams(passNr);
this.csSolveStretchConst2[passNr].dispatch(this.ih_max);
}
// solve shear constraints
for (let passNr = 0; passNr < 2; passNr++) {
this.updateParams(passNr);
this.csSolveShearConst1[passNr].dispatch(this.iw_max-1);
}
for (let passNr = 0; passNr < 2; passNr++) {
this.updateParams(passNr);
this.csSolveShearConst2[passNr].dispatch(this.ih_max-1);
}
// solve bending constraints
for (let ih = 0; ih < (this.ih_max-1) - 1; ih++) {
this.updateParams(ih);
this.csSolveBendConst1[ih].dispatch(this.iw_max);
}
for (let iw = 0; iw < (this.iw_max-1 - 1); iw++) {
this.updateParams(iw);
this.csSolveBendConst2[iw].dispatch(this.ih_max);
}
// update velocities
this.csUpdateVel.dispatch(Math.ceil(this.numParticles / 64));
}
updateMeshes() {
const meshVertexBuffer = this.customMesh.getVertexBuffer("position");
this.particleBuffers.read().then((res) => {
const resFloats = new Float32Array(res.buffer);
let pos = [];
for (let i = 0; i < this.numParticles; ++i) {
pos.push(
resFloats[16 * i],
resFloats[16 * i + 1],
resFloats[16 * i + 2]);
}
meshVertexBuffer.update(pos);
});
}
startGrab(pos) {
this.grab1UBuffer.updateVector3("pos", pos);
this.grab1UBuffer.update();
this.grabIdBuffer.update(new Int32Array([this.grabId]));
this.csStartGrab.dispatch(1);
this.grabIdBuffer.read().then((res) => {
const resFloats = new Int32Array(res.buffer);
console.log("res gId", resFloats[0]);
});
this.grabInvMassBuffer.read().then((res) => {
const resFloats = new Float32Array(res.buffer);
console.log("res gInvMass", resFloats[0]);
});
}
moveGrabbed(pos, vel) {
this.grab12UBuffer.updateVector3("pos", pos);
this.grab12UBuffer.update();
this.csMoveGrabbed.dispatch(1);
}
endGrab(pos, vel) {
this.grab2UBuffer.updateVector3("vel", vel);
this.grab2UBuffer.update();
this.csEndGrab.dispatch(1);
}
}
//-- Grabber --------
class Grabber {
constructor() {
this.physicsObject = null;
this.distance = 0.0;
this.prevPos = new BABYLON.Vector3();
this.vel = new BABYLON.Vector3();
this.time = 0.0;
}
increaseTime(dt) {
this.time += dt;
}
start(pickInfo) {
this.physicsObject = null;
let obj = pickInfo.pickedMesh.userData;
let indices = pickInfo.pickedMesh.getIndices();
console.log("faceId", pickInfo.faceId,
"index", indices[pickInfo.faceId * 3], indices[pickInfo.faceId * 3 + 1], indices[pickInfo.faceId * 3 + 2]);
this.physicsObject = obj;
this.distance = pickInfo.distance;
this.ray = pickInfo.ray;
let pos = pickInfo.pickedPoint;
let minD2 = Number.MAX_VALUE;
let gId = -1;
for (let j = 0; j < 3; j++) {
let idx = indices[pickInfo.faceId * 3 + j];
let obj_pos = this.physicsObject.points[idx].pos;
let d2 = BABYLON.Vector3.DistanceSquared(pos, obj_pos);
if (d2 < minD2) {
minD2 = d2;
gId = idx;
}
}
console.log("gId", gId);
obj.grabId = gId;
this.physicsObject.startGrab(pos);
this.prevPos.copyFrom(pos);
this.vel.set(0.0, 0.0, 0.0);
this.time = 0.0;
}
move(pickInfo) {
if (this.physicsObject) {
let pos = pickInfo.ray.origin.clone();
pos = pos.add(pickInfo.ray.direction.scale(this.distance));
this.vel.copyFrom(pos);
this.vel = this.vel.subtract(this.prevPos);
if (this.time > 0.0) {
this.vel.scale(1 / this.time);
}
else {
this.vel.set(0.0, 0.0, 0.0);
}
this.prevPos.copyFrom(pos);
this.time = 0.0;
this.physicsObject.moveGrabbed(pos, this.vel);
}
}
end(pickInfo) {
if (this.physicsObject) {
this.physicsObject.endGrab(this.prevPos, this.vel);
this.physicsObject = null;
}
}
}
//-- mouse event
function mouse_event() {
scene.onPointerDown = function (event, pickInfo){
console.log("mouseDown pickInfo", pickInfo);
console.log("mouseDown pickInfo", pickInfo.hit);
if (pickInfo.hit == true) {
console.log("mouseDown pickInfo id", pickInfo.pickedMesh.id);
console.log("mouseDown pickInfo faceId", pickInfo.faceId);
gGrabber.start(pickInfo);
camera.detachControl(canvas);
}
}
scene.onPointerMove = function (event, pickInfo){
gGrabber.move(pickInfo);
}
scene.onPointerUp = function (event, pickInfo){
gGrabber.end(pickInfo);
camera.attachControl(canvas, true);
}
}
function update() {
mouse_event();
for (let step = 0; step < numSubsteps; step++) {
cloth.simulate();
}
cloth.updateMeshes();
gGrabber.increaseTime(dt);
}
let clothMesh = createClothMesh(cwidth, cheight, cdivisions);
let cloth = new Cloth(clothMesh, scene);
let gGrabber = new Grabber();
scene.registerBeforeRender(function () {
update();
});
engine.runRenderLoop(() => {
scene.render();
});
};
init();
</script>
</body>
</html>customMesh3D.js
// create Mesh Vertices and Indices
function createClothMesh(width, height, subdivision) {
let xi, yi, numXY, spaceX, spaceY;
numXY = subdivision;
spaceX = width / subdivision;
spaceY = height / subdivision;
let positions = [];
let px, py;
for (yi = 0; yi < (numXY + 1); yi++)
for (xi = 0; xi < (numXY + 1); xi++) {
let id = xi + yi * (numXY + 1);
px = (-numXY * 0.5 + xi) * spaceX;
py = ( numXY * 0.5 - yi) * spaceY;
pz = 0.0;
//console.log("id position", id, px, py,pz);
positions.push(px, py, pz);
}
//console.log("positions", positions);
let indices = [];
let id0, id1, id2, id3;
for (yi = 0; yi < numXY; yi++)
for (xi = 0; xi < numXY; xi++) {
id0 = yi * (numXY + 1) + xi;
id1 = (yi + 1) * (numXY + 1) + xi;
id2 = (yi + 1) * (numXY + 1) + xi + 1;
id3 = yi * (numXY + 1) + xi + 1;
indices.push(id0, id1, id2);
indices.push(id0, id2, id3);
}
return { name: "custom", vertices: positions, indices: indices}
}
Chrome Stable WebGPU
2022/07/05/
chrome stable version 103
WebGPU: version 2022/05/13
WGSL: version 2022/04/11
chrome の安定版（現在のバージョンは 103）でも、WebGPU を利用できます。
WebGPU Samples https://austin-eng.com/webgpu-samples/
このサンプルでは、Next.js, Node.js などのライブラリを使用しています。
これらのライブラリを使用しないで、WebGPU を利用するには、コマンドラインから --enable-unsafe-webgpu を付けて起動します。（”chrome://flags" では表示されません。）
mac の例 open -a "Google Chrome" --args --enable-unsafe-webgpu
chrome canary と同様に、この時点では GPUDevice の読込みが不安定になっています。requestDevice の read error となる場合は、再読み込みや再起動を行います。
変更点
| chrome 103 | chrome canary 105 |
|---|---|
| WebGPU | |
| GPUTextureFormat | |
| context.getPreferredFormat(adapter) | navigator.gpu.getPreferredCanvasFormat() |
| context.configure()内 | |
| compositingAlphaMode: "premultiplied" | alphaMode: "premultiplied" |
| WGSL | |
| Pipeline Stage Attributes | |
| @stage(vertex) | @vertex |
| @stage(fragment) | @frafment |
| @stage(compute) | @compute |
WebGPU Triangle and Instancing (WebGPU, WGSL update)
06/28/2022
chrome canary version 105
WebGPU, WGSL 06/17/2022
chrome canary について
（WebGPU を利用するには、chrome://flags/#enable-unsafe-webgpu
を設定する。）
この時点では GPUDevice の読み込みが不安定になっています。
requestDevice の read error となる場合は、再読み込みをします。
これを何回か繰り返すと、エラーは無いが表示が出来なくなる場合があります。
この場合は、chrome を再起動します。
[実行結果]
Triangle
Instancing
変更箇所
Triangle
webGPU
error
1 canvas.getContext("gpupresent") -> canvas.getContext("webgpu")
2 context.configureSwapChain() -> context.configure()
3 context.getPreferredFormat(adapter) ->
navigator.gpu.getPreferredCanvasFormat()
4 swapChain.getCurrentTexture().createView() ->
context.getCurrentTexture().createView()
5 render pass: endPass() -> end()
warning
1 CanvasConfiguration
alphaMode default change ->
alphaMode を使用するときは、”premultiplied" を設定
2 render pipeline lauout
layout の設定が必要 -> implicit pipeline layout の場合は、’auto' を設定
3 renderPass ColorAttachments
attachment -> view
loadValue -> clearValue
loadOp -> 要設定
storeOp -> 要設定
WGSL
attributes の変更
1 [[stage(vertex)]] -> @vertex 2 [[builtin(vertex_index)]] -> @builtin(vertex_index) [[builtin(position)]] -> @builtin(position) 3 [[stage(fragment)]] -> @fragment 4 [[location(0)]] -> @location(0)
Instancing
上記以外の変更箇所
WebGPU
1 compute pass dispatch() -> dispatchWorkgroups()
2 compute pass endPass() -> end()
WGSL
1 structure
[[block]] struct Uniforms {} -> [[block]] 不要
struct Uniforms {
proj_view : mat4x4<f32>
}
メンバー区切りの変更 ";" -> ","2 attributes
[[binding(0), group(0)]] -> @binding(0) @group(0) [[stage(compute)]] -> @compute @workgroup_size(64)
プログラム
triangle.html
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8">
<title>WebGPU Triangle WGSL</title>
</head>
<body>
<canvas id="webgpu-canvas" width="400" height="400"></canvas>
<script>
(async () => {
const adapter = await navigator.gpu.requestAdapter();
const device = await adapter.requestDevice();
const canvas = document.getElementById("webgpu-canvas");
const context = canvas.getContext("webgpu");
const wgslShaders = {
vertex: `
@vertex
fn main(@builtin(vertex_index) VertexIndex : u32)
-> @builtin(position) vec4<f32> {
var pos = array<vec2<f32>, 3>(
vec2<f32>(0.0, 0.5),
vec2<f32>(-0.5, -0.5),
vec2<f32>(0.5, -0.5));
return vec4<f32>(pos[VertexIndex], 0.0, 1.0);
}`,
fragment: `
@fragment
fn main() -> @location(0) vec4<f32> {
return vec4<f32>(1.0, 0.0, 0.0, 1.0);
}`,
};
const presentationFormat = navigator.gpu.getPreferredCanvasFormat();
context.configure({
device,
format: presentationFormat,
alphaMode: "premultiplied",
});
const pipeline = device.createRenderPipeline({
layout: "auto",
vertex: {
module: device.createShaderModule({
code: wgslShaders.vertex
}),
entryPoint: "main",
},
fragment: {
module: device.createShaderModule({
code: wgslShaders.fragment
}),
entryPoint: "main",
targets: [{
format: "bgra8unorm",
}],
},
primitive: {
topology: "triangle-list"
},
});
// render
requestAnimationFrame(function frame() {
const commandEncoder = device.createCommandEncoder();
const textureView = context.getCurrentTexture().createView();
const renderPassDesc = {
colorAttachments: [{
view: textureView,
clearValue: { r: 0.0, g: 0.0, b: 0.0, a: 1.0 },
loadOp: 'clear',
storeOp: 'store',
}]
};
// render pass
const renderPass = commandEncoder.beginRenderPass(renderPassDesc);
renderPass.setPipeline(pipeline);
renderPass.draw(3, 1, 0, 0);
renderPass.end();
device.queue.submit([commandEncoder.finish()]);
requestAnimationFrame(frame);
});
})();
</script>
</body>
</html>instancing.html
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8">
<title>WebGPU Instancing</title>
<script src="../libs/gl-matrix-min.js"></script>
<script src="../libs/webgl-util.js"></script>
<script src="../libs/utils.js"></script>
</head>
<body>
<canvas id="webgpu-canvas" width="500" height="400"></canvas>
<script>
(async () => {
const adapter = await navigator.gpu.requestAdapter();
const device = await adapter.requestDevice();
const canvas = document.getElementById("webgpu-canvas");
const context = canvas.getContext("webgpu");
const NUM_PARTICLES = 100;
const POINT_SIZE = 0.1;
//// camera
const eye = [0, 0, 4];
const center = [0, 0, 0];
const up = [0, 1, 0];
//// initial particle data
var particleData = new Float32Array(8 * NUM_PARTICLES);
for (let i = 0; i < particleData.length; i += 8)
{
//// position data
particleData[i] = Math.random() * 2.0 - 1.0;
particleData[i + 1] = Math.random() * 2.0 - 1.0;
particleData[i + 2] = Math.random() * 2.0 - 1.0;
particleData[i + 3] = 1;
//// velocity data
particleData[i + 4] = (Math.random() * 2.0 - 1.0) * 0.005;
particleData[i + 5] = (Math.random() * 2.0 - 1.0) * 0.005;
particleData[i + 6] = (Math.random() * 2.0 - 1.0) * 0.005;
particleData[i + 7] = 1;
}
const particleBuffer = device.createBuffer({
size: particleData.byteLength,
usage: GPUBufferUsage.VERTEX | GPUBufferUsage.COPY_DST | GPUBufferUsage.STORAGE
});
device.queue.writeBuffer(particleBuffer, 0, particleData);
//// WGSL Shader (vertex, fragnebt, compute)
const wgslShaders = {
vertex:`
struct Uniforms {
proj_view : mat4x4<f32>
};
@binding(0) @group(0) var<uniform> uniforms : Uniforms;
@vertex
fn main(@location(0) vertexPosition : vec3<f32>,
@location(2) position : vec3<f32>)
-> @builtin(position) vec4<f32>{
var scale : f32 = ${POINT_SIZE};
var scaleMTX : mat4x4<f32> = mat4x4<f32>(
vec4<f32>(scale, 0.0, 0.0, 0.0),
vec4<f32>(0.0, scale, 0.0, 0.0),
vec4<f32>(0.0, 0.0, scale, 0.0),
vec4<f32>(position, 1.0)
);
return uniforms.proj_view * scaleMTX * vec4<f32>(vertexPosition, 1.0);
}
`,
fragment:`
@fragment
fn main() -> @location(0) vec4<f32> {
return vec4<f32>(1.0, 0.0, 0.0, 1.0);
}
`,
compute:`
struct Particle {
pos : vec3<f32>,
vel : vec3<f32>,
};
struct Particles {
particles : array<Particle, ${NUM_PARTICLES}>
};
@binding(0) @group(0) var<storage, read_write> particle : Particles;
@compute @workgroup_size(64)
fn main(@builtin(global_invocation_id) GlobalInvocationID : vec3<u32>) {
var index: u32 = GlobalInvocationID.x;
var position: vec3<f32> = particle.particles[index].pos;
var velocity: vec3<f32> = particle.particles[index].vel;
var new_position: vec3<f32> = position + velocity;
if (abs(new_position.x) >= 1.0) {
particle.particles[index].vel.x = - velocity.x;
}
if (abs(new_position.y) >= 1.0) {
particle.particles[index].vel.y = - velocity.y;
}
if (abs(new_position.z) >= 1.0) {
particle.particles[index].vel.z = - velocity.z;
}
particle.particles[index].pos = new_position;
}
`
};
const computePipeline = device.createComputePipeline({
layout: "auto",
compute: {
module: device.createShaderModule({
code: wgslShaders.compute}),
entryPoint: "main"
}
});
const computeBindGroup = device.createBindGroup({
layout: computePipeline.getBindGroupLayout(0),
entries: [
{
binding: 0,
resource: {
buffer: particleBuffer
}
},
]
});
const cubeData = utils.createCube();
const numVertices = cubeData.positions.length / 3;
const vertexBuffer = device.createBuffer({
size: cubeData.positions.byteLength,
usage: GPUBufferUsage.VERTEX | GPUBufferUsage.COPY_DST
});
device.queue.writeBuffer(vertexBuffer, 0, cubeData.positions);
//// Setup render outputs
const presentationFormat = navigator.gpu.getPreferredCanvasFormat();
context.configure({
device,
format: presentationFormat,
alphaMode: "premultiplied",
});
var pipeline_point = device.createRenderPipeline({
layout: "auto",
vertex: {
module: device.createShaderModule({
code: wgslShaders.vertex}),
entryPoint: "main",
buffers: [
{
arrayStride: 12,
attributes: [
{
shaderLocation: 0,
format: "float32x3",
offset: 0
}]
},
{
arrayStride: 16 * 2,
stepMode: "instance",
attributes: [{
shaderLocation: 2,
format: "float32x3",
offset: 0
}]
}
],
},
fragment: {
module: device.createShaderModule({
code: wgslShaders.fragment}),
entryPoint: "main",
targets: [{
format: "bgra8unorm",
}],
},
primitive: {
topology: "triangle-list"
},
});
//// Setup renderPassDesc
var renderPassDesc = {
colorAttachments: [{
view: undefined,
clearValue: { r: 0.2, g: 0.2, b: 0.2, a: 1.0 },
loadOp: 'clear',
storeOp: 'store',
}],
};
//// Create uniform buffer
var viewParamsBuffer = device.createBuffer({
size: 16 * 4,
usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST
});
//// Create bind group (setup uniform buffer)
var bindGroup = device.createBindGroup({
layout: pipeline_point.getBindGroupLayout(0),
entries: [
{
binding: 0,
resource: {
buffer: viewParamsBuffer
}
}
]
});
//// Create arcball camera and view projection matrix
var camera = new ArcballCamera(eye, center, up,
0.5, [canvas.width, canvas.height]);
var projection = mat4.perspective(mat4.create(), 50 * Math.PI / 180.0,
canvas.width / canvas.height, 0.1, 100);
//// projection_view matrix
var projView = mat4.create();
//// Controller
var controller = new Controller();
controller.mousemove = function(prev, cur, evt) {
if (evt.buttons == 1) {
camera.rotate(prev, cur);
} else if (evt.buttons == 2) {
camera.pan([cur[0] - prev[0], prev[1] - cur[1]]);
}
};
controller.wheel = function(amt) { camera.zoom(amt * 3.0); };
controller.registerForCanvas(canvas);
//// render
requestAnimationFrame(function frame() {
const commandEncoder = device.createCommandEncoder();
//// compute pass
var computePass = commandEncoder.beginComputePass();
computePass.setPipeline(computePipeline);
computePass.setBindGroup(0, computeBindGroup);
computePass.dispatchWorkgroups(NUM_PARTICLES);
computePass.end();
//// render pass setting
renderPassDesc.colorAttachments[0].view =
context.getCurrentTexture().createView();
projView = mat4.mul(projView, projection, camera.camera);
device.queue.writeBuffer(viewParamsBuffer, 0, projView);
const renderPass = commandEncoder.beginRenderPass(renderPassDesc);
renderPass.setPipeline(pipeline_point);
renderPass.setVertexBuffer(0, vertexBuffer);
renderPass.setVertexBuffer(1, particleBuffer);
renderPass.setBindGroup(0, bindGroup);
renderPass.draw(numVertices, NUM_PARTICLES, 0, 0);
renderPass.end();
device.queue.submit([commandEncoder.finish()]);
requestAnimationFrame(frame);
});
})();
</script>
</body>
</html>プログラムには、以下のライブラリを使用しています。
１ point(cubeで表現)を描画するためのcubeプログラム utils.js
WebGPU Examples https://github.com/tsherif/webgpu-examples
２ cameraプログラム webgl-util.min.js
WebGPU Experiments https://github.com/Twinklebear/webgpu-experiments
３ matrix and vector計算ライブラリ gl-matrix.js
https://github.com/toji/gl-matrix
物理シミュレーション 剛体の衝突 4 PositionBased ( ShapeMatching )
PositionBased ( ShapeMatching ) による、剛体の衝突シミュレーションです。
[ 実行結果]
参考にした文献とサイト
Meshless Deformations Based on Shape Matching
https://matthias-research.github.io/pages/publications/publications.html
藤澤 誠著「CGのための物理シミュレーションの基礎」( マイナビ )
実装用に参考にしたサイト
Physical animation with shape matching (Claudio Esperança)
https://observablehq.com/@esperanc/physical-animation-with-shape-matching
右上にある […] アイコンの Export からコードをダウンロードできる。
69eb2ec167e70941@1046.js がメインのプログラムである。
シェイプマッチング (Shape Matching) 法
シェイプマッチング法は、剛体を頂点座標で表し、剛体が変形した際の頂点の移動を、剛体の元の形状を保つように移動させる方法である。
（図は、「CGのための物理シミュレーションの基礎」より作成）
剛体の初期状態の頂点座標を 、変形後の頂点座標を
とすると、
は、回転行列
、変形の行列
、平行移動ベクトル
を使って、以下のように表される。
( ：初期状態の剛体の重心座標 )
最終的に移動する位置 (目標位置、goal position)は、剛体の元の形状を保つように移動させるため、回転行列と平行移動ベクトルを使って表す。
と
は、
と
から推定するが、これらは
と
の差が最小となるように、つまり以下の式が最小となるように決める。
(
：変形後の剛体の重心座標 ) であるので、
ここで、 、
、
(回転の行列に変形の行列を含めて拡張した行列)である。
を
の成分
(
、
：次元) で微分する。
より
は、
の右極分解から得ることができる。( 極分解については下参照）
これから、剛体の移動後の目標位置 が求められる。
極分解
任意の正方行列は、対称行列と直交行列の積に分解される。これが極分解である。
：正方行列、
：直交行列、
、
：対称行列 として、
右極分解
左極分解
、
は対称行列 なので、
は直交行列なので、
(
：単位行列）
である。これより、右極分解の場合、 と
が次のように求まる。
ここで、行列の平方根が必要となるが、これは以下のように計算することができる。
１）ケーリー・ハミルトン(Cayley-Hamilton)の定理を使う
The Square Roots of 2 x 2 Matrices
(https://www.maa.org/sites/default/files/pdf/cms_upload/Square_Roots-Sullivan13884.pdf) 参照
と
は２次正方行列であるとする。
であるとき、
と
にケーリー・ハミルトンの定理を用いると、
は次にように求まる。
２）対角化の方法を使う
行列 の固有値は
、固有ベクトルを
、
とすると、対角化行列
とその逆行列
は、次のようになる。
、
、
対角行列 は、
である。対角行列の平方根はすぐに求まる。
これより、行列 の平方根の行列
は、
、
に注意して式を変形すると、１）の
と同じ式が得られる。
mat2.js ライブラリーでの行列の平方根 (mat2sqrt) の計算は、この式を用いている。
シミュレーションの方法
１ 時間積分
（WorldShapeMatching.js class World, step メソッド）
1) 重力による速度の更新と位置、回転角の更新
2) 壁との衝突後に移動させる剛体の頂点座標(goal position, 目標位置)の計算
(class World, projectWalls メソッド
+ BodyShapeMatching.js class Body, shapeMatch メソッド )
3) 剛体同士の衝突後に移動させる剛体の頂点座標(goal position)の計算
(class World, projectCollision メソッド
+ class Body, shapeMatch メソッド )
4) 剛体の回転角と位置の更新
(class Body, updateFromShape メソッド )
2), 3) で得られた頂点座標から、剛体全体の回転角と位置を更新する。
通常のシェイプマッチングでは、剛性(stiffness)の効果を各頂点に施すが、
ここでは剛体全体の回転角と位置に施している。
(shapeMatchAttenuation の部分)
5) 剛体の速度と角速度の更新
２ 壁との衝突
(class World, projectWalls メソッド + class Body, shapeMatch メソッド )
剛体と壁との衝突
vertex で衝突する場合
edge で衝突する場合
剛体が壁と衝突した状態（上図変形の左側の図）を初期状態とし、その時の頂点座標を とする。変形後の状態（上図変形の右側の図）の頂点座標 を図のように、めり込んだ頂点を壁の表面まで押し返したようにとる。
プログラムでは、 、
を、それぞれ body.worldShape 、body.collisionShape としている。
を計算しているのが、world.projectWalls(body) メソッドである。
これらの値から、実際に移動させる剛体の回転行列 と平行移動行列
(平行移動ベクトルを行列にしている)を計算し、頂点座標
を求める。この部分が body.shapeMatch() メソッドである。
body.shapeMatch() メソッド
M = shapeMatch(this.worldShape, this.collisionShape)
-> と
の平行移動も含めた行列として
と
を計算し、
としている。
this.collisionShape = MT.mulVectors(this.worldShape, 1)
-> ( MT = M ) を計算し、
を新たに body.collisionShape
としている。 の計算をしている。
３ 剛体同士の衝突
(class World, projectCollision メソッド + class Body, shapeMatch メソッド )
剛体同士の衝突
vertex と edge で衝突
edge と edge で衝突
衝突判定には、前回と同様に、sap法と gjk-epa法を用いている。
剛体同士が衝突した状態を初期状態とし、その時の頂点座標が である。変形後の状態の頂点座標
は、gjk-epa から得られる貫通の深さの半分を、衝突している頂点を衝突の法線方向に押し戻したようにとる。
壁との衝突と同様に、この と
から、回転行列
と平行移動行列
を計算し、移動する頂点座標
を求める。
プログラム
collision-shapeMatching.html
<!DOCTYPE html>
<meta charset="utf-8">
<title>BJS Physical Animation ( ShapeMatching )</title>
<script src="https://preview.babylonjs.com/babylon.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.5/require.min.js"></script>
<body>
<canvas id="canvas" width="600" height="430" style="border: 1px solid gray;"></canvas>
<script>
function init() {
const canvas = document.getElementById("canvas");
const engine = new BABYLON.Engine(canvas);
var width = canvas.width;
var height = canvas.height;
var scene = new BABYLON.Scene(engine);
scene.clearColor = new BABYLON.Color3(1.0, 1.0, 1.0);
var camera = new BABYLON.ArcRotateCamera("camera1",
3 * Math.PI / 2, Math.PI / 2, 25, new BABYLON.Vector3(0, 8, 0), scene);
camera.attachControl(canvas, true);
var light1 = new BABYLON.HemisphericLight("light1", new BABYLON.Vector3(0, 1, 0), scene);
light1.intensity = 0.8;
var light2 = new BABYLON.HemisphericLight("light2", new BABYLON.Vector3(0, -1, 0), scene);
light2.intensity = 0.5;
require(['./WorldShapeMatching', './BodyShapeMatching', './math'], function(World, Body, math) {
let walls = [];
let wall_angle1 = -0.3;
let wall_angle2 = 0.5;
walls.push(new Wall(new math.Vector2( 0, 0),
new math.Vector2(-Math.sin(wall_angle1), Math.cos(wall_angle1))));
walls.push(new Wall(new math.Vector2( 0, 0),
new math.Vector2(-Math.sin(wall_angle2), Math.cos(wall_angle2))));
let bodies = [];
// body 0
bodies.push(
new Body(
createPolygonPoints(new math.Vector2(0, 0), 2, 4),
new math.Vector2(7, 14), // world position(pos)
Math.PI / 4, // angle
1 // mass
)
);
bodies[0].vel = new math.Vector2(0, 0);
bodies[0].angVel = 0.0;
// body 1
bodies.push(
new Body(
createPolygonPoints(new math.Vector2(0, 0), 1, 4),
new math.Vector2(-12, 14),
0,
1
)
);
bodies[1].vel = new math.Vector2(0, 0);
bodies[1].angVel = 0.0;
// body 2
bodies.push(
new Body(
createPolygonPoints(new math.Vector2(0, 0), 1.5, 4),
new math.Vector2(-8, 14),
0,
1
)
);
bodies[2].vel = new math.Vector2(2, 0);
bodies[2].angVel = 0.0;
//
console.log("bodies", Array.from(bodies));
var d, p1, p2;
for (let i = 0; i < walls.length; i++) {
d = new math.Vector2(400 * walls[i].normal.y, -400 * walls[i].normal.x);
p1 = walls[i].position.addV(d);
p2 = walls[i].position.subV(d);
let wallpath = [
new BABYLON.Vector3(p1.x, p1.y, 0),
new BABYLON.Vector3(p2.x, p2.y, 0)
];
let wallLines = BABYLON.Mesh.CreateLines("lines" + i, wallpath, scene);
wallLines.color = BABYLON.Color3.Black();
}
var i, j;
var objects, object;
var colors = [
BABYLON.Color3.Blue(),
BABYLON.Color3.Green(),
BABYLON.Color3.Black(),
];
var objectMesh = [];
objects = bodies;
for (j = 0; j < objects.length; ++j) {
object = objects[j];
var worldVertices = object.worldShape;
console.log("object index", j, worldVertices);
const abodyMat = new BABYLON.StandardMaterial("abodyMat" + j, scene);
abodyMat.diffuseColor = colors[j];
var points = [];
for (i = 0; i < worldVertices.length; ++i) {
const apoint = BABYLON.MeshBuilder.CreateSphere("apoint" + i,
{diameter: 0.3, segments: 8}, scene);
apoint.position.x = worldVertices[i].x;
apoint.position.y = worldVertices[i].y;
apoint.position.z = 0;
apoint.material = abodyMat;
points.push(apoint);
}
var path = [];
for (i = 0; i < worldVertices.length; i++)
path.push(
new BABYLON.Vector3(worldVertices[i].x, worldVertices[i].y, 0)
);
path.push(new BABYLON.Vector3(worldVertices[0].x, worldVertices[0].y, 0));
const options = {points: path, updatable: true};
var aLine = BABYLON.MeshBuilder.CreateLines("aLine" + j, options, scene);
aLine.color = colors[j];
objectMesh.push({points: points, line: aLine});
}
function update_mesh() {
var i, j
for (j = 0; j < objects.length; j++) {
var object = objects[j];
var worldVertices = object.worldShape;
//// points
var apoint = objectMesh[j].points;
for (i = 0; i < worldVertices.length; ++i) {
apoint[i].position.x = worldVertices[i].x;
apoint[i].position.y = worldVertices[i].y;
apoint[i].position.z = 0.0;
}
//// line
var aline = objectMesh[j].line;
var path = [];
for (i = 0; i < worldVertices.length; i++)
path.push(
new BABYLON.Vector3(worldVertices[i].x, worldVertices[i].y, 0)
);
path.push(new BABYLON.Vector3(worldVertices[0].x, worldVertices[0].y, 0));
var options = {points: path, instance: aline};
aline = BABYLON.MeshBuilder.CreateLines("aline", options);
}
}
let world = new World();
world.walls = walls;
world.bodies = bodies;
scene.registerBeforeRender(function () {
world.step();
update_mesh();
});
function createPolygonPoints(center = new math.Vector2(0, 0), radius = 1, nsides=3) {
let poly = [];
let delta = Math.PI * 2 / nsides;
for (let i = 0; i < nsides; i++) {
poly.push(
new math.Vector2(center.x+radius*Math.cos(delta*i),
center.y+radius*Math.sin(delta*i))
)
}
// inertia
let [a, b] = [(1 + Math.cos(delta)) / 2, (1 - Math.cos(delta)) / 2];
let inertia = radius * radius / 6 * ( 3*(a**3 + b**3) + 9*a*b - 2*b );
return {shape:poly, inertia:inertia};
}
function Wall(position, normal) {
this.position = position;
this.normal = normal;
this.dist = function(position) {
return this.normal.dot(position.subV(this.position));
};
}
});
engine.runRenderLoop(() => {
scene.render();
});
}
init();
</script>
WorldShapeMatching.js
define(['./math', './gjk_epa'], function(math, Gjk) {
function centerOfMass(pts) {
let sum = pts.reduce((a,b) => new math.Vector2(a.x+b.x, a.y+b.y));
return new math.Vector2(sum.x/pts.length, sum.y/pts.length);
}
function supportEdge(pts, v) {
let maxdot = Number.NEGATIVE_INFINITY;
let maxdot2 = Number.NEGATIVE_INFINITY;
let best, best2;
for (let p of pts) {
let dot = p.dot(v);
if (dot > maxdot) {
[maxdot2, best2] = [maxdot, best];
[maxdot, best] = [dot, p];
} else if (dot > maxdot2) {
[maxdot2, best2] = [dot, p];
}
}
if (Math.abs(maxdot - maxdot2) < 0.01) return [best2, best];
return [best];
}
function closestSegmentPoint(p, q, r) {
let qr = r.subV(q);
let s = qr.getLength();
if (s < 0.00001) return q; // Degenerate line segment
let v = qr.normalized();
let u = p.subV(q);
let d = u.dot(v);
if (d < 0) return q;
if (d > s) return r;
return lerp(q, r, d/s);
}
function lerp(a, b, t) {
var ax = a.x,
ay = a.y;
out_x = ax + t * (b.x - ax);
out_y = ay + t * (b.y - ay);
return new math.Vector2(out_x, out_y);
}
function sap(polygons, v) {
let n = polygons.length;
let pairs = [];
let proj = [];
polygons.forEach((poly, i) => {
let min = Number.POSITIVE_INFINITY;
let max = Number.NEGATIVE_INFINITY;
for (let p of poly) {
let dot = p.dot(v);
min = Math.min(min, dot);
max = Math.max(max, dot);
}
proj.push([min, i], [max, i]);
});
proj.sort((a, b) => a[0] - b[0]);
let inside = new Set();
for (let [_, i] of proj) {
if (inside.has(i)) {
inside.delete(i);
} else {
pairs[i] = [];
for (let j of inside) {
if (i < j) pairs[j].push(i);
else pairs[i].push(j);
}
inside.add(i);
}
}
return pairs;
}
const g = 9.8;
const dt = 0.02;
const wall_restitution = 0.9;
const wall_friction = 0.02;
class World {
constructor() {
this.nIterations = 4;
this.bodies = [];
this.walls = [];
}
step() {
//// Apply external force
this.bodies.forEach(body => {
body.vel = body.vel.addV(new math.Vector2(0, -g * dt));
body.pos = body.oldPos.addV(body.vel.mulS(dt));
body.ang = body.oldAng + body.angVel * dt;
});
//// collisions
for (let irelax = 0; irelax < this.nIterations; irelax++) {
//// walls
this.bodies.forEach(body =>
this.projectWalls(body)
);
this.bodies.forEach(body => body.shapeMatch());
//// bodies
let pairs = sap(this.bodies.map(body => body.worldShape), new math.Vector2(1, 0));
this.bodies.forEach((body, i) => {
pairs[i].forEach(j => this.projectCollision(body, this.bodies[j]), irelax);
});
this.bodies.forEach(body => body.shapeMatch());
} // for irelax
this.bodies.forEach(body => {
if (body.collisionShape) body.updateFromShape(body.collisionShape);
});
this.bodies.forEach(body => {
body.vel = body.pos.subV(body.oldPos).divS(dt);
body.oldPos = body.pos;
body.angVel = (body.ang - body.oldAng) / dt;
body.oldAng = body.ang;
});
}
//// projectWalls
projectWalls(body) {
let projected = false;
let poly = body.collisionShape || body.worldShape.slice();
poly.forEach((p, i) => {
for (let wall of this.walls) {
let dist = wall.dist(p);
if (dist <= 0) {
var velocity_normal = wall.normal.mulS(wall.normal.dot(body.vel));
var velocity_tangent = body.vel.subV(velocity_normal);
velocity_normal = velocity_normal.mulS(-wall_restitution);
velocity_tangent = velocity_tangent.mulS(-wall_friction);
p = p.addV(velocity_normal.addV(velocity_tangent).mulS(dt));
dist = wall.dist(p);
if (dist < 0)
p = p.addV(wall.normal.mulS(-dist));
projected = true;
}
}
poly[i] = p;
});
if (projected) body.collisionShape = poly;
}
//// projectCollision
projectCollision(bodyA, bodyB, gravityBias = 0) {
let a = bodyA.collisionShape || bodyA.worldShape.slice();
let b = bodyB.collisionShape || bodyB.worldShape.slice();
let hit = Gjk.gjk(a, b);
if (hit) {
let { p, q, n } = Gjk.epa(a, b, ...hit);
let aPoints = supportEdge(a, n);
let bPoints = supportEdge(b, n.mulS(-1));
let [massA, massB] = [bodyA.mass, bodyB.mass];
if (gravityBias) {
if (bodyA.pos.y > bodyB.pos.y) massA += massB * gravityBias;
else massB += massA * gravityBias;
}
let aContact, bContact, penetration;
let aContactDisplaced1, aContactDisplaced2, bContactDisplaced1, bContactDisplaced2;
if (aPoints.length + bPoints.length == 4) {
//// Edge-edge collision
let center = centerOfMass([...aPoints, ...bPoints]);
aContact = closestSegmentPoint(center, ...aPoints);
bContact = closestSegmentPoint(center, ...bPoints);
penetration = aContact.subV(bContact).getLength();
aContactDisplaced1 =
aPoints[0].addV(
n.mulS((-penetration * massA) / (massA + massB))
);
aContactDisplaced2 =
aPoints[1].addV(
n.mulS((-penetration * massA) / (massA + massB))
);
bContactDisplaced1 =
bPoints[0].addV(
n.mulS((penetration * massB) / (massA + massB))
);
bContactDisplaced2 =
bPoints[1].addV(
n.mulS((penetration * massB) / (massA + massB))
);
a.splice(a.lastIndexOf(aPoints[0]), 1, aContactDisplaced1);
a.splice(a.lastIndexOf(aPoints[1]), 1, aContactDisplaced2);
b.splice(b.lastIndexOf(bPoints[0]), 1, bContactDisplaced1);
b.splice(b.lastIndexOf(bPoints[1]), 1, bContactDisplaced2);
} else {
//// Vertex-edge collision
if (aPoints.length + bPoints.length != 3) {
throw "Weird collision";
}
if (aPoints.length == 2) {
aContact = closestSegmentPoint(bPoints[0], ...aPoints);
penetration = aContact.subV(bPoints[0]).getLength();
aContactDisplaced1 =
aPoints[0].addV(
n.mulS((-penetration * massA) / (massA + massB))
);
aContactDisplaced2 =
aPoints[1].addV(
n.mulS((-penetration * massA) / (massA + massB))
);
a.splice(a.lastIndexOf(aPoints[0]), 1, aContactDisplaced1);
a.splice(a.lastIndexOf(aPoints[1]), 1, aContactDisplaced2);
bContactDisplaced1 =
bPoints[0].addV(
n.mulS((penetration * massB) / (massA + massB))
);
b.splice(b.lastIndexOf(bPoints[0]), 1, bContactDisplaced1);
} else {
//// bPoints.length == 2!
bContact = closestSegmentPoint(aPoints[0], ...bPoints);
penetration = aPoints[0].subV(bContact).getLength();
bContactDisplaced1 =
bPoints[0].addV(
n.mulS((penetration * massB) / (massA + massB))
);
bContactDisplaced2 =
bPoints[1].addV(
n.mulS((penetration * massB) / (massA + massB))
);
b.splice(b.lastIndexOf(bPoints[0]), 1, bContactDisplaced1);
b.splice(b.lastIndexOf(bPoints[1]), 1, bContactDisplaced2);
aContactDisplaced1 =
aPoints[0].addV(
n.mulS((-penetration * massA) / (massA + massB))
);
a.splice(a.lastIndexOf(aPoints[0]), 1, aContactDisplaced1);
}
}
bodyA.collisionShape = a;
bodyB.collisionShape = b;
}
}
}
return World;
});
BodyShapeMatching.js
define(['./math', './mat2'], function(math, mat2) {
function centerOfMass(pts) {
let sum = pts.reduce((a,b) => new math.Vector2(a.x+b.x, a.y+b.y));
return new math.Vector2(sum.x/pts.length, sum.y/pts.length);
}
function lerp(a, b, t) {
var ax = a.x,
ay = a.y;
out_x = ax + t * (b.x - ax);
out_y = ay + t * (b.y - ay);
return new math.Vector2(out_x, out_y);
}
function shapeMatch(srcPoints, dstPoints) {
//// srcPoints
let srcCenter = centerOfMass(srcPoints);
let srcVectors = srcPoints.map(p => p.subV(srcCenter));
//// dstPoints
let dstCenter = centerOfMass(dstPoints);
let dstVectors = dstPoints.map (p => p.subV(dstCenter));
//// Compute rotation and compose with the two translations
let [[a,b],[c,d]] = shapeMatchRotation(srcVectors,dstVectors);
let rot = new math.Matrix3();
rot = rot.fromValues(a, b, 0, c, d, 0, dstCenter.x, dstCenter.y, 1);
let trans = new math.Matrix3();
trans = trans.fromValues(1, 0, 0, 0, 1, 0, -srcCenter.x, -srcCenter.y, 1);
let transformMatrix = trans.mulMatrix(rot);
return transformMatrix;
}
function shapeMatchRotation(srcVectors, dstVectors) {
//// srcVectors, dstVectors: convert element Vector2 -> element Array
let srcVectors2 = srcVectors.map(p => [p.x, p.y]);
let dstVectors2 = dstVectors.map(p => [p.x, p.y]);
//// function that computes p x q^T
let pqT = (p,q) => [[p[0]*q[0],p[0]*q[1]],
[p[1]*q[0],p[1]*q[1]]];
let Apq = srcVectors2.map((p,i) => pqT(p,dstVectors2[i])).reduce(mat2.mat2sum)
let ApqTxApq = mat2.mat2mul(mat2.mat2transpose(Apq),Apq);
let S = mat2.mat2sqrt(ApqTxApq)
let Sinv = mat2.mat2inv(S);
return mat2.mat2mul(Apq,Sinv)
}
const shapeMatchAttenuation = 0.98;
class Body {
constructor(shape, pos = new Vector2(0, 0), ang = 0, mass = 1) {
console.log("ceate Bodies", shape);
this.mass = mass;
this.pos = this.oldPos = pos;
this.ang = this.oldAng = ang;
this.vel = new math.Vector2(0, 0);
this.angVel = 0;
let center = centerOfMass(shape.shape);
this.shape = shape.shape.map(p => new math.Vector2(p.x - center.x, p.y - center.y));
this.collisionShape = null; // Non-null if collision was detected
}
//// Returns the shape points in world space
get worldShape() {
let [c, s] = [Math.cos(this.ang), Math.sin(this.ang)];
let transf = new math.Matrix3();
transf = transf.fromValues(c, -s, this.pos.x, s, c, this.pos.y, 0, 0, 0);
this.curWorldShape = transf.mulVectors(this.shape, 1);
return this.curWorldShape;
}
//// ShapeMatching
shapeMatch() {
if (this.collisionShape) {
let M = shapeMatch(this.worldShape, this.collisionShape);
let MT = new math.Matrix3();
MT = MT.fromValues(
M.m00, M.m10, M.m20,
M.m01, M.m11, M.m21,
M.m02, M.m12, M.m22
);
this.collisionShape = MT.mulVectors(this.worldShape, 1);
}
}
//// Updates rotation and position
updateFromShape(shape) {
let center = centerOfMass(shape);
//// Rotation component
let rot = shapeMatchRotation(
this.worldShape.map(p => p.subV(this.pos)),
shape.map(p => p.subV(center))
);
if (!Number.isNaN(rot[0][0])) {
//// Avoid degenerate projections
let dang = Math.atan2(rot[0][1], rot[0][0]);
this.ang += dang * shapeMatchAttenuation;
}
//// Translation component
this.pos = lerp(this.pos, center, shapeMatchAttenuation);
//// new worldShape
this.curWorldShape = shape;
this.collisionShape = null;
}
}
return Body;
});
mat2.js
// 2x2 matrix utilities
define( function() {
return {
mat2mul: function(a,b) {
let prod = (i,j) => a[i][0]*b[0][j]+a[i][1]*b[1][j];
return [[prod(0,0),prod(0,1)],
[prod(1,0),prod(1,1)]]
},
mat2sum: function(m1,m2) {
return [[m1[0][0]+m2[0][0],m1[0][1]+m2[0][1]],
[m1[1][0]+m2[1][0],m1[1][1]+m2[1][1]]]
},
mat2transpose: function(m) {
return [[m[0][0],m[1][0]],
[m[0][1],m[1][1]]]
},
mat2det: function([[a,b],[c,d]]) {
return a*d-b*c
},
mat2sqrt: function([[a,b],[c,d]]) {
let s = Math.sqrt(a*d-b*c);
let t = Math.sqrt(a+d+2*s);
return [[(a+s)/t,b/t],[c/t,(d+s)/t]]
},
mat2inv: function([[a,b],[c,d]]) {
let det = a*d-b*c;
return [[d/det,-c/det],[-b/det,a/det]]
}
}
});gjk-epa.js、math.js は「剛体の衝突 １」に掲載してある。
物理シミュレーション 剛体の衝突 ３ PositionBased ( Constraint )
PositionBased(Constraint) による、剛体の衝突シミュレーションです。
[ 実行結果]
２Dシミュレーション
メッシュの描画には、 BabylonJS を使用しています。
参考にした文献とサイト
Position Based Dynamics
https://matthias-research.github.io/pages/publications/publications.html
Modeling and Solving Constraints
https://ubm-twvideo01.s3.amazonaws.com/o1/vault/gdc09/slides/04-GDC09_Catto_Erin_Solver.pdf
物理エンジンの作り方
http://cedec.cesa.or.jp/2008/archives/file/pg07.pdf
物理エンジンの作り方 その２
http://cedec.cesa.or.jp/2009/ssn_archive/pdf/sep3rd/PG75.pdf
実装用に参考にしたサイト
Inequality Constraints
http://myselph.de/gamePhysics/inequalityConstraints.html
2D Rigid Body Physics Engine
https://github.com/Sopiro/Physics
(実際に参照したのは、古いバージョン:(2021/11/06 ダウンロード)
の game.js。新バージョンの contact.js に相当。)
シミュレーションの方法
１ 時間積分
（WorldConstraint.js class World, step メソッド）
1) 重力による速度の更新
2) 壁との衝突による速度と角速度の更新
(projectWalls メソッド)
衝突時の拘束を解いて、速度と角速度を更新する。
3) 剛体同士の衝突による速度と角速度の更新
(projectCollision メソッド)
壁との衝突と同様に、衝突時の拘束を解いて、速度と角速度を更新する。
4) 位置と角度の更新
新しく得られた速度と角速度を数値積分して、位置と角度を更新する。
２ 壁との衝突
（WorldConstraint.js class World, projectWalls メソッド）
剛体と壁との衝突の様子
剛体と壁との衝突点は、前回の ForceBased の場合と同様にして、知ることが
できる。
衝突前後の速度と角速度を、それぞれ 、
、
、
とすると、
衝突点の速度は、並進の速度と回転による速度の合成速度で、次のようになる。
衝突前：
衝突後：
：剛体の重心から衝突点へのベクトル（重心-衝突点ベクトル）
貫通なし拘束は、速度を用いて表される。
(
は衝突面の法線方向の単位ベクトル )
一般化速度 ベクトル（T：転置）を用いると、拘束の式は次のようになる。
はヤコビアンで、
拘束の式に、反発の効果を bias として取り入れる。
反発係数(restitution)を とすると
となるので、拘束の式を次のようにする。
,
次に拘束力を取り入れる。
拘束力 は、運動量の変化で表される。
は質量マトリックスで、
（
：単位マトリックス、
：慣性モーメントマトリックス）である。
は力の作用時間である。
拘束力の方向は となるので、大きさを
（未定乗数）として、
を
の中に含めると、次のようになる。
新たな はインパルスに相当する。
拘束の式は、次のようになる。
この式から を求める。
具体的には、以下の式を繰返し法を用いて解いていく。
,
このとき、繰返し毎にインパルス()を加算し、それがゼロとなるようにクランプ(clamp)する。
これで、 (拘束力が作用しない) 、
(剛体が壁から離れた状態)となり、拘束が解けたことになる。
実際のシミュレーションでは、安定性のため、次のような項を取り入れる。
１つ目は、貫通に関する安定項である。
：貫通の長さ、
：貫通に対するslop(遊び)で適当に設定、
：係数で適当に設定
２つ目は、反発項に関する修正である。
：反発に対するslopで適当に設定
bias項 を
にする。
摩擦がある場合には、摩擦による拘束を追加する。
この場合の拘束は、次のように表される。
(
は衝突面の接線方向の単位ベクトル )
一般化速度 を用いると、拘束の式は次のようになる。
,
次に拘束力を取り入れる。摩擦力は、
：衝突面での垂直抗力
：摩擦係数
である。インパルスで表すと、次のようになる。
は上記の法線方向のインパルスである。加算したインパルスを用いる。
拘束の式は、インパルスを取り入れると次のようになる。
この式とインパルスの式から を求める。
具体的には、以下の式を、上述した法線方向の拘束の式の後に追加して、
繰返し法を用いて解いていく。
このとき、繰返し毎にインパルス()を加算し、それが
と
の間の値に
なるようにクランプ(clamp)する。これで拘束が解けたことになる。
繰返しの最後に得られた速度 が、拘束によって更新された速度となる。
３ 剛体同士の衝突
（WorldConstraint.js class World, projectCollision メソッド）
剛体同士の衝突の様子
衝突判定には、ForceBased と同様な方法、sap法（sap() 関数）と gjk-epa 法（gjk-epa.js）を用いている。
gjk-epa から、貫通の深さ、衝突の方向（衝突面の法線方向）、衝突点の座標等の情報が得られる。
衝突前後の剛体A、Bの速度と角速度を、それぞれ 、
、
、
、
、
、
、
とすると、
剛体A、Bの衝突点の速度は、次のようになる。
剛体A
衝突前：
衝突後：
剛体B
衝突前：
衝突後：
、
：剛体A、Bの重心から衝突点へのベクトル
貫通なし拘束は、剛体の相対速度を用いて表される。
(
は衝突面の法線方向の単位ベクトル )
一般化速度 ベクトルを用いると、拘束の式は次のようになる。
拘束力をインパルス( )として取り入れると、拘束の式は次のようになる。
：質量
、
、慣性モーメント
、
の質量マトリックス
：bias
摩擦がある場合の拘束の式は、次のようになる。
(
は衝突面の接線方向の単位ベクトル )
（
：摩擦係数）
拘束は、壁との衝突の場合と同様に、繰返し法で解いていく。
プログラム
collision-constraint.html
<!-- collision-constraint.html <= BJS-gjk-constraint-0-1.html -->
<!DOCTYPE html>
<meta charset="utf-8">
<title>BJS Physical Animation ( constraint )</title>
<script src="https://preview.babylonjs.com/babylon.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.5/require.min.js"></script>
<body>
<canvas id="canvas" width="600" height="430" style="border: 1px solid gray;"></canvas>
<script>
function init() {
const canvas = document.getElementById("canvas");
const engine = new BABYLON.Engine(canvas);
var width = canvas.width;
var height = canvas.height;
var scene = new BABYLON.Scene(engine);
scene.clearColor = new BABYLON.Color3(1.0, 1.0, 1.0);
var camera = new BABYLON.ArcRotateCamera("camera1",
3 * Math.PI / 2, Math.PI / 2, 25, new BABYLON.Vector3(0, 8, 0), scene);
camera.attachControl(canvas, true);
var light1 = new BABYLON.HemisphericLight("light1", new BABYLON.Vector3(0, 1, 0), scene);
light1.intensity = 0.8;
var light2 = new BABYLON.HemisphericLight("light2", new BABYLON.Vector3(0, -1, 0), scene);
light2.intensity = 0.5;
require(['./WorldConstraint', './math'], function(World, math) {
let walls = [];
let wall_angle1 = -0.3;
let wall_angle2 = 0.5;
walls.push(new Wall(new math.Vector2( 0, 0),
new math.Vector2(-Math.sin(wall_angle1), Math.cos(wall_angle1))));
walls.push(new Wall(new math.Vector2( 0, 0),
new math.Vector2(-Math.sin(wall_angle2), Math.cos(wall_angle2))));
let bodies = [];
// body 0
bodies.push(
new Body(
createPolygonPoints(new math.Vector2(0, 0), 2, 4),
new math.Vector2(7, 14), // world position(pos)
Math.PI / 4, // angle
1 // mass
)
);
bodies[0].vel = new math.Vector2(0, 0);
bodies[0].angVel = 0.0;
// body 1
bodies.push(
new Body(
createPolygonPoints(new math.Vector2(0, 0), 1, 4),
new math.Vector2(-12, 14),
0,
1
)
);
bodies[1].vel = new math.Vector2(0, 0);
bodies[1].angVel = 0.0;
// body 2
bodies.push(
new Body(
createPolygonPoints(new math.Vector2(0, 0), 1.5, 4),
new math.Vector2(-8, 14),
0,
1
)
);
bodies[2].vel = new math.Vector2(2, 0);
bodies[2].angVel = 0.0;
//
console.log("bodies", Array.from(bodies));
var d, p1, p2;
for (let i = 0; i < walls.length; i++) {
d = new math.Vector2(400 * walls[i].normal.y, -400 * walls[i].normal.x);
p1 = walls[i].position.addV(d);
p2 = walls[i].position.subV(d);
let wallpath = [
new BABYLON.Vector3(p1.x, p1.y, 0),
new BABYLON.Vector3(p2.x, p2.y, 0)
];
let wallLines = BABYLON.Mesh.CreateLines("lines" + i, wallpath, scene);
wallLines.color = BABYLON.Color3.Black();
}
var i, j;
var objects, object;
var colors = [
BABYLON.Color3.Blue(),
BABYLON.Color3.Green(),
BABYLON.Color3.Black(),
];
var objectMesh = [];
objects = bodies;
for (j = 0; j < objects.length; ++j) {
object = objects[j];
var worldVertices = object.getWorldShape();
console.log("object index", j, worldVertices);
const abodyMat = new BABYLON.StandardMaterial("abodyMat" + j, scene);
abodyMat.diffuseColor = colors[j];
var points = [];
for (i = 0; i < worldVertices.length; ++i) {
const apoint = BABYLON.MeshBuilder.CreateSphere("apoint" + i,
{diameter: 0.3, segments: 8}, scene);
apoint.position.x = worldVertices[i].x;
apoint.position.y = worldVertices[i].y;
apoint.position.z = 0;
apoint.material = abodyMat;
points.push(apoint);
}
var path = [];
for (i = 0; i < worldVertices.length; i++)
path.push(
new BABYLON.Vector3(worldVertices[i].x, worldVertices[i].y, 0)
);
path.push(new BABYLON.Vector3(worldVertices[0].x, worldVertices[0].y, 0));
const options = {points: path, updatable: true};
var aLine = BABYLON.MeshBuilder.CreateLines("aLine" + j, options, scene);
aLine.color = colors[j];
objectMesh.push({points: points, line: aLine});
}
function update_mesh() {
var i, j
for (j = 0; j < objects.length; j++) {
var object = objects[j];
var worldVertices = object.getWorldShape();
//// points
var apoint = objectMesh[j].points;
for (i = 0; i < worldVertices.length; ++i) {
apoint[i].position.x = worldVertices[i].x;
apoint[i].position.y = worldVertices[i].y;
apoint[i].position.z = 0.0;
}
//// line
var aline = objectMesh[j].line;
var path = [];
for (i = 0; i < worldVertices.length; i++)
path.push(
new BABYLON.Vector3(worldVertices[i].x, worldVertices[i].y, 0)
);
path.push(new BABYLON.Vector3(worldVertices[0].x, worldVertices[0].y, 0));
var options = {points: path, instance: aline};
aline = BABYLON.MeshBuilder.CreateLines("aline", options);
}
}
let world = new World();
world.walls = walls;
world.bodies = bodies;
scene.registerBeforeRender(function () {
world.step();
update_mesh();
});
function createPolygonPoints(center = new math.Vector2(0, 0), radius = 1, nsides=3) {
let poly = [];
let delta = Math.PI * 2 / nsides;
for (let i = 0; i < nsides; i++) {
poly.push(
new math.Vector2(center.x+radius*Math.cos(delta*i),
center.y+radius*Math.sin(delta*i))
)
}
// inertia
let [a, b] = [(1 + Math.cos(delta)) / 2, (1 - Math.cos(delta)) / 2];
let inertia = radius * radius / 6 * ( 3*(a**3 + b**3) + 9*a*b - 2*b );
return {shape:poly, inertia:inertia};
}
function Wall(position, normal) {
//// position, normal: Vector2
this.position = position;
this.normal = normal;
this.dist = function(position) {
return this.normal.dot(position.subV(this.position));
};
}
function Body(shape, pos = new Vector2(0, 0), ang = 0, mass = 1) {
console.log("ceate Bodies", shape);
this.mass = mass;
this.invMass = 1 / this.mass;
this.inertia = this.mass * shape.inertia;
this.invInertia = 1 / this.inertia;
this.pos = pos;
this.ang = ang;
this.vel = new math.Vector2(0, 0);
this.angVel = 0;
let center = centerOfMass(shape.shape);
this.shape = shape.shape.map(p => new math.Vector2(p.x - center.x, p.y - center.y));
//// Returns the shape points in world space
this.getWorldShape = function () {
let [c, s] = [Math.cos(this.ang), Math.sin(this.ang)];
let transf = new math.Matrix3();
transf = transf.fromValues(c, -s, this.pos.x, s, c, this.pos.y, 0, 0, 0);
this.curWorldShape = transf.mulVectors(this.shape, 1);
return this.curWorldShape;
}
}
function centerOfMass(pts) {
//// pts: Array, element Vector2
let sum = pts.reduce((a,b) => new math.Vector2(a.x+b.x, a.y+b.y));
return new math.Vector2(sum.x/pts.length, sum.y/pts.length);
}
});
engine.runRenderLoop(() => {
scene.render();
});
}
init();
</script>
WorldConstraint.js
define(['./math', './gjk_epa'], function(math, Gjk) {
function centerOfMass(pts) {
let sum = pts.reduce((a,b) => new math.Vector2(a.x+b.x, a.y+b.y));
return new math.Vector2(sum.x/pts.length, sum.y/pts.length);
}
function supportEdge(pts, v) {
let maxdot = Number.NEGATIVE_INFINITY;
let maxdot2 = Number.NEGATIVE_INFINITY;
let best, best2;
for (let p of pts) {
let dot = p.dot(v);
if (dot > maxdot) {
[maxdot2, best2] = [maxdot, best];
[maxdot, best] = [dot, p];
} else if (dot > maxdot2) {
[maxdot2, best2] = [dot, p];
}
}
if (Math.abs(maxdot - maxdot2) < 0.01) return [best2, best];
return [best];
}
function closestSegmentPoint(p, q, r) {
let qr = r.subV(q);
let s = qr.getLength();
if (s < 0.00001) return q;
//// Degenerate line segment
//// both endpoints at approximately the same distance
let v = qr.normalized();
let u = p.subV(q);
let d = u.dot(v);
if (d < 0) return q;
if (d > s) return r;
return lerp(q, r, d/s);
}
function sap(polygons, v) {
let n = polygons.length;
let pairs = [];
let proj = [];
polygons.forEach((poly, i) => {
let min = Number.POSITIVE_INFINITY;
let max = Number.NEGATIVE_INFINITY;
for (let p of poly) {
let dot = p.dot(v);
min = Math.min(min, dot);
max = Math.max(max, dot);
}
proj.push([min, i], [max, i]);
});
proj.sort((a, b) => a[0] - b[0]);
let inside = new Set();
for (let [_, i] of proj) {
if (inside.has(i)) {
inside.delete(i);
} else {
pairs[i] = [];
for (let j of inside) {
if (i < j) pairs[j].push(i);
else pairs[i].push(j);
}
inside.add(i);
}
}
return pairs;
}
function lerp(a, b, t) {
//// glmatrix vec2 lerp: linear interpolation between two vec2's
//// a, b: Vector2, t: float (0 - 1)
var ax = a.x,
ay = a.y;
out_x = ax + t * (b.x - ax);
out_y = ay + t * (b.y - ay);
return new math.Vector2(out_x, out_y);
}
const g = 9.8;
const dt = 0.02;
class World {
constructor() {
this.nIterations = 4;
this.bodies = [];
this.walls = [];
}
step() {
//// Apply external force
this.bodies.forEach(body => {
body.vel = body.vel.addV(new math.Vector2(0, -g * dt));
});
//// Solve Constraints
//// Constraint Walls
this.projectWalls();
//// Constraint Collisions
this.projectCollision();
//// Update position and angle
this.bodies.forEach(body => {
body.pos = body.pos.addV(body.vel.mulS(dt));
body.ang += body.angVel * dt;
});
}
//// projectWalls (Resolve violated velocity constraint)
projectWalls() {
let normalImpulseSum = [];
let tangentImpulseSum = [];
let cnt2 = [];
for (let i = 0; i < this.bodies.length; i++) cnt2[i] = [];
let contact = [];
let contactIdx = 0;
this.bodies.forEach((body, j) => {
this.walls.forEach((wall, k) => {
let worldVertices = body.getWorldShape().slice();
let contactPoint = new math.Vector2(0, 0);
cnt2[j][k] = 0;
let dist_min = Number.POSITIVE_INFINITY;
for (let i = 0; i < worldVertices.length; i++) {
let dist = wall.dist(worldVertices[i]);
if (dist < dist_min) dist_min = dist;
if (dist <= 0) {
contactPoint = contactPoint.addV(worldVertices[i]);
cnt2[j][k]++;
}
}
if (cnt2[j][k] != 0) {
contactPoint = contactPoint.divS(cnt2[j][k]);
contact.push({contactPoint: contactPoint, dist_min: dist_min});
normalImpulseSum[contactIdx] = 0.0;
tangentImpulseSum[contactIdx] = 0.0;
contactIdx++;
}
});
});
for (let i = 0; i < this.nIterations; i++) {
let contactIdx = 0;
this.bodies.forEach((body, j) => {
this.walls.forEach((wall, k) => {
if (cnt2[j][k] != 0) {
applyImpulseWall(contactIdx, contact, body, wall);
contactIdx++;
}
});
});
}
function applyImpulseWall(cIdx, contact, body, wall) {
let contactNormal = wall.normal.inverted();
let ra = contact[cIdx].contactPoint.subV(body.pos);
//// Jacobian for non-penetration constraint
let j_va = contactNormal.inverted();
let j_wa = -ra.cross(contactNormal);
let beta = 0.4;
let restitution = 0.5;
//// Velocity at contact point
let vRot = new math.Vector2(
-body.angVel * ra.y,
body.angVel * ra.x);
let vp = body.vel.addV(vRot);
vp.invert();
let vpn = vp.dot(contactNormal);
let penetration_slop = 0.05;
let restitution_slop = 5.0;
let bias = -(beta / dt) * Math.max(-contact[cIdx].dist_min - penetration_slop, 0)
+ restitution * Math.max(vpn - restitution_slop, 0);
let k = body.invMass
+ j_wa * body.invInertia * j_wa;
let massNormal = 1.0 / k;
//// Jacobian * velocity vector
let jv = j_va.dot(body.vel) + j_wa * body.angVel;
let lambda = massNormal * -(jv + bias);
let previousTotalLambda = normalImpulseSum[cIdx];
normalImpulseSum[cIdx] = Math.max(0.0, normalImpulseSum[cIdx] + lambda);
lambda = normalImpulseSum[cIdx] - previousTotalLambda;
body.vel = body.vel.addV(j_va.mulS(body.invMass * lambda));
body.angVel += body.invInertia * j_wa * lambda;
//// Jacobian for friction constraint
let contactTangent = new math.Vector2(-contactNormal.y, contactNormal.x);
contactTangent.normalize();
j_va = contactTangent.inverted();
j_wa = -ra.cross(contactTangent);
jv = j_va.dot(body.vel) + j_wa * body.angVel;
k = body.invMass + j_wa * body.invInertia * j_wa;
let massTangent = 1.0 / k;
lambda = massTangent * -jv;
let friction = 0.01;
previousTotalLambda = tangentImpulseSum[cIdx];
let maxFriction = friction * normalImpulseSum[cIdx];
tangentImpulseSum[cIdx] = math.clamp(tangentImpulseSum[cIdx] + lambda, -maxFriction, maxFriction);
lambda = tangentImpulseSum[cIdx] - previousTotalLambda;
body.vel = body.vel.addV(j_va.mulS(body.invMass * lambda));
body.angVel = body.angVel + body.invInertia * j_wa * lambda;
} // function applyImpluseWall
} // projectWalls
//// projectCollision
projectCollision() {
let pairs = sap(this.bodies.map(body => body.getWorldShape()), new math.Vector2(1, 0));
let normalImpulseSum = [];
let tangentImpulseSum = [];
let nContacts = 0;
for (let j = 0; j < this.bodies.length; j++) {
nContacts += pairs[j].length;
}
for (let i = 0; i < nContacts; i++) {
normalImpulseSum[i] = 0.0;
tangentImpulseSum[i] = 0.0;
}
for (let i = 0; i < this.nIterations; i++) {
let contactIdx = 0;
this.bodies.forEach((body, j) => {
pairs[j].forEach(k => {
applyImpulse(contactIdx, body, this.bodies[k]);
contactIdx++;
});
});
} // for iteration
function applyImpulse(cIdx, bodyA, bodyB) {
let a = bodyA.getWorldShape().slice();
let b = bodyB.getWorldShape().slice();
let hit = Gjk.gjk(a, b);
if (hit) {
let { p, q, dist, n } = Gjk.epa(a, b, ...hit);
let aPoints = supportEdge(a, n);
let bPoints = supportEdge(b, n.mulS(-1));
let [massA, massB] = [bodyA.mass, bodyB.mass];
let aContact, bContact, aContactDisplaced, bContactDisplaced, penetration;
if (aPoints.length + bPoints.length == 4) {
//// Edge-edge collision
let center = centerOfMass([...aPoints, ...bPoints]);
aContact = closestSegmentPoint(center, ...aPoints);
bContact = closestSegmentPoint(center, ...bPoints);
penetration = aContact.subV(bContact).getLength();
aContactDisplaced =
aContact.addV(
n.mulS((-penetration * massA) / (massA + massB))
);
bContactDisplaced =
bContact.addV(
n.mulS((penetration * massB) / (massA + massB))
);
} else {
//// Vertex-edge collision
if (aPoints.length + bPoints.length != 3) {
throw "Weird collision";
}
if (aPoints.length == 2) {
aContact = closestSegmentPoint(bPoints[0], ...aPoints);
penetration = aContact.subV(bPoints[0]).getLength();
aContactDisplaced =
aContact.addV(
n.mulS((-penetration * massA) / (massA + massB))
);
bContactDisplaced =
bPoints[0].addV(
n.mulS((penetration * massB) / (massA + massB))
);
} else {
//// bPoints.length == 2!
bContact = closestSegmentPoint(aPoints[0], ...bPoints);
penetration = aPoints[0].subV(bContact).getLength();
bContactDisplaced =
bContact.addV(
n.mulS((penetration * massB) / (massA + massB))
);
aContactDisplaced =
aPoints[0].addV(
n.mulS((-penetration * massA) / (massA + massB))
);
}
}
let ra = aContactDisplaced.subV(bodyA.pos);
let rb = bContactDisplaced.subV(bodyB.pos);
//// Jacobian for non-penetration constraint
let j_va = n.inverted();
let j_wa = -ra.cross(n);
let j_vb = n;
let j_wb = rb.cross(n);
let beta = 0.2;
let restitution = 0.3;
//// Relative velocity at contact point
let vaRot = new math.Vector2(
-bodyA.angVel * ra.y, bodyA.angVel * ra.x);
let vbRot = new math.Vector2(
-bodyB.angVel * rb.y, bodyB.angVel * rb.x);
let vpa = bodyA.vel.addV(vaRot);
let vpb = bodyB.vel.addV(vbRot);
let vp = vpb.subV(vpa);
let vpn = vp.dot(n);
let penetration_slop = 0.1;
let restitution_slop = 4.0;
let bias = -(beta / dt) * Math.max(penetration - penetration_slop, 0) +
restitution * Math.max(vpn - restitution_slop, 0);
let k = +bodyA.invMass
+ j_wa * bodyA.invInertia * j_wa
+ bodyB.invMass
+ j_wb * bodyB.invInertia * j_wb;
let massNormal = 1.0 / k;
//// Jacobian * velocity vector
let jv = +j_va.dot(bodyA.vel)
+ j_wa * bodyA.angVel
+ j_vb.dot(bodyB.vel)
+ j_wb * bodyB.angVel;
let lambda = massNormal * -(jv + bias);
let previousTotalLambda = normalImpulseSum[cIdx];
normalImpulseSum[cIdx] = Math.max(0.0, normalImpulseSum[cIdx] + lambda);
lambda = normalImpulseSum[cIdx] - previousTotalLambda;
bodyA.vel = bodyA.vel.addV(j_va.mulS(bodyA.invMass * lambda));
bodyA.angVel += bodyA.invInertia * j_wa * lambda;
bodyB.vel = bodyB.vel.addV(j_vb.mulS(bodyB.invMass * lambda));
bodyB.angVel += bodyB.invInertia * j_wb * lambda;
//// Jacobian for friction constraint
let vTangent = new math.Vector2(-n.y, n.x);
vTangent.normalize();
j_va = vTangent.inverted();
j_wa = -ra.cross(vTangent);
j_vb = vTangent;
j_wb = rb.cross(vTangent);
jv = j_va.dot(bodyA.vel)
+ j_wa * bodyA.angVel;
+ j_vb.dot(bodyB.vel)
+ j_wb * bodyB.angVel;
k = bodyA.invMass
+ j_wa * bodyA.invInertia * j_wa;
+ bodyB.invMass
+ j_wb * bodyB.invInertia * j_wb;
let massTangent = 1.0 / k;
lambda = massTangent * -jv;
let friction = 0.2;
previousTotalLambda = tangentImpulseSum[cIdx];
let maxFriction = friction * normalImpulseSum[cIdx];
tangentImpulseSum[cIdx] = math.clamp(tangentImpulseSum[cIdx] + lambda, -maxFriction, maxFriction);
lambda = tangentImpulseSum[cIdx] - previousTotalLambda;
bodyA.vel = bodyA.vel.addV(j_va.mulS(bodyA.invMass * lambda));
bodyA.angVel = bodyA.angVel + bodyA.invInertia * j_wa * lambda;
bodyB.vel = bodyB.vel.addV(j_vb.mulS(bodyB.invMass * lambda));
bodyB.angVel = bodyB.angVel + bodyB.invInertia * j_wb * lambda;
} // end-if(hit)
} // end applyImplus function
} // end projectCollision()
}
return World;
});
物理シミュレーション 剛体の衝突 2 ( Force Based )
Force Based による、剛体の衝突シミュレーションです。
[ 実行結果 ]
２Dシミュレーション
メッシュの描画には、 BabylonJS を使用しています。
シミュレーションの方法
壁との衝突、剛体同士の衝突での力、トルクの計算は、以下を参照しています。
「WebGLによる 物理シミュレーション」( 酒井幸市、工学社 ）
第４章 剛体の衝突
１ 時間積分
剛体に作用する力
・重力 force_g
・壁との衝突時
壁から受ける力 force_wall
その力によるトルク（モーメント）torque_wall
・剛体同士の衝突時
他の剛体から受ける力 force_collision
その力によるトルク（モーメント）torque_collision
これらの力の総和(force)とトルクの総和(torque)から、加速度(acc)、
角加速度(angAcc)、速度(vel)、角速度(angVel)、そして位置(pos)、角度(ang)と
数値積分していきます。(html simulate() 関数）
２Dを考えているので、トルクはz成分のみとなる。
２ 壁との衝突
（BodyForce.js class Body, projectWalls メソッド）
衝突点に作用する力とそれによって生じるトルクを計算する。
まず、衝突点(図の赤丸)を調べる。四角や三角など多角形の剛体では、衝突の
仕方は２通り起こる。edge または vertex での衝突である。
これを判定するため、多角形の各頂点と壁との距離を計算し、距離がゼロ以下
となる頂点数を数える。１ならvertex、2 なら edge で衝突している。これで
衝突点を知ることができる。
衝突のような撃力が働く場合は、力積から力を求める。力積は、力とその力が
作用する時間との積である。また、力積は、力が働く前後の運動量の変化として
表すこともできる。
(１) 摩擦を考慮しない場合
この場合は、剛体は壁の法線方向の力を受ける。
力積を 、力を
、作用する時間を
、剛体の質量を
、衝突前後の速度を
、
として、
衝突点での力 によって生じるトルクは、
,
は剛体の重心から衝突点へのベクトル（重心-衝突点ベクトル）
(2Dを考えているので、トルクはz成分のみであるが、ベクトルで表している。)
角力積 は、角運動量の変化に等しいので、慣性モーメントを
、衝突前後の
角速度を 、
として、
2Dを考えているので、慣性モーメントはスカラー量である。
剛体が四角の場合、 (
四角の
方向と
方向の辺の
長さ)。（プログラムでは、多角形の慣性モーメントを用いている。）
剛体の衝突点の速度は、並進の速度と回転による速度の合成速度で、
衝突前：
衝突後：
力積と角力積の式、反発係数(restitution) の式
(
は衝突面の法線方向の単位ベクトル)
から、
力積 は法線方向の成分のみなので、以下のようになる。
これから、力とトルクが以下のように求まる。
(２) 摩擦を考慮した場合
この場合は、壁の接線方向の力(摩擦力)も働く。
衝突面の接線方向の単位ベクトルは、
動摩擦係数を 、法線方向の力を
((１)で求めた力)として、
摩擦力 は、
摩擦力によるトルク は、
摩擦力による力積は、
摩擦力による角力積は、
である。
衝突点の、衝突後の合成速度の接線成分 は、
合成速度 は、以下のとき、ゼロとなる。
この時、剛体は、滑らずに転がるようになる。
衝突している時間 内に
がゼロとなる動摩擦係数（臨界動摩擦係数）は、
・ の場合
となるので、衝突している時間
内に
がゼロとなる。
力積と角力積は、
とする。
・ の場合
なので、
とする。
これから、求める力とトルクは以下のようになる。
の場合
の場合
３ 剛体同士の衝突
（BodyForce.js class Body, projectCollision メソッド）
剛体同士の衝突でも、着目剛体が他の剛体から受ける力とトルクは、壁との衝突の場合と同様に計算できる。
衝突判定には、sweep-and-prune(sap) 法（html sap() 関数）と gjk-epa 法 (前回記事、gjk-epa.js)を用いている。
衝突は、edge 同士 または edge と vertex で起こる。sapで衝突ペアを検出し、gjk-epa から、edge で衝突しているのか、vertex で衝突しているのか、penetration の深さ、衝突の方向、衝突点の座標等の情報が得られる。
各衝突ペアについて、それぞれの剛体に作用する力積と角力積を計算し、それから力とトルクを求めます。
(１) 摩擦を考慮しない場合
剛体１の質量を 、衝突前後の速度を
、
、剛体２の質量を
、衝突前後の速度を
、
とする。また、剛体１に働く力を
、作用する時間を
とすると、
剛体１に作用する力積：
剛体２に作用する力積：
剛体１の慣性モーメントを 、衝突前後の角速度を
、
、また剛体２の慣性モーメントを
、衝突前後の角速度を
、
として、
剛体１に作用する角力積：
剛体２に作用する角力積：
、
、
、
は、それぞれ剛体１、２に作用するトルクと重心-衝突点ベクトルである。
衝突点での合成速度は、以下のようになる。
剛体１ 衝突前：
衝突後：
剛体２ 衝突前：
衝突後：
力積、角力積、反発係数(restitution) の式
(
は衝突面の法線方向の単位ベクトル)
から、剛体１に作用する力 とトルク
は、
剛体２に作用する力 とトルク
は、
(２) 摩擦を考慮した場合
壁との衝突の場合と同様に、臨界動摩擦係数を以下のように定義する。
は衝突面の接線方向の単位ベクトルである。
摩擦力による、剛体１、２の力とトルクは、以下のようになる。
の場合
剛体１
剛体２
の場合
剛体１
剛体２
プログラム
collision-force.html
<!DOCTYPE html>
<meta charset="utf-8">
<title>BJS Physical Animation ( force based )</title>
<script src="https://preview.babylonjs.com/babylon.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.5/require.min.js"></script>
<body>
<canvas id="canvas" width="600" height="430" style="border: 1px solid gray;"></canvas>
<script>
function init() {
const canvas = document.getElementById("canvas");
const engine = new BABYLON.Engine(canvas);
var width = canvas.width;
var height = canvas.height;
var scene = new BABYLON.Scene(engine);
scene.clearColor = new BABYLON.Color3(1.0, 1.0, 1.0);
var camera = new BABYLON.ArcRotateCamera("camera1",
3 * Math.PI / 2, Math.PI / 2, 25, new BABYLON.Vector3(0, 8, 0), scene);
camera.attachControl(canvas, true);
var light1 = new BABYLON.HemisphericLight("light1", new BABYLON.Vector3(0, 1, 0), scene);
light1.intensity = 0.8;
var light2 = new BABYLON.HemisphericLight("light2", new BABYLON.Vector3(0, -1, 0), scene);
light2.intensity = 0.5;
require(['./BodyForce', './math'], function(Body, math) {
let walls = [];
let wall_angle1 = -0.3;
let wall_angle2 = 0.5;
walls.push(new Wall(new math.Vector2( 0, 0),
new math.Vector2(-Math.sin(wall_angle1), Math.cos(wall_angle1))));
walls.push(new Wall(new math.Vector2( 0, 0),
new math.Vector2(-Math.sin(wall_angle2), Math.cos(wall_angle2))));
let bodies = [];
// body [0]
bodies.push(
new Body(
createPolygonPoints(new math.Vector2(0, 0), 1, 4),
new math.Vector2(-8, 10), // world position(pos)
Math.PI / 4, // angle
1 // mass
)
);
bodies[0].vel = new math.Vector2(2, 0);
bodies[0].angVel = 0.0;
// body [1]
bodies.push(
new Body(
createPolygonPoints(new math.Vector2(0, 0), 1, 3),
new math.Vector2(8, 10),
//Math.PI / 4, // for nsides = 4
0, // for nsides = 3
1
)
);
bodies[1].vel = new math.Vector2(0, 0);
bodies[1].angVel = 0.0;
//
console.log("bodies", Array.from(bodies));
var d, p1, p2;
for (let i = 0; i < walls.length; i++) {
d = new math.Vector2(400 * walls[i].normal.y, -400 * walls[i].normal.x);
p1 = walls[i].position.addV(d);
p2 = walls[i].position.subV(d);
let wallpath = [
new BABYLON.Vector3(p1.x, p1.y, 0),
new BABYLON.Vector3(p2.x, p2.y, 0)
];
let wallLines = BABYLON.Mesh.CreateLines("lines" + i, wallpath, scene);
wallLines.color = BABYLON.Color3.Black();
}
var i, j;
var objects, object;
var colors = [
BABYLON.Color3.Blue(),
BABYLON.Color3.Green(),
BABYLON.Color3.Black(),
];
var objectMesh = [];
objects = bodies;
for (j = 0; j < objects.length; ++j) {
object = objects[j];
console.log("object index", j, object.worldShape);
const abodyMat = new BABYLON.StandardMaterial("abodyMat" + j, scene);
abodyMat.diffuseColor = colors[j];
var points = [];
for (i = 0; i < object.worldShape.length; ++i) {
const apoint = BABYLON.MeshBuilder.CreateSphere("apoint" + i,
{diameter: 0.3, segments: 8}, scene);
apoint.position.x = object.worldShape[i].x;
apoint.position.y = object.worldShape[i].y;
apoint.position.z = 0;
apoint.material = abodyMat;
points.push(apoint);
}
var path = [];
for (i = 0; i < object.worldShape.length; i++)
path.push(
new BABYLON.Vector3(object.worldShape[i].x, object.worldShape[i].y, 0)
);
path.push(new BABYLON.Vector3(object.worldShape[0].x, object.worldShape[0].y, 0));
const options = {points: path, updatable: true};
var aLine = BABYLON.MeshBuilder.CreateLines("aLine" + j, options, scene);
aLine.color = colors[j];
objectMesh.push({points: points, line: aLine});
}
function update_mesh() {
var i, j
for (j = 0; j < objects.length; j++) {
var object = objects[j];
var worldVertices = object.worldShape;
//// points
var apoint = objectMesh[j].points;
for (i = 0; i < worldVertices.length; ++i) {
apoint[i].position.x = worldVertices[i].x;
apoint[i].position.y = worldVertices[i].y;
apoint[i].position.z = 0.0;
}
//// line
var aline = objectMesh[j].line;
var path = [];
for (i = 0; i < worldVertices.length; i++)
path.push(
new BABYLON.Vector3(worldVertices[i].x, worldVertices[i].y, 0)
);
path.push(new BABYLON.Vector3(worldVertices[0].x, worldVertices[0].y, 0));
var options = {points: path, instance: aline};
aline = BABYLON.MeshBuilder.CreateLines("aline", options);
}
}
scene.registerBeforeRender(function () {
simulate();
update_mesh();
});
function simulate() {
var g = 9.8;
var dt = 0.02;
//// Integrate
bodies.forEach(body => {
body.force = body.force0.addV(new math.Vector2(0, -body.mass * g));
body.torque = 0.0;
body.projectWalls(walls);
});
let pairs = sap(bodies.map(body => body.worldShape), new math.Vector2(1, 0));
bodies.forEach((body, i) => {
pairs[i].forEach(j => body.projectCollision(bodies[j]));
});
bodies.forEach(body => {
var acc = body.force.divS(body.mass);
body.vel = body.vel.addV(acc.mulS(dt));
body.pos = body.pos.addV(body.vel.mulS(dt));
var angAcc = body.invInertia * body.torque;
body.angVel += angAcc * dt;
body.ang += body.angVel * dt;
});
}
function createPolygonPoints(center = new math.Vector2(0, 0), radius = 1, nsides=3) {
let poly = [];
let delta = Math.PI * 2 / nsides;
for (let i = 0; i < nsides; i++) {
poly.push(
new math.Vector2(center.x+radius*Math.cos(delta*i),
center.y+radius*Math.sin(delta*i))
)
}
// inertia
let [a, b] = [(1 + Math.cos(delta)) / 2, (1 - Math.cos(delta)) / 2];
let inertia = radius * radius / 6 * ( 3*(a**3 + b**3) + 9*a*b - 2*b );
return {shape:poly, inertia:inertia};
}
function sap(polygons, v) {
let n = polygons.length;
let pairs = [];
let proj = [];
polygons.forEach((poly, i) => {
let min = Number.POSITIVE_INFINITY;
let max = Number.NEGATIVE_INFINITY;
for (let p of poly) {
let dot = p.dot(v);
min = Math.min(min, dot);
max = Math.max(max, dot);
}
proj.push([min, i], [max, i]);
});
proj.sort((a, b) => a[0] - b[0]);
let inside = new Set();
for (let [_, i] of proj) {
if (inside.has(i)) {
inside.delete(i);
} else {
pairs[i] = [];
for (let j of inside) {
if (i < j) pairs[j].push(i);
else pairs[i].push(j);
}
inside.add(i);
}
}
return pairs;
}
function Wall(position, normal) {
//// position, normal: Vector2
this.position = position;
this.normal = normal;
this.dist = function(position) {
return this.normal.dot(position.subV(this.position));
};
}
});
engine.runRenderLoop(() => {
scene.render();
});
}
init();
</script>
BodyForce.js
define(['./math', './gjk_epa'], function(math, Gjk) {
function centerOfMass(pts) {
//// pts: Array, element Vector2
let sum = pts.reduce((a,b) => new math.Vector2(a.x+b.x, a.y+b.y));
return new math.Vector2(sum.x/pts.length, sum.y/pts.length);
}
function supportEdge(pts, v) {
let maxdot = Number.NEGATIVE_INFINITY;
let maxdot2 = Number.NEGATIVE_INFINITY;
let best, best2;
for (let p of pts) {
let dot = p.dot(v);
if (dot > maxdot) {
[maxdot2, best2] = [maxdot, best];
[maxdot, best] = [dot, p];
} else if (dot > maxdot2) {
[maxdot2, best2] = [dot, p];
}
}
if (Math.abs(maxdot - maxdot2) < 0.01) return [best2, best];
return [best];
}
function closestSegmentPoint(p, q, r) {
//// p, q, r: Vector2
let qr = r.subV(q);
let s = qr.getLength();
if (s < 0.00001) return q;
//// Degenerate line segment
//// both endpoints at approximately the same distance
let v = qr.normalized();
let u = p.subV(q);
let d = u.dot(v);
if (d < 0) return q;
if (d > s) return r;
return lerp(q, r, d/s);
}
function lerp(a, b, t) {
//// glmatrix vec2 lerp: linear interpolation between two vec2's
//// a, b: Vector2, t: float (0 - 1)
var ax = a.x,
ay = a.y;
out_x = ax + t * (b.x - ax);
out_y = ay + t * (b.y - ay);
return new math.Vector2(out_x, out_y);
}
const g = 9.8;
const dt = 0.02;
const restitution = 0.5;
const restitution2 = 0.4;
class Body {
constructor(shape, pos = new Vector2(0, 0), ang = 0, mass = 1) {
console.log("ceate Bodies", shape);
this.mass = mass;
this.inertia = this.mass * shape.inertia;
this.invInertia = 1 / this.inertia;
this.pos = pos;
this.ang = ang;
this.vel = new math.Vector2(0, 0);
this.angVel = 0;
this.force = new math.Vector2(0, 0);
this.force0 = new math.Vector2(0, 0);
this.torque = 0; // z成分のみ
this.vGravityToPoint = new math.Vector2(0, 0);
let center = centerOfMass(shape.shape);
this.shape = shape.shape.map(p => new math.Vector2(p.x - center.x, p.y - center.y));
}
//// Returns the shape points in world space
get worldShape() {
let [c, s] = [Math.cos(this.ang), Math.sin(this.ang)];
let transf = new math.Matrix3();
transf = transf.fromValues(c, -s, this.pos.x, s, c, this.pos.y, 0, 0, 0);
this.curWorldShape = transf.mulVectors(this.shape, 1);
return this.curWorldShape;
}
//// projectWalls
projectWalls(walls) {
//// walls: Array, element: wall info
//// 注 wall.normal の方向
let muK = 0.1;
var worldVertices = this.worldShape.slice();
for (let wall of walls) {
var cnt = 0;
var vPointCollision = new math.Vector2(0, 0);
var dist_min = Number.POSITIVE_INFINITY;
for (let i = 0; i < worldVertices.length; i++) {
let dist = wall.dist(worldVertices[i]);
if (dist < dist_min) dist_min = dist;
if (dist <= 0) {
vPointCollision = vPointCollision.addV(worldVertices[i]);
cnt++;
}
}
var flagCollisionF = false;
if (cnt != 0) {
vPointCollision = vPointCollision.divS(cnt);
this.vGravityToPoint = vPointCollision.subV(this.pos);
if (dist_min <= 0 && this.vel.y < -0.3)
{
var c = 0.4;
var vRot = new math.Vector2(
-this.angVel * this.vGravityToPoint.y,
this.angVel * this.vGravityToPoint.x);
var vp = this.vel.addV(vRot);
var vpxn = vp.cross(wall.normal);
var vTangent = new math.Vector2(wall.normal.y * vpxn, -wall.normal.x * vpxn);
vTangent.normalize();
var rxn = this.vGravityToPoint.cross(wall.normal)
var a2 = this.invInertia * rxn;
var a3 = new math.Vector2(-a2 * this.vGravityToPoint.y,
a2 * this.vGravityToPoint.x);
var rikiseki = - (restitution + 1.0) * wall.normal.dot(vp) / (1.0/this.mass + wall.normal.dot(a3));
this.force = this.force.addV(wall.normal.mulS(rikiseki / dt));
var torq2 = rxn * c * rikiseki / dt;
this.torque += torq2;
////摩擦を考慮
var rxt = this.invInertia * this.vGravityToPoint.cross(vTangent);
var rxt2 = this.vGravityToPoint.cross(vTangent);
var c2 = vTangent.dot(
new math.Vector2(-rxt * this.vGravityToPoint.y,
rxt * this.vGravityToPoint.x));
var B = -vTangent.dot(vp) / (1.0/this.mass + c2);
var muC = Math.abs(B / rikiseki);
if(muK >= muC){
this.force = this.force.addV(vTangent.mulS(B/dt));
this.torque += rxt2 * c * B/dt;
}
else
{
this.force = this.force.addV(vTangent.mulS(-muK * rikiseki / dt));
this.torque += -rxt2 * c*muK * rikiseki / dt;
}
flagCollisionF = true;
}
if (!flagCollisionF) {
////重力と抗力によるトルク
var torq = this.vGravityToPoint.cross(new math.Vector2(0.0, this.mass * g)) * 0.5;
this.torque += torq;
}
}
if (dist_min < 0.0) this.pos.y -= dist_min;
}
}
//// projectCollision
projectCollision(other, gravityBias = 0) {
let a = this.worldShape.slice();
let b = other.worldShape.slice();
let hit = Gjk.gjk(a, b);
//---- 重心-衝突点ベクトル vGravityToPoint を求める
if (hit) {
let { p, q, dist, n } = Gjk.epa(a, b, ...hit);
let aPoints = supportEdge(a, n);
let bPoints = supportEdge(b, n.mulS(-1));
let [massA, massB] = [this.mass, other.mass];
let aContact, bContact, aContactDisplaced, bContactDisplaced, penetration;
if (aPoints.length + bPoints.length == 4) {
//// Edge-edge collision
let center = centerOfMass([...aPoints, ...bPoints]);
aContact = closestSegmentPoint(center, ...aPoints);
bContact = closestSegmentPoint(center, ...bPoints);
penetration = aContact.subV(bContact).getLength();
aContactDisplaced =
aContact.addV(
n.mulS((-penetration * massA) / (massA + massB))
);
bContactDisplaced =
bContact.addV(
n.mulS((penetration * massB) / (massA + massB))
);
} else {
//// Vertex-edge collision
if (aPoints.length + bPoints.length != 3) {
throw "Weird collision";
}
if (aPoints.length == 2) {
aContact = closestSegmentPoint(bPoints[0], ...aPoints);
penetration = aContact.subV(bPoints[0]).getLength();
aContactDisplaced =
aContact.addV(
n.mulS((-penetration * massA) / (massA + massB))
);
bContactDisplaced =
bPoints[0].addV(
n.mulS((penetration * massB) / (massA + massB))
);
} else {
//// bPoints.length == 2!
bContact = closestSegmentPoint(aPoints[0], ...bPoints);
penetration = aPoints[0].subV(bContact).getLength();
bContactDisplaced =
bContact.addV(
n.mulS((penetration * massB) / (massA + massB))
);
aContactDisplaced =
aPoints[0].addV(
n.mulS((-penetration * massA) / (massA + massB))
);
}
}
this.vGravityToPoint = aContactDisplaced.subV(this.pos);
other.vGravityToPoint = bContactDisplaced.subV(other.pos);
//---- 力積の計算
let muK = 0.5;
var rikiseki;
var vn1, vn2;
var muc,B, m1, m2, mI1, mI2;
var vVelocityP1, vVelocityP2;
var vTangential; //接線方向
var vA1 = new math.Vector2(0.0, 0.0);
var vA2 = new math.Vector2(0.0, 0.0);
var vRg1 = new math.Vector2(0.0, 0.0);
var vRg2 = new math.Vector2(0.0, 0.0);
this.force = new math.Vector2(0.0, 0.0);
this.torque = 0.0;
other.force = new math.Vector2(0.0, 0.0);
other.torque = 0.0;
vRg1 = this.vGravityToPoint.copy();
vRg2 = other.vGravityToPoint.copy();
m1 = this.mass;
m2 = other.mass;
mI1 = this.invInertia;
mI2 = other.invInertia;
var vRot1 = new math.Vector2(
-this.angVel * vRg1.y, this.angVel * vRg1.x);
vVelocityP1 = this.vel.addV(vRot1);
var vNormal = n;
vn1 = vVelocityP1.dot(vNormal);
var vRot2 = new math.Vector2(
-other.angVel * vRg2.y, other.angVel * vRg2.x);
vVelocityP2 = other.vel.addV(vRot2);
vn2 = vVelocityP2.dot(vNormal);
////衝突応答
var vp12 = vVelocityP1.subV(vVelocityP2);
var vp12xn = vp12.cross(vNormal);
vTangential = new math.Vector2(vNormal.y * vp12xn, -vNormal.x * vp12xn);
vTangential.normalize();
////力積
var a11 = vRg1.cross(vNormal);
var a12 = mI1 * a11;
var a13 = new math.Vector2(-a12 * vRg1.y, a12 * vRg1.x);
var a21 = vRg2.cross(vNormal);
var a22 = mI2 * a21;
var a23 = new math.Vector2(-a22 * vRg2.y, a22 * vRg2.x);
rikiseki = - (restitution2 + 1.0) * (vn1 - vn2) / (1.0/m1 + 1.0/m2
+ vNormal.dot(a13) + vNormal.dot(a23));
////力の総和
this.force = this.force.addV(vNormal.mulS(rikiseki / dt));
other.force = other.force.subV(vNormal.mulS(rikiseki / dt));
////トルクの総和
this.torque += a11 * rikiseki / dt;
other.torque -= a21 * rikiseki / dt;
////摩擦を考慮
var rxt1 = vRg1.cross(vTangential);
var rxt2 = vRg2.cross(vTangential);
vA1 = new math.Vector2(-mI1 * rxt1 * vRg1.y, mI1 * rxt1 * vRg1.x);
vA2 = new math.Vector2(-mI2 * rxt2 * vRg2.y, mI2 * rxt2 * vRg2.x);
B = -vTangential.dot(vp12) / (1.0/m1+1.0/m2+ vTangential.dot(vA1.addV(vA2)));
muc = Math.abs(B / rikiseki);
if(muK >= muc)
{
this.force = this.force.addV(vTangential.mulS(B/dt));
this.torque += rxt1 * B/dt;
other.force = other.force.subV(vTangential.mulS(B/dt));
other.torque -= rxt2 * B/dt;
} else
{
this.force = this.force.addV(vTangential.mulS(muK * rikiseki / dt));
this.torque += rxt1 * muK * rikiseki / dt;
other.force = other.force.subV(vTangential.mulS(muK * rikiseki / dt));
other.torque -= rxt2 * muK * rikiseki / dt;
}
////衝突時にめり込んだ分引き離す
////0.01は分離を確実にするため
this.pos = this.pos.subV(vNormal.mulS(penetration/4.0+0.01));
other.pos = other.pos.addV(vNormal.mulS(penetration/4.0+0.01));
}
}
}
return Body;
});
