| 109 | } |
| 110 | |
| 111 | void simulate(af::array &mass, vector<af::array> &pos, vector<af::array> &vels, |
| 112 | vector<af::array> &forces, float dt) { |
| 113 | for (int i = 0; i < (int)pos.size(); ++i) { |
| 114 | pos[i] += vels[i] * dt; |
| 115 | pos[i].eval(); |
| 116 | } |
| 117 | |
| 118 | // calculate forces to each particle |
| 119 | vector<af::array> diff(pos.size()); |
| 120 | af::array dist = af::constant(0, pos[0].dims(0), pos[0].dims(0)); |
| 121 | |
| 122 | for (int i = 0; i < (int)pos.size(); ++i) { |
| 123 | diff[i] = tile(pos[i], 1, pos[i].dims(0)) - |
| 124 | transpose(tile(pos[i], 1, pos[i].dims(0))); |
| 125 | dist += (diff[i] * diff[i]); |
| 126 | } |
| 127 | |
| 128 | dist = sqrt(dist); |
| 129 | dist = af::max(min_dist, dist); |
| 130 | dist *= dist * dist; |
| 131 | |
| 132 | for (int i = 0; i < (int)pos.size(); ++i) { |
| 133 | // calculate force vectors |
| 134 | forces[i] = diff[i] / dist; |
| 135 | forces[i].eval(); |
| 136 | |
| 137 | // af::array idx = af::where(af::isNaN(forces[i])); |
| 138 | // if(idx.elements() > 0) |
| 139 | // forces[i](idx) = 0; |
| 140 | // forces[i] = sum(forces[i]).T(); |
| 141 | forces[i] = matmul(forces[i].T(), mass); |
| 142 | |
| 143 | // update force scaled to time, magnitude constant |
| 144 | forces[i] *= (gravity_constant); |
| 145 | forces[i].eval(); |
| 146 | |
| 147 | // update velocities from forces |
| 148 | vels[i] += forces[i] * dt; |
| 149 | vels[i].eval(); |
| 150 | |
| 151 | // noise |
| 152 | // forces[i] += 0.1 * af::randn(forces[i].dims(0)); |
| 153 | |
| 154 | // dampening |
| 155 | // vels[i] *= 1 - (0.005*dt); |
| 156 | } |
| 157 | } |
| 158 | |
| 159 | void collisions(vector<af::array> &pos, vector<af::array> &vels, bool is3D) { |
| 160 | // clamp particles inside screen border |