examples/xpbd.html

<!doctype html>
<html class="no-js" lang="en">
<head>
	<meta charset="utf-8">
	<style>
		body {font-family: Helvetica, sans-serif;}
		table {background-color:#CCDDEE;text-align:left}
	</style>
	<script type="text/x-mathjax-config">
	  MathJax.Hub.Config({
		extensions: ["tex2jax.js"],
		jax: ["input/TeX", "output/HTML-CSS"],
		tex2jax: {
		  inlineMath: [ ['$','$'], ["\\(","\\)"] ],
		  displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
		  processEscapes: true
		},
		"HTML-CSS": { fonts: ["TeX"] }
	  });
	</script>
	<script type="text/javascript" aync src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.4/MathJax.js"></script>
	<title>eXtended Position-Based Dynamics (XPBD)</title>
</head>
<body>
<main>
	<h1 style="text-align:center">eXtended Position-Based Dynamics (XPBD)</h1>
	<table style="align_center;border-radius: 20px;padding: 20px;margin:auto">
		<col width="1100"> 
		<col width="400"> 
		<tr>
			<td>
				<canvas id="simCanvas" width="1024" height="768" style="border:2px solid #000000;border-radius: 20px;background-color:#EEEEEE">Your browser does not support the HTML5 canvas tag.</canvas>
			</td>
			<td>
				<table>		
				<col width="180" style="padding-right:10px"> 
				<col width="100"> 
				<tr>
					<td><label>Current time</label></td>
					<td><span id="time">0.00</span> s</td>
				</tr>
				<tr>
					<td><label>Time per sim. step</label></td>
					<td><span id="timePerStep">0.00</span> ms</td>
				</tr>
				<tr>
					<td><label># particles</label></td>
					<td><span id="numParticles">0</span></td>
				</tr>
				<tr>
					<td><label># constraints</label></td>
					<td><span id="numConstraints">0</span></td>
				</tr>
				<tr>
					<td><label for="widthInput">Width</label></td>
					<td><input onchange="gui.restart()" id="widthInput" type="number" value="40" step="1"></td>
				</tr>
				<tr>
					<td><label for="heightInput">Height</label></td>
					<td><input onchange="gui.restart()" id="heightInput" type="number" value="30" step="1"></td>
				</tr>
				<tr>
					<td><label for="fixedParticlesInput"># fixed particles</label></td>
					<td><select onchange="gui.restart()" id="fixedParticlesInput">
					  <option value="1">1</option>
					  <option value="2" selected="selected">2</option>
					  <option value="4">4</option>
					</select></td>
				</tr>
				<tr>
					<td><label for="timeStepSizeInput">Time step size</label></td>
					<td><input onchange="gui.sim.timeStepSize=parseFloat(value)" id="timeStepSizeInput" type="number" value="0.01" step="0.001"></td>
				</tr>
				<tr>
					<td><label for="iterationsInput">Iterations</label></td>
					<td><input onchange="gui.sim.numIterations=parseInt(value)" id="iterationsInput" type="number" value="5" step="1"></td>
				</tr>
				<tr>
					<td><label for="stiffnessInput">Stiffness</label></td>
					<td><input onchange="gui.sim.stiffness=parseFloat(value)" id="stiffnessInput" type="number" value="1000.0" step="10.0"></td>
				</tr>
				<tr>
					<td><label for="gravityInput">Gravity</label></td>
					<td><input onchange="gui.sim.gravity=parseFloat(value)" id="gravityInput" type="number" value="-9.81" step="0.01"></td>
				</tr>
				<tr>
					<td><label for="massInput">Mass</label></td>
					<td><input onchange="gui.sim.mass=parseFloat(value)" id="massInput" type="number" value="0.5" step="0.01"></td>
				</tr>
				<tr>
					<td></td>
					<td><button onclick="gui.restart()" type="button" id="restart">Restart</button></td>
				</tr>
				<tr>
					<td></td>
					<td><button onclick="gui.doPause()" type="button" id="Pause">Pause</button></td>
				</tr>
				</table>
			</td>
		</tr>
		<tr><td>
			<h2>XPBD algorithm:</h2>
				This example shows the eXtended Position-Based Dynamics (XPBD) method introduced by Macklin et al. [MMC16,BMM17]. When simulating deformable solids using the original PBD method, we have the problem that the constraint stiffness depends on the time step size and the iteration count. XPBD solves this problem by using a compliance formulation. The method performs the following steps:
				<ol>
					<li>time integration to predict particle positions</li>
					<li class="nostyle"><b>loop</b></li>
					<li style="margin-left:40px">compute Lagrange multipliers</li>
					<li style="margin-left:40px">determine position correction</li>
					<li>update velocities</li>
				</ol>
				In this example we use distance constraints 
				$$C_i(\mathbf{x}_{i_1}, \mathbf{x}_{i_2}) = \| \mathbf x_{i_1} -\mathbf x_{i_2} \|-d,$$
				where $d$ is the rest length between particles $\mathbf{x}_{i_1}$ and  $\mathbf{x}_{i_2}$. Moreover, we define a stiffness $k$ for each distance constraint. So finally each constraint behaves like a spring with stiffness $k$.
			
			<h3>1. Time integration</h3>
				In the first step the particle positions are advected using a symplectic Euler method to obtain predicted positions $\mathbf x^*$:
				$$\begin{align*}
					\mathbf v^* &= \mathbf v(t) + \Delta t \mathbf a_\text{ext}(t) \\
					\mathbf x^* &= \mathbf x(t) + \Delta t \mathbf v^*,
				\end{align*}$$
				where $\mathbf a_\text{ext}$ are accelerations due to external forces.
				
			<h3>2. Solver loop</h3>
				Position corrections are computed and applied in a loop using a fixed number of iterations. 
				
			<h3>3. Compute Lagrange multipliers</h3>
			
				<h4>PBD</h4>
				<p>The general form to compute the Lagrange multipliers in position-based dynamics is
				$$\mathbf J \mathbf M^{-1} \mathbf J^T \boldsymbol \lambda = -\mathbf C(\mathbf p),$$
				where $\mathbf J = (\partial C / \partial \mathbf x)^T$ is the Jacobian of the constraint function and $\mathbf M$ denotes the mass matrix.</p>
				
				<p>Since the distance constraint is a scalar function the Lagrange multiplier is determined as
				$$\lambda_i = - \frac{C_i(\mathbf x)}{\sum_j \frac{1}{m_j} \| \partial C_i(\mathbf x) / \partial \mathbf x_j\|^2},$$
				where $j$ denotes the indices of the particles which are influenced by the constraint.
				</p>
				
				<h4>XPBD</h4>
				<p>When using XPBD, we do not directly compute a Lagrange multiplier in each iteration. Instead we start with $\lambda_i = 0$ and then in each iteration we compute a change of the Lagrange multiplier as
				</p>
				<p>$$\Delta \lambda_i = \frac{-C_i(\mathbf x) - \tilde \alpha_i \lambda_i}{\sum_j \frac{1}{m_j} \| \partial C_i(\mathbf x) / \partial \mathbf x_j\|^2 + \tilde \alpha_i},$$
				where $\tilde \alpha = \frac{\alpha}{\Delta t^2}$ and $\alpha$ is the compliance value which is in our example $\alpha = \frac 1 k$. </p>
				
				<p>Note, in this example we set $\alpha = 0$ if $k$ is set to 0 by the user. In this special case the algorithm is equivalent to PBD.</p>
						
				<h4>Constraint gradients:</h4>
				<p>To compute the Lagrange multipliers, the constraint gradients are required which are computed as: </p>		
					$$\begin{align*}
					 \frac{\partial C_i}{\partial \mathbf x_{i_1}} &= \frac{\mathbf x_{i_1} -\mathbf x_{i_2}}{\| \mathbf x_{i_1} -\mathbf x_{i_2} \|}  \\
					 \frac{\partial C_i}{\partial \mathbf x_{i_2}} &= - \frac{\mathbf x_{i_1} -\mathbf x_{i_2}}{\| \mathbf x_{i_1} -\mathbf x_{i_2} \|}.
					\end{align*}$$
					
			<h3>4. Position correction:</h3>
				The position correction for a particle is determined using the change of the Lagrange multiplier:
					$$\begin{align*}
					\Delta \mathbf x_{i_1} &= \frac{1}{m_{i_1}} \frac{\partial C_i}{\partial \mathbf x_{i_1}} \Delta \lambda_i  \\
					\Delta \mathbf x_{i_2} &= \frac{1}{m_{i_2}} \frac{\partial C_i}{\partial \mathbf x_{i_2}} \Delta \lambda_i 
					\end{align*}$$
				and the correction is applied to update $\mathbf x^*$ as
					$$\mathbf x^* := \mathbf x^* + \Delta \mathbf x.$$
					
			<h3>5. Velocity update:</h3>
				The final positions and velocities are computed as:
				$$\begin{align*}
					\mathbf x_i(t+\Delta t) &= \mathbf x_i^* \\
					\mathbf v_i(t+\Delta t) &= \frac{1}{\Delta t} (\mathbf x_i(t+\Delta t) - \mathbf x_i(t))
				\end{align*}
				$$	
			
			<h3>References</h3>
			<ul>
			<li>[MMC16] Miles Macklin, Matthias Müller, Nuttapong Chentanez. XPBD: Position-Based Simulation of Compliant Constrained Dynamics. In Proceedings of ACM Motion in Games, 2016</li>
			<li>[BMM17] Jan Bender, Matthias Müller, Miles Macklin. A Survey on Position Based Dynamics, 2017. Eurographics Tutorials, 2017</li>
			</ul>
			
		</td></tr>
	</table>	
	
</main>

<script id="simulation_code" type="text/javascript">
	class Particle 
	{
		constructor (x, y) 
		{
			this.x = x;
			this.y = y;
			this.xOld = x;
			this.yOld = y;
			this.vx = 0;
			this.vy = 0;
		}
	}
	
	class Constraint 
	{
		constructor (index1, index2, restLength) 
		{
			this.index1 = index1;
			this.index2 = index2;
			this.length = restLength;
			this.lambda = 0.0;
		}
	}
	
	class Simulation
	{
		constructor(width, height, fixed)
		{
			this.constraints = [];
			this.particles = [];
			this.width = width;
			this.height = height;
			this.numIterations = 2;
			this.mass = 0.5;
			this.gravity = -9.81;
			this.timeStepSize = 0.05;
			this.stiffness = 1000;
			this.time = 0;
			this.numFixedParticles = fixed;

			this.init();
		}
		
		init()
		{			
			// create particles  
			let i;
			let j;
			let w = this.width;
			let h = this.height;
			let s = 0.1;
			for (i = 0; i < h; i++) 
				for (j = 0; j < w; j++) 
					this.particles.push(new Particle(s*2*j, -s*2*i))
			
			// create horizontal constraints
			for (i = 0; i < h; i++)
				for (j = 0; j < w - 1; j++)
					this.constraints.push(new Constraint(i*w + j, i*w + j + 1, s*2.0));
			
			// create vertical constraints
			for (i = 0; i < h-1; i++)
				for (j = 0; j < w; j++)
					this.constraints.push(new Constraint(i*w + j, i*w + j + w, s*2.0));
			
			// create diagonal constraints
			for (i = 0; i < h-1; i++)
			{
				for (j = 0; j < w-1; j++)
				{
					this.constraints.push(new Constraint(i*w + j, i*w + j + w + 1, s*Math.sqrt(8)));			
					this.constraints.push(new Constraint(i*w + j +1, i*w + j + w, s*Math.sqrt(8)));
				}
			}
		}
		
		isActive(i)
		{
			// particle 0 and the selected particle are fixed
			if (this.numFixedParticles == 1)
				return (i != 0) && (i != gui.selectedParticle);
			else if (this.numFixedParticles == 2)
				return (i != 0) && (i != this.width-1) && (i != gui.selectedParticle);
			else 
				return (i != 0) && (i != this.width-1) && (i != (this.height-1)*this.width) && (i != this.height*this.width-1) && (i != gui.selectedParticle);
		}
		
		// simulation step
		simulationStep()
		{	
			let dt = this.timeStepSize;
			let alpha = 0.0;
			if (this.stiffness != 0.0)
				alpha = 1.0/this.stiffness;
			let alpha_tilde = alpha / (dt*dt);
		
			// compute preview using symplectic Euler
			for (let i = 0; i < this.particles.length; i++) 
			{
				if (this.isActive(i))
				{
					let p = this.particles[i];
					
					// store old position for velocity update after the position correction
					p.xOld = p.x;	
					p.yOld = p.y;
				
					// integrate velocity considering gravitational acceleration 
					p.vy = p.vy + dt * this.gravity
				
					// integrate position
					p.x = p.x + dt * p.vx;
					p.y = p.y + dt * p.vy;
				}
			}
			
			// reset lambda
			for (let i = 0; i < this.constraints.length; i++) 
				this.constraints[i].lambda = 0.0;
					
			// constraint handling
			for (let iter = 0; iter < this.numIterations; iter++) 
			{
				for (let i = 0; i < this.constraints.length; i++) 
				{
					let index1 = this.constraints[i].index1;
					let index2 = this.constraints[i].index2;
					let p1 = this.particles[index1];
					let p2 = this.particles[index2];
				
					// compute length of constraint
					let dx = p1.x - p2.x;
					let dy = p1.y - p2.y;
					let dl = Math.sqrt(dx*dx + dy*dy)
					
					// compute displacement: C = ||x1-x2|| - restLength
					let C = dl - this.constraints[i].length;
										
					// determine gradient
					let gradC_x = 0.0;
					let gradC_y = 0.0;
					if (dl > 0.001)
					{
						gradC_x = dx/dl; 
						gradC_y = dy/dl; 
						let K = (gradC_x*gradC_x + gradC_y*gradC_y) / this.mass;
						if (i != 0)
							K = 2*K;
						K += alpha_tilde;
						
						// comute Lagrange multiplier
						let delta_lambda = -(C + alpha_tilde*this.constraints[i].lambda)/ K;	
						this.constraints[i].lambda += delta_lambda;
						
						// determine position correction
						if (this.isActive(index1))
						{
							p1.x = p1.x + 1.0/this.mass*delta_lambda * gradC_x;
							p1.y = p1.y + 1.0/this.mass*delta_lambda * gradC_y;
						}
						if (this.isActive(index2))
						{
							p2.x = p2.x - 1.0/this.mass*delta_lambda * gradC_x;
							p2.y = p2.y - 1.0/this.mass*delta_lambda * gradC_y;
						}
					}
				}
			}
				
			// update velocities
			for (let i = 0; i < this.particles.length; i++) 
			{
				if (this.isActive(i))
				{
					let p = this.particles[i];
					p.vx = 1.0/dt * (p.x - p.xOld);
					p.vy = 1.0/dt * (p.y - p.yOld);
				}
			}
			
			// update simulation time
			this.time = this.time + dt;
		}
	}
	
	class GUI
	{
		constructor()
		{
			this.canvas = document.getElementById("simCanvas");
			this.c = this.canvas.getContext("2d");
			this.requestID = -1;
			
			this.timeSum = 0.0;
			this.counter = 0;
			this.pause = false;
			this.origin = { x : this.canvas.width / 2, y : this.canvas.height/2-200};
			this.zoom = 50;
			this.particleRadius = 0.025;
			this.selectedParticle = -1;
			
			// register mouse event listeners (zoom/selection)
			this.canvas.addEventListener("mousedown", this.mouseDown.bind(this), false);
			this.canvas.addEventListener("mousemove", this.mouseMove.bind(this), false);
			this.canvas.addEventListener("mouseup", this.mouseUp.bind(this), false);	
			this.canvas.addEventListener("wheel", this.wheel.bind(this), false);		
		}
		
		// set simulation parameters from GUI and start mainLoop
		restart()
		{
			window.cancelAnimationFrame(this.requestID);
			let w = parseInt(document.getElementById('widthInput').value);
			let h = parseInt(document.getElementById('heightInput').value);
			let f = parseInt(document.getElementById('fixedParticlesInput').value);
			delete this.sim;
			this.sim = new Simulation(w, h, f);
			
			this.timeSum = 0.0;
			this.counter = 0;

			this.sim.numIterations = parseInt(document.getElementById('iterationsInput').value);
			this.sim.gravity = parseFloat(document.getElementById('gravityInput').value);
			this.sim.timeStepSize = parseFloat(document.getElementById('timeStepSizeInput').value);
			this.sim.mass = parseFloat(document.getElementById('massInput').value);
			this.sim.stiffness = parseFloat(document.getElementById('stiffnessInput').value);
			
			document.getElementById("numParticles").innerHTML = this.sim.particles.length;
			document.getElementById("numConstraints").innerHTML = this.sim.constraints.length;
			this.mainLoop();
		}
		
		drawCoordinateSystem()
		{
			// draw x-axis
			this.c.strokeStyle = "#FF0000";
			this.c.beginPath();
			this.c.moveTo(this.origin.x, this.origin.y);
			this.c.lineTo(this.origin.x+1*this.zoom, this.origin.y);
			this.c.stroke();	

			// draw y-axis
			this.c.strokeStyle = "#00FF00";	
			this.c.beginPath();
			this.c.moveTo(this.origin.x, this.origin.y);
			this.c.lineTo(this.origin.x, this.origin.y-1*this.zoom);
			this.c.stroke();		
		}
		
		draw()
		{
			this.c.clearRect(0, 0, this.canvas.width, this.canvas.height);
			
			this.drawCoordinateSystem();
			
			// draw constraints as lines
			this.c.strokeStyle = "#666666";
			for (let i = 0; i < this.sim.constraints.length; i++) 
			{
				let index1 = this.sim.constraints[i].index1;
				let index2 = this.sim.constraints[i].index2;
				let p1 = this.sim.particles[index1];
				let p2 = this.sim.particles[index2];
								
				this.c.beginPath();
				this.c.moveTo(this.origin.x + p1.x*this.zoom, this.origin.y - p1.y*this.zoom);
				this.c.lineTo(this.origin.x + p2.x*this.zoom, this.origin.y - p2.y*this.zoom);
				this.c.stroke();	
			}
			
			// draw particles as circles
			for (let i = 0; i < this.sim.particles.length; i++) 
			{
				let p = this.sim.particles[i];
				let r = this.particleRadius;
				if (i == this.selectedParticle)
				{
					// draw selected particle in red with larger radius
					this.c.fillStyle = "#FF0000";	
					r = 3*r;
				}
				else 
					this.c.fillStyle = "#0000FF";	
				let px = this.origin.x + p.x * this.zoom;
				let py = this.origin.y - p.y * this.zoom;
				
				this.c.beginPath();			
				this.c.arc(px, py, r * this.zoom, 0, Math.PI*2, true); 
				this.c.closePath();
				this.c.fill();
			}
		}
		
		mainLoop()
		{
			// perform multiple sim steps per render step
			for (let i=0; i < 8; i++)
			{
				let t0 = performance.now();
				
				this.sim.simulationStep();

				let t1 = performance.now();
				this.timeSum += t1 - t0;
				this.counter += 1;
					
				if (this.counter % 50 == 0)				
				{
					this.timeSum /= this.counter;
					document.getElementById("timePerStep").innerHTML = this.timeSum.toFixed(2);		
					this.timeSum = 0.0;
					this.counter = 0;
				}	
				document.getElementById("time").innerHTML = this.sim.time.toFixed(2);	
			}			

			this.draw();
			if (!this.pause)
				this.requestID = window.requestAnimationFrame(this.mainLoop.bind(this));		
		}
		
		doPause()
		{
			this.pause = !this.pause;
			if (!this.pause)
				this.mainLoop();
		}
		
		mouseDown(event)
		{
			// left mouse button down
			if (event.which == 1)
			{
				let mousePos = this.getMousePos(this.canvas, event);
				for (let i = 0; i < this.sim.particles.length; i++) 
				{
					let p = this.sim.particles[i];
					let px = this.origin.x + p.x * this.zoom;
					let py = this.origin.y - p.y * this.zoom;
					let dx = px - mousePos.x
					let dy = py - mousePos.y
					let dist2 = Math.sqrt(dx * dx + dy * dy)
					if (dist2 < 10)
					{
						this.selectedParticle = i;
						break;
					}				
				}			
			}
		}
		
		getMousePos(canvas, event) 
		{
			let rect = canvas.getBoundingClientRect();
			return {
				x: event.clientX - rect.left,
				y: event.clientY - rect.top
			};
		}
	
		mouseMove(event)
		{
			if (this.selectedParticle != -1)
			{
				let mousePos = this.getMousePos(this.canvas, event);		
				this.sim.particles[this.selectedParticle].x = (mousePos.x - this.origin.x) / this.zoom;
				this.sim.particles[this.selectedParticle].y = -(mousePos.y - this.origin.y) / this.zoom;
			}
		}
		
		mouseUp(event) 
		{
			this.selectedParticle = -1;
		}
		
		wheel(event)
		{
			event.preventDefault();
			this.zoom += event.deltaY * -0.05;
			if (this.zoom < 1) 
				this.zoom = 1;
		}
	}
	
	gui = new GUI();
	gui.restart();
</script>

</body>
</html>