diff --git a/src/03-simulation/Simulation.ts b/src/03-simulation/Simulation.ts index ca4346b52becdf59a16ceb319e12b723a03dfe7d..340ece80290ba376ac9781fc762ce579ac19669b 100644 --- a/src/03-simulation/Simulation.ts +++ b/src/03-simulation/Simulation.ts @@ -3,10 +3,50 @@ import { PipelineObserver } from "../core/Pipeline"; import { Particle } from './physics/particle'; import { SimulationData } from "./SimulationData"; + +function adaptive(x: Array<number>, h: number, eMax = 1e-10, n: number = 5): number { + const rk4 = SimulationData.INTEGRATORS[2]; + + // 1. calculate approximate + let xFull = rk4.step([...x], h, Simulation.derivation); + let xHalf = rk4.step([...x], h/2, Simulation.derivation); + xHalf = rk4.step(xHalf, h/2, Simulation.derivation); + + // 2. calculate error + let eList = []; + for (let i = 0; i < x.length; i++) { + eList.push(Math.abs(xFull[i] - xHalf[i])); + } + let e = 0; + for (let i = 0; i < eList.length; i++) { + e += Math.pow(eList[i], 2); + } + + if (e < 1e-8) { + return h; + } + + e = Math.sqrt(e); // total error + + // 3. e = k * h^n + // e.g. h = 0.001 mit n = 5 => e = k * 0.001^5 + // => k = e / h^n + // => k * h^n = eMax + // <=> h = pow(eMax / k, 1/n) + // => h = pow(eMax / (e / h^n), 1/n) + + let newH = Math.pow(eMax / (e / Math.pow(h, n)), 1 / n); + if (isNaN(newH) || !isFinite(newH)) { + newH = h; + } + return newH; +} + export class Simulation implements PipelineObserver { public static readonly G = new THREE.Vector3(0, -9.807, 0); public static current: Particle; + public static h = .001; constructor() { SimulationData.generateTextile(); @@ -14,20 +54,34 @@ export class Simulation implements PipelineObserver { update(_deltaTime: number): void { - const h = _deltaTime / 100; - for (let t = 0; t < _deltaTime; t += h) { + + if (SimulationData.paused) { + return; + } + + for (let t = 0; t < _deltaTime; t += Simulation.h) { SimulationData.pList?.forEach(pSub => { pSub.forEach(p => { if (p.usePhysics) { Simulation.current = p; let x = [p.position.x, p.position.y, p.position.z, p.velocity.x, p.velocity.y, p.velocity.z]; - x = SimulationData.INTEGRATORS[SimulationData.integratorType].step(x, h, Simulation.derivation); + + if (SimulationData.stepSizeMethod == 1) { + Simulation.h = adaptive(x, Simulation.h); + } + + x = SimulationData.INTEGRATORS[SimulationData.integrationMethod] + .step(x, Simulation.h, Simulation.derivation); + p.setState(x[0], x[1], x[2], x[3], x[4], x[5]); } }); }) } + + console.log(Simulation.h); + } diff --git a/src/03-simulation/SimulationData.ts b/src/03-simulation/SimulationData.ts index efa529423f11cab8bdd686dc7de7f450e6ac2d2c..03a81b9ce1315a81b392c325cc70d96e251dc70e 100644 --- a/src/03-simulation/SimulationData.ts +++ b/src/03-simulation/SimulationData.ts @@ -9,12 +9,15 @@ document.addEventListener('keydown', (event: KeyboardEvent) => { if (event.key === 'e') { SimulationData.switchParticlePhysics(); + } else if (event.key === 'w') { + SimulationData.paused = !SimulationData.paused; } }); class SimulationData implements PipelineGUI, PipelineRenderable { public static container: THREE.Group = new THREE.Group(); + public static paused: boolean = false; // cloth settings public static readonly MAX_TOTAL_MASS = 10; @@ -52,19 +55,15 @@ class SimulationData implements PipelineGUI, PipelineRenderable { public static sColorBend: number = 0x0000ff; // integrator - public static readonly STEP_SIZE = {Fixed:0, Adaptive:1}; public static readonly INTEGRATORS: Integrator[] = [ new EulerIntegrator(), new MidpointIntegrator(), - new RungeKuttaIntegrator + new RungeKuttaIntegrator() ] - public static integratorType: number = 2; - public static stepSize: number = 0; - - public static setEuler = () => { SimulationData.integratorType = 0; } - public static setMidpoint = () => { SimulationData.integratorType = 1; } - public static setRungeKutta = () => { SimulationData.integratorType = 2; } + public static integrationMethod: number = 2; + public static stepSizeMethod: number = 0; + public static delta: number = 1e-3; public static resetExternalForce(): void { SimulationData.externalForce.x = 0; @@ -347,6 +346,7 @@ class SimulationData implements PipelineGUI, PipelineRenderable { const general = gui.addFolder('General'); + general.add(SimulationData, 'paused').listen().name("Pause Simulation"); general.add(SimulationData, 'totalMass', SimulationData.MIN_TOTAL_MASS, SimulationData.MAX_TOTAL_MASS) .step(0.1).name('Total Mass').onChange(SimulationData.changedParticleSize); general.add(SimulationData, 'width', SimulationData.MIN_WIDTH, SimulationData.MAX_WIDTH) @@ -359,9 +359,9 @@ class SimulationData implements PipelineGUI, PipelineRenderable { const integration = gui.addFolder('Integrator'); - integration.add(SimulationData, 'setEuler').name('Euler').onChange(()=>{console.log(SimulationData.integratorType)}); - integration.add(SimulationData, 'setMidpoint').name('Midpoint').onChange(()=>{console.log(SimulationData.integratorType)}); - integration.add(SimulationData, 'setRungeKutta').name('Runge Kutta (4th order)').onChange(()=>{console.log(SimulationData.integratorType)}); + integration.add(SimulationData, 'integrationMethod', { 'Euler': 0, 'Midpoint': 1, 'Runge-Kutta': 2 }).name('Method'); + integration.add(SimulationData, 'stepSizeMethod', { 'Fixed': 0, 'Adaptive': 1 }).name('Step Size').onChange((e) => { SimulationData.stepSizeMethod = e; console.log(e);}); + integration.add(SimulationData, 'delta', 1e-10, 1e-6, 1e-10).name('Error Max'); integration.open(); const extForce = general.addFolder('External Force'); diff --git a/src/03-simulation/physics/Integrators.ts b/src/03-simulation/physics/Integrators.ts index 39da71ccb6ecd91bee3371acb85b322fdba2b1cd..f47ee3cda63caf5337b070b86aa24579cd34e37c 100644 --- a/src/03-simulation/physics/Integrators.ts +++ b/src/03-simulation/physics/Integrators.ts @@ -1,18 +1,6 @@ abstract class Integrator { abstract step(x: Array<number>, h: number, f: (x: Array<number>) => Array<number>): Array<number>; - - public static adaptiveStep(integrator: Integrator, x: Array<number>, h: number, - f: (x: Array<number>) => Array<number>, errorMax: number): number { - let xFull = integrator.step(x, h, f); - let xHalf = integrator.step(x, h / 2, f); - let error: Array<number> = []; - for (let i = 0; i < xFull.length; i++) { - error.push(Math.abs(xFull[i] - xHalf[i])); - } - let errorApproxMax = error.reduce((a, b) => Math.max(a, b)); - return h * Math.sqrt(errorMax / errorApproxMax); - } } class EulerIntegrator extends Integrator {