1 /** 2 Inochi2D Physics Subsystem 3 4 Copyright: 5 Copyright © 2020-2026, Inochi2D Project 6 7 License: 8 $(LINK2 https://github.com/Inochi2D/inochi2d/blob/main/LICENSE, BSD 2-clause License) 9 10 Authors: 11 Hoshino Lina 12 */ 13 module inochi2d.core.phys.system; 14 import inochi2d; 15 import numath; 16 import numem; 17 import nulib; 18 19 abstract 20 class PhysicsSystem : NuObject { 21 private: 22 @nogc: 23 map!(float*, size_t) variableMap; 24 vector!(float*) refs; 25 26 // Derivatives 27 float[] kstate; 28 float[] k0; 29 float[] k1; 30 float[] k2; 31 float[] k3; 32 float[] k4; 33 34 float t; 35 36 /// Helper that resizes all of the temporary arrays. 37 void resize(size_t length) { 38 kstate = kstate.nu_resize(length); 39 k0 = k0.nu_resize(length); 40 k1 = k1.nu_resize(length); 41 k2 = k2.nu_resize(length); 42 k3 = k3.nu_resize(length); 43 k4 = k4.nu_resize(length); 44 } 45 46 protected: 47 /** 48 Add a float variable to the simulation 49 */ 50 size_t addVariable(float* var) { 51 size_t index = refs.length; 52 53 variableMap[var] = index; 54 refs ~= var; 55 56 this.resize(refs.length); 57 return index; 58 } 59 60 /** 61 Add a vec2 variable to the simulation 62 */ 63 size_t addVariable(vec2* var) { 64 size_t index = addVariable(&var.data[0]); 65 addVariable(&var.data[1]); 66 return index; 67 } 68 69 /** 70 Set the derivative of a variable (solver input) by index 71 */ 72 void setD(size_t index, float value) { 73 k0[index] = value; 74 } 75 76 /** 77 Set the derivative of a float variable (solver input) 78 */ 79 void setD(ref float var, float value) { 80 setD(variableMap[&var], value); 81 } 82 83 /** 84 Set the derivative of a vec2 variable (solver input) 85 */ 86 void setD(ref vec2 var, vec2 value) { 87 setD(var.data[0], value.x); 88 setD(var.data[1], value.y); 89 } 90 91 float[] getState() { 92 foreach (idx, ptr; refs) 93 kstate[idx] = *ptr; 94 return kstate; 95 } 96 97 void setState(float[] vals) { 98 foreach (idx, ptr; refs) { 99 *ptr = vals[idx]; 100 } 101 } 102 103 /** 104 Evaluate the simulation at a given time 105 */ 106 abstract void eval(float t); 107 108 public: 109 110 ~this() { 111 variableMap.clear(); 112 refs.clear(); 113 114 // Free derivatives. 115 nu_freea(kstate); 116 nu_freea(k0); 117 nu_freea(k1); 118 nu_freea(k2); 119 nu_freea(k3); 120 nu_freea(k4); 121 } 122 123 /** 124 Run a simulation tick (Runge-Kutta method) 125 */ 126 void tick(float h) { 127 float[] cur = getState(); 128 k0[0 .. $] = 0; 129 130 eval(t); 131 k1[0 .. $] = k0[0 .. $]; 132 133 foreach (i; 0 .. cur.length) 134 *refs[i] = cur[i] + h * k1[i] / 2f; 135 eval(t + h / 2f); 136 k2[0 .. $] = k0[0 .. $]; 137 138 foreach (i; 0 .. cur.length) 139 *refs[i] = cur[i] + h * k2[i] / 2f; 140 eval(t + h / 2f); 141 k3[0 .. $] = k0[0 .. $]; 142 143 foreach (i; 0 .. cur.length) 144 *refs[i] = cur[i] + h * k3[i]; 145 eval(t + h); 146 k4[0 .. $] = k0[0 .. $]; 147 148 foreach (i; 0 .. cur.length) { 149 *refs[i] = cur[i] + h * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6f; 150 if (!isFinite(*refs[i])) { 151 // Simulation failed, revert 152 foreach (j; 0 .. cur.length) 153 *refs[j] = cur[j]; 154 break; 155 } 156 } 157 158 t += h; 159 } 160 161 /** 162 Updates the anchor for the physics system 163 */ 164 abstract void updateAnchor(); 165 }