feat: adaptive fixed

src/03-simulation/DemoSimulation.ts

 import { PipelineData } from "../core/Pipeline";
 import { Simulation } from "./Simulation";
-import { SimulationData } from './SimulationData';
+import { SimData } from './SimData';
 
 class SimulationDemo extends PipelineData {
     public constructor() { super(); }
 
     data(): void {
-        const data = new SimulationData();
+        const data = new SimData();
         this.addGUIElement(data);
         this.addRenderable(data);
 

src/03-simulation/Simulation.ts

 import * as THREE from 'three';
 import { PipelineObserver } from "../core/Pipeline";
 import { Particle } from './physics/particle';
-import { SimulationData } from "./SimulationData";
+import { SimData } from "./SimData";
 
+export class Simulation implements PipelineObserver {
 
-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);
+    public static readonly G = new THREE.Vector3(0, -9.807, 0);
+    public static current: Particle;
 
-    // 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);
-    }
+    public static timeAccumulator = 0;
 
-    if (e < 1e-8) {
-        return h;
+    constructor() { 
+        SimData.generateTextile();
     }
 
-    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)
+    update(_deltaTime: number): void {
 
-    let newH = Math.pow(eMax / (e / Math.pow(h, n)), 1 / n);
-    if (isNaN(newH) || !isFinite(newH)) {
-        newH = h;
-    }
-    return newH;
+        if (SimData.paused) {
+            return;
         }
 
-export class Simulation implements PipelineObserver {
+        const integrator = SimData.INTEGRATORS[SimData.integrationMethod];
+        Simulation.timeAccumulator += _deltaTime;
 
-    public static readonly G = new THREE.Vector3(0, -9.807, 0);
-    public static current: Particle;
-    public static h = .001;
+        // 1. find smallest step siz
+        let h = SimData.stepSize;
+        if (SimData.stepSizeMethod == 1) {
+            let hMin = Number.MAX_VALUE;
+            SimData.pList?.forEach(pSub => {
+                pSub.forEach(p => {
+                    let x = [p.position.x, p.position.y, p.position.z,
+                        p.velocity.x, p.velocity.y, p.velocity.z];
 
-    constructor() { 
-        SimulationData.generateTextile();
+                    let hTmp = integrator.adaptive(x, h, Simulation.derivation, SimData.eOrder,
+                        SimData.eFactor * Math.pow(10, SimData.eMaxExponent),
+                        SimData.hFactor * Math.pow(10, SimData.hMaxExponent),
+                        SimData.hFactor * Math.pow(10, SimData.hMinExponent));
+                    if (hTmp < hMin) {
+                        hMin = hTmp;
                     }
+                });
+            })
 
+            h = hMin;
 
-    update(_deltaTime: number): void {
-
-        if (SimulationData.paused) {
-            return;
+            if (Simulation.timeAccumulator > 1) {
+                Simulation.timeAccumulator = 0;
+                console.log(h);
+            }
         }
     
-        for (let t = 0; t < _deltaTime; t += Simulation.h) {
-            SimulationData.pList?.forEach(pSub => {
+        // 2. apply simualtion
+        for (let t = 0, i = 0; t < _deltaTime && i < 10000; t += h, i++) {
+            SimData.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];
 
-                        if (SimulationData.stepSizeMethod == 1) {
-                            Simulation.h = adaptive(x, Simulation.h);
-                        }
-
-                        x = SimulationData.INTEGRATORS[SimulationData.integrationMethod]
-                            .step(x, Simulation.h, Simulation.derivation);
+                        x = integrator.step(x, h, Simulation.derivation);
 
                         p.setState(x[0], x[1], x[2], x[3], x[4], x[5]);
                     }
                 });
             })
         }
-
-        console.log(Simulation.h);
-
     }
 
-
     public static derivation(x: Array<number>): Array<number> {
 
         // 1. create new array to store derivatives
@@ -99,10 +82,13 @@ export class Simulation implements PipelineObserver {
 
         // 3.1 accumulate forces
         let totalForce = new THREE.Vector3(0, 0, 0);
-        totalForce.add(SimulationData.externalForce);
-        totalForce.divideScalar(Math.pow(SimulationData.resolution, 2));
-        totalForce.divideScalar(Math.sqrt(Math.pow(x[0], 2) + Math.pow(x[1], 2) + Math.pow(x[2], 2)));
-        // totalForce.multiplyScalar(Math.sin(2 * Math.PI * Math.random()) + .5);
+        totalForce.add(SimData.forceFieldStrength);
+        totalForce.divideScalar(Math.pow(SimData.resolution, 2));
+
+        const dist = SimData.forceFieldPosition.distanceTo(new THREE.Vector3(x[0], x[1], x[2]));
+        if (dist > 0) {
+            totalForce.divideScalar(dist);
+        }
 
         // 3.2 accumulate spring forces
         Simulation.current.references.forEach(spring => {
@@ -111,7 +97,7 @@ export class Simulation implements PipelineObserver {
         });
 
         // 3.3 damping force
-        const dampingForce = Simulation.current.velocity.clone().multiplyScalar(-SimulationData.damping);
+        const dampingForce = Simulation.current.velocity.clone().multiplyScalar(-SimData.damping);
         totalForce.add(dampingForce);
 
         // 3.3 acceleration = force / mass

src/03-simulation/physics/Integrators.ts

 
 abstract class Integrator {
     abstract step(x: Array<number>, h: number, f: (x: Array<number>) => Array<number>): Array<number>;
+
+    /**
+     * 
+     * @param x current state
+     * @param h current step size
+     * @param f how to deriviate the state
+     * @param eMax tolerance for error
+     * @param eOrder error class of the used method
+     * @param hMax maximum step size
+     * @param hMin minimum step size
+     * @returns the next step size between hMax and hMin
+     */
+    public adaptive(x: Array<number>, h: number, f: (x: Array<number>) => Array<number>,
+        eOrder: number, eMax: number, hMax: number = 0.05, hMin: number = 0.00005): number {
+        // 1. calculate approximate
+        let xFull = this.step([...x], h, f);
+        let xHalf = this.step([...x], h / 2, f);
+        xHalf = this.step(xHalf, h / 2, f);
+
+        // 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);
+        }
+
+        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, eOrder)), 1 / eOrder);
+        if (isNaN(newH) || !isFinite(newH)) {
+            newH = h;
+        }
+
+        let clampedH = Math.min(hMax, Math.max(newH, hMin));
+        return clampedH;
+    }
 }
 
 class EulerIntegrator extends Integrator {
-    step(x: Array<number>, h: number, f: (x: Array<number>) => Array<number>): Array<number>
-    {
+    step(x: Array<number>, h: number, f: (x: Array<number>) => Array<number>): Array<number> {
         let y = f(x);
         let x_new = x.map((_, i) => x[i] + h * y[i]);
         return x_new;
@@ -21,8 +66,7 @@ class MidpointIntegrator extends Integrator {
         this.euler = new EulerIntegrator();
     }
 
-    step(x: Array<number>, h: number, f: (x: Array<number>) => Array<number>): Array<number>
-    {
+    step(x: Array<number>, h: number, f: (x: Array<number>) => Array<number>): Array<number> {
         let euler_step = this.euler.step(x, h / 2, f);
         let y = f(euler_step);
         let x_new = x.map((_, i) => x[i] + h * y[i]);
@@ -31,8 +75,7 @@ class MidpointIntegrator extends Integrator {
 }
 
 class RungeKuttaIntegrator extends Integrator {
-    step(x: Array<number>, h: number, f: (x: Array<number>) => Array<number>): Array<number>
-    {
+    step(x: Array<number>, h: number, f: (x: Array<number>) => Array<number>): Array<number> {
         let k1 = f(x);
         let k2 = f(x.map((_, i) => x[i] + h / 2 * k1[i]));
         let k3 = f(x.map((_, i) => x[i] + h / 2 * k2[i]));

src/03-simulation/physics/Particle.ts

@@ -14,12 +14,13 @@ export interface ParticleParameters {
     references: Spring[];
 
     // visual parameters
-    radius: number;
     color: number;
     wireframe: boolean;
 }
 
 export class Particle implements PipelineRenderable {
+    private static counter = 0;
+
     public position: THREE.Vector3;
     public velocity: THREE.Vector3;
 
@@ -37,10 +38,12 @@ export class Particle implements PipelineRenderable {
         this.references = params.references;
 
         this.mesh = new THREE.Mesh(
-            new THREE.SphereGeometry(params.radius, 4, 2),
+            new THREE.SphereGeometry(),
             new THREE.MeshBasicMaterial({ color: params.color, wireframe: params.wireframe })
         );
         this.mesh.position.copy(this.position);
+        this.mesh.name = `Particle ${Particle.counter++}`;
+        this.mesh.geometry.name = this.mesh.name;
     }
 
     public setColor(color: number): void {
@@ -60,7 +63,8 @@ export class Particle implements PipelineRenderable {
 
     public setRadius(radius: number): void {
         this.mesh.geometry.dispose();
-        this.mesh.geometry = new THREE.SphereGeometry(radius, 4, 2);
+        this.mesh.geometry = new THREE.SphereGeometry(Math.abs(radius), 4, 2);
+        this.mesh.geometry.name = this.mesh.name;
     }
 
     public object(): THREE.Object3D<THREE.Event> {

src/03-simulation/physics/Springs.ts

@@ -13,6 +13,7 @@ export interface SpringParameters {
 }
 
 export abstract class Spring {
+    private static counter = 0;
 
     public positionA: THREE.Vector3;
     public positionB: THREE.Vector3;
@@ -25,10 +26,19 @@ export abstract class Spring {
         this.positionB = params.positionB;
         this.restingLength = this.positionA.distanceTo(this.positionB);
 
+        const bufferGeometry = new THREE.BufferGeometry();
+        bufferGeometry.setAttribute('position', new THREE.BufferAttribute(new Float32Array(6), 3));
+        bufferGeometry.attributes.position.setXYZ(0, this.positionA.x, this.positionA.y, this.positionA.z);
+        bufferGeometry.attributes.position.setXYZ(1, this.positionB.x, this.positionB.y, this.positionB.z);
+        bufferGeometry.attributes.position.needsUpdate = true;
+
         this.line = new THREE.Line(
-            new THREE.BufferGeometry().setFromPoints([this.positionA, this.positionB]),
+            bufferGeometry,
             new THREE.LineBasicMaterial({ color: params.color })
         );
+
+        this.line.name = `Spring ${Spring.counter++}`;
+        this.line.geometry.name = this.line.name;
     }
 
     abstract get springConstant(): number;
@@ -55,6 +65,7 @@ export abstract class Spring {
         this.line.geometry.attributes.position.setXYZ(0, this.positionA.x, this.positionA.y, this.positionA.z);
         this.line.geometry.attributes.position.setXYZ(1, this.positionB.x, this.positionB.y, this.positionB.z);
         this.line.geometry.attributes.position.needsUpdate = true;
+        this.line.geometry.name = this.line.name;
     }
 }