1 /* 2 Copyright (c) 2013-2017 Timur Gafarov 3 4 Boost Software License - Version 1.0 - August 17th, 2003 5 6 Permission is hereby granted, free of charge, to any person or organization 7 obtaining a copy of the software and accompanying documentation covered by 8 this license (the "Software") to use, reproduce, display, distribute, 9 execute, and transmit the Software, and to prepare derivative works of the 10 Software, and to permit third-parties to whom the Software is furnished to 11 do so, all subject to the following: 12 13 The copyright notices in the Software and this entire statement, including 14 the above license grant, this restriction and the following disclaimer, 15 must be included in all copies of the Software, in whole or in part, and 16 all derivative works of the Software, unless such copies or derivative 17 works are solely in the form of machine-executable object code generated by 18 a source language processor. 19 20 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 21 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 22 FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT 23 SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE 24 FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, 25 ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 26 DEALINGS IN THE SOFTWARE. 27 */ 28 29 module dmech.rigidbody; 30 31 import std.math; 32 33 import dlib.core.memory; 34 import dlib.container.array; 35 import dlib.math.vector; 36 import dlib.math.matrix; 37 import dlib.math.quaternion; 38 import dlib.math.transformation; 39 import dlib.math.utils; 40 41 import dmech.shape; 42 import dmech.contact; 43 44 interface CollisionDispatcher 45 { 46 void onNewContact(RigidBody rb, Contact c); 47 } 48 49 class RigidBody: Freeable 50 { 51 Vector3f position; 52 Quaternionf orientation; 53 54 Vector3f linearVelocity; 55 Vector3f angularVelocity; 56 57 Vector3f pseudoLinearVelocity; 58 Vector3f pseudoAngularVelocity; 59 60 Vector3f linearAcceleration; 61 Vector3f angularAcceleration; 62 63 float mass; 64 float invMass; 65 66 Matrix3x3f inertiaTensor; 67 Matrix3x3f invInertiaTensor; //TODO: rename to invInertiaWorld 68 69 Vector3f centerOfMass; 70 71 Vector3f forceAccumulator; 72 Vector3f torqueAccumulator; 73 74 float bounce; 75 float friction; 76 77 bool active = true; 78 bool useGravity = true; 79 bool enableRotation = true; 80 81 DynamicArray!ShapeComponent shapes; 82 83 bool dynamic; 84 85 float damping = 0.5f; 86 float stopThreshold = 0.15f; //0.15f 87 float stopThresholdPV = 0.0f; //0.01f 88 89 bool useOwnGravity = false; 90 Vector3f gravity = Vector3f(0, 0, 0); 91 92 DynamicArray!CollisionDispatcher collisionDispatchers; 93 94 void contactEvent(Contact c) 95 { 96 foreach(d; collisionDispatchers.data) 97 { 98 d.onNewContact(this, c); 99 } 100 } 101 102 bool raycastable = true; 103 uint numContacts = 0; 104 bool useFriction = true; 105 106 float maxSpeed = float.max; 107 108 this() 109 { 110 position = Vector3f(0.0f, 0.0f, 0.0f); 111 orientation = Quaternionf(0.0f, 0.0f, 0.0f, 1.0f); 112 113 linearVelocity = Vector3f(0.0f, 0.0f, 0.0f); 114 angularVelocity = Vector3f(0.0f, 0.0f, 0.0f); 115 116 pseudoLinearVelocity = Vector3f(0.0f, 0.0f, 0.0f); 117 pseudoAngularVelocity = Vector3f(0.0f, 0.0f, 0.0f); 118 119 linearAcceleration = Vector3f(0.0f, 0.0f, 0.0f); 120 angularAcceleration = Vector3f(0.0f, 0.0f, 0.0f); 121 122 mass = 1.0f; 123 invMass = 1.0f; 124 125 inertiaTensor = matrixf( 126 mass, 0, 0, 127 0, mass, 0, 128 0, 0, mass 129 ); 130 invInertiaTensor = matrixf( 131 0, 0, 0, 132 0, 0, 0, 133 0, 0, 0 134 ); 135 136 centerOfMass = Vector3f(0.0f, 0.0f, 0.0f); 137 138 forceAccumulator = Vector3f(0.0f, 0.0f, 0.0f); 139 torqueAccumulator = Vector3f(0.0f, 0.0f, 0.0f); 140 141 bounce = 0.0f; 142 friction = 0.5f; 143 144 dynamic = true; 145 } 146 147 Matrix4x4f transformation() 148 { 149 Matrix4x4f t; 150 t = translationMatrix(position); 151 t *= orientation.toMatrix4x4(); 152 return t; 153 } 154 155 void addShapeComponent(ShapeComponent shape) 156 { 157 shape.transformation = transformation() * translationMatrix(shape.centroid); 158 shapes.append(shape); 159 160 if (!dynamic) 161 return; 162 163 centerOfMass = Vector3f(0, 0, 0); 164 float m = 0.0f; 165 166 foreach (sh; shapes.data) 167 { 168 m += sh.mass; 169 centerOfMass += sh.mass * sh.centroid; 170 } 171 172 float invm = 1.0f / m; 173 centerOfMass *= invm; 174 175 mass += shape.mass; 176 invMass = 1.0f / mass; 177 178 // Compute inertia tensor using Huygens-Steiner theorem 179 inertiaTensor = Matrix3x3f.zero; 180 foreach (sh; shapes.data) 181 { 182 Vector3f r = centerOfMass - sh.centroid; 183 inertiaTensor += 184 sh.geometry.inertiaTensor(sh.mass) + 185 (Matrix3x3f.identity * dot(r, r) - tensorProduct(r, r)) * sh.mass; 186 } 187 188 invInertiaTensor = inertiaTensor.inverse; 189 } 190 191 void integrateForces(float dt) 192 { 193 if (!active || !dynamic) 194 return; 195 196 linearAcceleration = forceAccumulator * invMass; 197 if (enableRotation) 198 angularAcceleration = torqueAccumulator * invInertiaTensor; 199 200 linearVelocity += linearAcceleration * dt; 201 if (enableRotation) 202 angularVelocity += angularAcceleration * dt; 203 } 204 205 void integrateVelocities(float dt) 206 { 207 if (!active || !dynamic) 208 return; 209 210 float d = clamp(1.0f - dt * damping, 0.0f, 1.0f); 211 linearVelocity *= d; 212 angularVelocity *= d; 213 214 float speed = linearVelocity.length; 215 216 if (speed > maxSpeed) 217 linearVelocity = linearVelocity.normalized * maxSpeed; 218 219 if (speed > stopThreshold /* || numContacts < 3 */) 220 { 221 position += linearVelocity * dt; 222 } 223 224 if (enableRotation) 225 if (angularVelocity.length > 0.2f /* || numContacts < 3 */) //stopThreshold 226 { 227 orientation += 0.5f * Quaternionf(angularVelocity, 0.0f) * orientation * dt; 228 orientation.normalize(); 229 } 230 } 231 232 void integratePseudoVelocities(float dt) 233 { 234 if (!active || !dynamic) 235 return; 236 237 float d = clamp(1.0f - dt * damping, 0.0f, 1.0f); 238 239 pseudoLinearVelocity *= d; 240 pseudoAngularVelocity *= d; 241 242 if (pseudoLinearVelocity.length > stopThresholdPV) 243 { 244 position += pseudoLinearVelocity; 245 } 246 247 if (enableRotation) 248 if (pseudoAngularVelocity.length > stopThresholdPV) 249 { 250 orientation += 0.5f * Quaternionf(pseudoAngularVelocity, 0.0f) * orientation; 251 orientation.normalize(); 252 } 253 254 pseudoLinearVelocity = Vector3f(0.0f, 0.0f, 0.0f); 255 pseudoAngularVelocity = Vector3f(0.0f, 0.0f, 0.0f); 256 } 257 258 @property Vector3f worldCenterOfMass() 259 { 260 return position + orientation.rotate(centerOfMass); 261 } 262 263 void updateInertia() 264 { 265 if (active && dynamic) 266 { 267 auto rot = orientation.toMatrix3x3; 268 //invInertiaTensor = (rot * inertiaTensor * rot.transposed).inverse; 269 invInertiaTensor = (rot * inertiaTensor.inverse * rot.transposed); 270 } 271 } 272 273 void updateShapeComponents() 274 { 275 foreach (sh; shapes.data) 276 { 277 sh.transformation = transformation() * translationMatrix(sh.centroid); 278 } 279 } 280 281 void resetForces() 282 { 283 forceAccumulator = Vector3f(0, 0, 0); 284 torqueAccumulator = Vector3f(0, 0, 0); 285 } 286 287 void applyForce(Vector3f force) 288 { 289 if (!active || !dynamic) 290 return; 291 292 forceAccumulator += force; 293 } 294 295 void applyTorque(Vector3f torque) 296 { 297 if (!active || !dynamic) 298 return; 299 300 torqueAccumulator += torque; 301 } 302 303 void applyImpulse(Vector3f impulse, Vector3f point) 304 { 305 if (!active || !dynamic) 306 return; 307 308 linearVelocity += impulse * invMass; 309 Vector3f angularImpulse = cross(point - worldCenterOfMass, impulse); 310 angularVelocity += angularImpulse * invInertiaTensor; 311 } 312 313 void applyForceAtPos(Vector3f force, Vector3f pos) 314 { 315 if (!active || !dynamic) 316 return; 317 318 forceAccumulator += force; 319 torqueAccumulator += cross(pos - worldCenterOfMass, force); 320 } 321 322 void applyPseudoImpulse(Vector3f impulse, Vector3f point) 323 { 324 if (!active || !dynamic) 325 return; 326 327 pseudoLinearVelocity += impulse * invMass; 328 Vector3f angularImpulse = cross(point - worldCenterOfMass, impulse); 329 pseudoAngularVelocity += angularImpulse * invInertiaTensor; 330 } 331 332 @property float linearKineticEnergy() 333 { 334 if (!active || !dynamic) 335 return float.infinity; 336 337 // 0.5 * m * v^2 338 float v = linearVelocity.length; 339 return 0.5f * mass * v * v; 340 } 341 342 @property float angularKineticEnergy() 343 { 344 if (!active || !dynamic) 345 return float.infinity; 346 347 // 0.5 * w * I * w 348 Vector3f w = angularVelocity; 349 return 0.5f * dot(w * inertiaTensor, w); 350 } 351 352 ~this() 353 { 354 shapes.free(); 355 collisionDispatchers.free(); 356 } 357 358 void free() 359 { 360 Delete(this); 361 } 362 }