inochi2d.phys.system source code

1 /*
2     Copyright © 2022, Inochi2D Project
3     Distributed under the 2-Clause BSD License, see LICENSE file.
4 
5     Author Asahi Lina
6 */
7 module inochi2d.phys.system;
8 import inochi2d;
9 import std.math : isFinite;
10 
11 abstract class PhysicsSystem {
12 private:
13     ulong[float*] variableMap;
14     float*[] refs;
15     float[] derivative;
16 
17     float t;
18 
19 protected:
20     /**
21         Add a float variable to the simulation
22     */
23     ulong addVariable(float* var) {
24         ulong index = refs.length;
25 
26         variableMap[var] = index;
27         refs ~= var;
28 
29         return index;
30     }
31 
32     /**
33         Add a vec2 variable to the simulation
34     */
35     ulong addVariable(vec2* var) {
36         ulong index = addVariable(&(var.vector[0]));
37         addVariable(&(var.vector[1]));
38         return index;
39     }
40 
41     /**
42         Set the derivative of a variable (solver input) by index
43     */
44     void setD(ulong index, float value) {
45         derivative[index] = value;
46     }
47 
48     /**
49         Set the derivative of a float variable (solver input)
50     */
51     void setD(ref float var, float value) {
52         ulong index = variableMap[&var];
53         setD(variableMap[&var], value);
54     }
55 
56     /**
57         Set the derivative of a vec2 variable (solver input)
58     */
59     void setD(ref vec2 var, vec2 value) {
60         setD(var.vector[0], value.x);
61         setD(var.vector[1], value.y);
62     }
63 
64     float[] getState() {
65         float[] vals;
66 
67         foreach(idx, ptr; refs) {
68             vals ~= *ptr;
69         }
70 
71         return vals;
72     }
73 
74     void setState(float[] vals) {
75         foreach(idx, ptr; refs) {
76             *ptr = vals[idx];
77         }
78     }
79 
80     /**
81         Evaluate the simulation at a given time
82     */
83     abstract void eval(float t) {
84     }
85 
86 public:
87     /**
88         Run a simulation tick (Runge-Kutta method)
89     */
90     void tick(float h) {
91         float[] cur = getState();
92         float[] tmp;
93         tmp.length = cur.length;
94         derivative.length = cur.length;
95         foreach(i; 0..derivative.length)
96             derivative[i] = 0;
97 
98         eval(t);
99         float[] k1 = derivative.dup;
100 
101         foreach(i; 0..cur.length)
102             *refs[i] = cur[i] + h * k1[i] / 2f;
103         eval(t + h / 2f);
104         float[] k2 = derivative.dup;
105 
106         foreach(i; 0..cur.length)
107             *refs[i] = cur[i] + h * k2[i] / 2f;
108         eval(t + h / 2f);
109         float[] k3 = derivative.dup;
110 
111         foreach(i; 0..cur.length)
112             *refs[i] = cur[i] + h * k3[i];
113         eval(t + h);
114         float[] k4 = derivative.dup;
115 
116         foreach(i; 0..cur.length) {
117             *refs[i] = cur[i] + h * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6f;
118             if (!isFinite(*refs[i])) {
119                 // Simulation failed, revert
120                 foreach(j; 0..cur.length)
121                     *refs[j] = cur[j];
122                 break;
123             }
124         }
125 
126         t += h;
127     }
128 
129 public:
130 
131     abstract void drawDebug(mat4 trans = mat4.identity);
132 }
133