examples/ibds.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>Impulse-Based Dynamic Simulation</title>
</head>
<body>
<main>
	<h1 style="text-align:center">Impulse-Based Dynamic Simulation</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.005" 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="2" step="1"></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>IBDS algorithm:</h2>
				This example shows the impulse-based dynamic simulation (IBDS) method [BFS05, BS06]. 
				IBDS is a predictor-corrector method to perform a simulation step for a constrained system. 
				For each position constraint $C(\mathbf x) = 0$ a corresponding velocity constraint is defined as the time derivative $\dot C(\mathbf x) = 0$.
				<ol>
					<li>position constraint solver</li>
					<li>time integration</li>
					<li>velocity constraint</li>
				</ol>
				Position constraint solver:
				<ol>
					<li class="nostyle"><b>loop</b></li>
					<li style="margin-left:40px">time integration to predict particle positions</li>
					<li style="margin-left:40px">compute Lagrange multiplier</li>
					<li style="margin-left:40px">determine correction impulse</li>
					<li style="margin-left:40px">update velocities</li>
				</ol>
				Velocity constraint solver:
				<ol>
					<li class="nostyle"><b>loop</b></li>
					<li style="margin-left:40px">compute Lagrange multiplier</li>
					<li style="margin-left:40px">determine correction impulse</li>
					<li style="margin-left:40px">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 = 0,$$
				where $d$ is the rest length between particles $\mathbf{x}_{i_1}$ and  $\mathbf{x}_{i_2}$.
				The corresponding velocity constraints are defined as 
				$$\dot C_i(\mathbf{x}_{i_1}, \mathbf{x}_{i_2}) = \frac{\mathbf x_{i_1} -\mathbf x_{i_2}}{\| \mathbf x_{i_1} -\mathbf x_{i_2} \|}(\mathbf v_{i_1} -\mathbf v_{i_2}) = 0.$$
				
				Note that no inversion handling is implemented in this example. So the model might get tangled. 
			
			<h3>1. Position constraint solver loop</h3>
				Correction impulses are computed and applied in a loop using a fixed number of iterations or until the error is within a user-defined tolerance. 

			<h3>2. Time integration</h3>
				In the first step of the position constraint solver 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>3. Compute Lagrange multipliers</h3>
			
				<p>The general form to compute the Lagrange multipliers in IBDS is
				$$\mathbf J \mathbf M^{-1} \mathbf J^T \boldsymbol \lambda = -\frac{1}{\Delta t} \mathbf C(\mathbf x^*),$$
				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)}{\Delta t \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>
				
				<p>To compute the Lagrange multipliers, the constraint gradients of the distance constraints are required which are computed as: </p>				
				<h4>Constraint gradients:</h4>
					$$\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. Correction impulse:</h3>

				The correction impulse is determined as:
				$$\mathbf p = \mathbf J^T \lambda$$

			<h3>5. Velocity update:</h3>
				Finally, the velocity of a constrained particle is updated:
					$$\Delta \mathbf v = \frac{1}{m} \mathbf p$$
					
			<h3>Velocity constraint solver loop</h3>
				The velocity constraints $\dot C(\mathbf x) = 0$ are solved analogously to the position constraints.
				The Lagrange multipliers are determined by solving: 
				$$\mathbf J \mathbf M^{-1} \mathbf J^T \boldsymbol \lambda = -\mathbf {\dot C}(\mathbf x).$$ 
				However, note that no position prediction is required since the solve is performed after the time integration step. 

			<h3>Time integration</h3>
				The particles are advected by numerical time integration. In our case we use a symplectic Euler method:
				$$\begin{align*}
					\mathbf v(t + \Delta t) &= \mathbf v(t) + \frac{\Delta t}{m} \mathbf f^{\text{ext}}(t) \\
					\mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t + \Delta t),
				\end{align*}$$
				where $\mathbf f^{\text{ext}}$ are the external forces.
			
			<h3>References</h3>
			<ul>
				<li>[BS06] Jan Bender and Alfred Schmitt. Fast Dynamic Simulation of Multi-Body Systems Using Impulses.  In Proceedings of Virtual Realitiy, Interactions, and Physical Simulation, 2006. Eurographics Association.</li>
				<li>[BFS05] Jan Bender, Dieter Finkenzeller and Alfred Schmitt. An impulse-based dynamic simulation system for VR applications. In Proceedings of Virtual Concept 2005 </li>
			</ul>
			
		</td></tr>
	</table>	
	
</main>

<script id="simulation_code" type="text/javascript">
	class Particle 
	{
		constructor (x, y) 
		{
			// position
			this.x = x;
			this.y = y;
			// predicted position
			this.px = x;
			this.py = y;
			// velocity
			this.vx = 0;
			this.vy = 0;
		}
	}
	
	class Constraint 
	{
		constructor (index1, index2, restLength) 
		{
			this.index1 = index1;
			this.index2 = index2;
			this.length = restLength;
			this.gradC_x = 0.0;
			this.gradC_y = 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.005;
			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);
		}
		
		// predict position of a particle by numerical integration using a symplectic Euler
		predictPosition(index, dt)
		{
			let p = this.particles[index];
		
			// point is static
			if (!this.isActive(index))
			{
				p.px = p.x;
				p.py = p.y;
			}
			else 
			{
				// compute preview using symplectic Euler
				p.px = p.x + dt * p.vx;
				p.py = p.y + dt * (p.vy + dt * this.gravity);
			}
		}
		
		// predictor-corrector solver for all constraints C
		positionConstraintSolver()
		{	
			// 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];
					
					// predict positions
					this.predictPosition(index1, this.timeStepSize);
					this.predictPosition(index2, this.timeStepSize);
					
					// compute length of constraint
					let dx = p1.px - p2.px;
					let dy = p1.py - p2.py;
					let dl = Math.sqrt(dx*dx + dy*dy)
					
					// predict displacement: C = ||x1-x2|| - restLength
					let C = dl - this.constraints[i].length;

					// constraint gradient 
					let gradC_x = p1.px - p2.px;
					let gradC_y = p1.py - p2.py;
					let l = Math.sqrt(gradC_x*gradC_x + gradC_y*gradC_y);
					if (l > 0.001)
					{
						gradC_x = gradC_x/l;
						gradC_y = gradC_y/l;
						
						// K = J M^-1 J^T, where J is the Jacobian of the constraint
						// in our case J J^T = identity since J is a unit vector
						let K = 1.0 / this.mass;						
						if (this.isActive(index1) && this.isActive(index2)) 
							K = 2*K;
							
						// compute Lagrange multiplier
						let lambda = -C/this.timeStepSize * (1/K);						
						
						// determine position correction
						if (this.isActive(index1))
						{
							p1.vx += 1.0/this.mass*lambda * gradC_x;
							p1.vy += 1.0/this.mass*lambda * gradC_y;
						}
						if (this.isActive(index2))
						{
							p2.vx -= 1.0/this.mass*lambda * gradC_x;
							p2.vy -= 1.0/this.mass*lambda * gradC_y;
						}
					}
				}
			}
		}
		
		// solver for all velocity constraints dC/dt
		velocityConstraintSolver()
		{
			// 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 velocity difference
					let dvx = p1.vx - p2.vx;
					let dvy = p1.vy - p2.vy;
					
					// constraint gradient
					let gradC_x = p1.x - p2.x;
					let gradC_y = p1.y - p2.y;
					let l = Math.sqrt(gradC_x*gradC_x + gradC_y*gradC_y);
					if (l > 0.001)
					{
						gradC_x = gradC_x/l;
						gradC_y = gradC_y/l;

						// C = (x1-x2)/|x1-x2| * (v1-v2)
						let C = gradC_x * dvx + gradC_y * dvy;
						
						let K = 1.0 / this.mass;
						if (this.isActive(index1) && this.isActive(index2)) 
							K = 2*K;
						
						// compute Lagrange multiplier
						let lambda = -C * (1/K);	
												
						// determine position correction
						if (this.isActive(index1))
						{
							p1.vx += 1.0/this.mass*lambda * gradC_x;
							p1.vy += 1.0/this.mass*lambda * gradC_y;
						}
						if (this.isActive(index2))
						{
							p2.vx -= 1.0/this.mass*lambda * gradC_x;
							p2.vy -= 1.0/this.mass*lambda * gradC_y;
						}
					}
				}
			}
		}
		
		// simulation step
		simulationStep()
		{	
			let dt = this.timeStepSize;
			
			// solve position constraints
			this.positionConstraintSolver();
		
			// perform time integration using symplectic Euler
			for (let i = 0; i < this.particles.length; i++) 
			{
				let p = this.particles[i];
				if (this.isActive(i))
				{
					// integrate velocity considering gravitational acceleration 
					p.vy = p.vy + dt * this.gravity;

					p.x = p.x + dt * p.vx;
					p.y = p.y + dt * p.vy;
				}
				else 
				{
					p.vx = 0.0;
					p.vy = 0.0;
				}
			}		
					
			// solve velocity constraints
			this.velocityConstraintSolver();
			
			// 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);
			
			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>