1 /* 2 Copyright (c) 2013-2014 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.world; 30 31 import std.math; 32 import std.range; 33 34 import dlib.math.vector; 35 import dlib.math.matrix; 36 import dlib.math.affine; 37 import dlib.geometry.triangle; 38 import dlib.geometry.sphere; 39 import dlib.geometry.ray; 40 41 import dmech.rigidbody; 42 import dmech.geometry; 43 import dmech.shape; 44 import dmech.contact; 45 import dmech.solver; 46 import dmech.pairhashtable; 47 import dmech.collision; 48 import dmech.pcm; 49 import dmech.constraint; 50 import dmech.bvh; 51 import dmech.mpr; 52 import dmech.raycast; 53 54 /* 55 * World object stores bodies and constraints and performs 56 * simulation cycles on them. 57 */ 58 59 class PhysicsWorld 60 { 61 RigidBody[] staticBodies; 62 RigidBody[] dynamicBodies; 63 Constraint[] constraints; 64 65 Vector3f gravity; 66 67 protected uint maxShapeId = 1; 68 69 PairHashTable!PersistentContactManifold manifolds; 70 71 bool broadphase = false; 72 bool warmstart = false; 73 74 uint positionCorrectionIterations = 10; 75 uint constraintIterations = 20; 76 77 BVHNode!Triangle bvhRoot = null; 78 79 // Proxy triangle to deal with BVH data 80 RigidBody proxyTri; 81 ShapeComponent proxyTriShape; 82 GeomTriangle proxyTriGeom; 83 84 this(size_t maxCollisions = 1000) 85 { 86 gravity = Vector3f(0.0f, -9.80665f, 0.0f); // Earth -9.80665f 87 88 manifolds = new PairHashTable!PersistentContactManifold(maxCollisions); 89 90 // Create proxy triangle 91 proxyTri = new RigidBody(); 92 proxyTri.position = Vector3f(0, 0, 0); 93 proxyTriGeom = new GeomTriangle( 94 Vector3f(-1.0f, 0.0f, -1.0f), 95 Vector3f(+1.0f, 0.0f, 0.0f), 96 Vector3f(-1.0f, 0.0f, +1.0f)); 97 proxyTriShape = new ShapeComponent(proxyTriGeom, Vector3f(0, 0, 0), 1); 98 proxyTriShape.id = maxShapeId; 99 maxShapeId++; 100 proxyTriShape.transformation = 101 proxyTri.transformation() * translationMatrix(proxyTriShape.centroid); 102 proxyTri.shapes ~= proxyTriShape; 103 proxyTri.mass = float.infinity; 104 proxyTri.invMass = 0.0f; 105 proxyTri.inertiaTensor = matrixf( 106 float.infinity, 0, 0, 107 0, float.infinity, 0, 108 0, 0, float.infinity 109 ); 110 proxyTri.invInertiaTensor = matrixf( 111 0, 0, 0, 112 0, 0, 0, 113 0, 0, 0 114 ); 115 proxyTri.dynamic = false; 116 } 117 118 RigidBody addDynamicBody(Vector3f pos, float mass = 0.0f) 119 { 120 auto b = new RigidBody(); 121 b.position = pos; 122 b.mass = mass; 123 b.invMass = 1.0f / mass; 124 b.inertiaTensor = matrixf( 125 mass, 0, 0, 126 0, mass, 0, 127 0, 0, mass 128 ); 129 b.invInertiaTensor = matrixf( 130 0, 0, 0, 131 0, 0, 0, 132 0, 0, 0 133 ); 134 b.dynamic = true; 135 dynamicBodies ~= b; 136 return b; 137 } 138 139 RigidBody addStaticBody(Vector3f pos) 140 { 141 auto b = new RigidBody(); 142 b.position = pos; 143 b.mass = float.infinity; 144 b.invMass = 0.0f; 145 b.inertiaTensor = matrixf( 146 float.infinity, 0, 0, 147 0, float.infinity, 0, 148 0, 0, float.infinity 149 ); 150 b.invInertiaTensor = matrixf( 151 0, 0, 0, 152 0, 0, 0, 153 0, 0, 0 154 ); 155 b.dynamic = false; 156 staticBodies ~= b; 157 return b; 158 } 159 160 ShapeComponent addShapeComponent(RigidBody b, Geometry geom, Vector3f position, float mass) 161 { 162 auto shape = new ShapeComponent(geom, position, mass); 163 shape.id = maxShapeId; 164 maxShapeId++; 165 b.addShapeComponent(shape); 166 return shape; 167 } 168 169 Constraint addConstraint(Constraint c) 170 { 171 constraints ~= c; 172 return c; 173 } 174 175 void update(double dt) 176 { 177 if (dynamicBodies.length == 0) 178 return; 179 180 foreach(ref m; manifolds) 181 { 182 m.update(); 183 } 184 185 foreach(b; dynamicBodies) 186 { 187 b.updateInertia(); 188 if (b.useGravity) 189 { 190 if (b.useOwnGravity) 191 b.applyForce(b.gravity * b.mass); 192 else 193 b.applyForce(gravity * b.mass); 194 } 195 b.integrateForces(dt); 196 b.resetForces(); 197 } 198 199 if (broadphase) 200 findDynamicCollisionsBroadphase(); 201 else 202 findDynamicCollisionsBruteForce(); 203 204 findStaticCollisionsBruteForce(); 205 206 solveConstraints(dt); 207 208 foreach(b; dynamicBodies) 209 { 210 b.integrateVelocities(dt); 211 } 212 213 foreach(iteration; 0..positionCorrectionIterations) 214 foreach(ref m; manifolds) 215 foreach(i; 0..m.numContacts) 216 { 217 auto c = &m.contacts[i]; 218 solvePositionError(c, m.numContacts); 219 } 220 221 foreach(b; dynamicBodies) 222 { 223 b.integratePseudoVelocities(dt); 224 b.updateShapeComponents(); 225 } 226 } 227 228 bool raycast( 229 Vector3f rayStart, 230 Vector3f rayDir, 231 float maxRayDist, 232 out CastResult castResult, 233 bool checkAgainstBodies = true, 234 bool checkAgainstBVH = true) 235 { 236 bool res = false; 237 float bestParam = float.max; 238 239 CastResult cr; 240 241 if (checkAgainstBodies) 242 foreach(b; chain(staticBodies, dynamicBodies)) 243 if (b.raycastable) 244 foreach(shape; b.shapes) 245 { 246 bool hit = convexRayCast(shape, rayStart, rayDir, maxRayDist, cr); 247 if (hit) 248 { 249 if (cr.param < bestParam) 250 { 251 bestParam = cr.param; 252 castResult = cr; 253 castResult.rbody = b; 254 res = true; 255 } 256 } 257 } 258 259 if (!checkAgainstBVH) 260 return res; 261 262 Ray ray = Ray(rayStart, rayStart + rayDir * maxRayDist); 263 264 if (bvhRoot !is null) 265 bvhRoot.traverseByRay(ray, (ref Triangle tri) 266 { 267 Vector3f ip; 268 bool hit = ray.intersectTriangle(tri.v[0], tri.v[1], tri.v[2], ip); 269 if (hit) 270 { 271 float param = distance(rayStart, ip); 272 if (param < bestParam) 273 { 274 bestParam = param; 275 castResult.fact = true; 276 castResult.param = param; 277 castResult.point = ip; 278 castResult.normal = tri.normal; 279 castResult.rbody = proxyTri; 280 castResult.shape = null; 281 res = true; 282 } 283 } 284 }); 285 286 return res; 287 } 288 289 void findDynamicCollisionsBruteForce() 290 { 291 for (int i = 0; i < dynamicBodies.length - 1; i++) 292 { 293 auto body1 = dynamicBodies[i]; 294 foreach(shape1; body1.shapes) 295 { 296 for (int j = i + 1; j < dynamicBodies.length; j++) 297 { 298 auto body2 = dynamicBodies[j]; 299 foreach(shape2; body2.shapes) 300 { 301 Contact c; 302 c.body1 = body1; 303 c.body2 = body2; 304 checkCollisionPair(shape1, shape2, c); 305 } 306 } 307 } 308 } 309 } 310 311 void findDynamicCollisionsBroadphase() 312 { 313 for (int i = 0; i < dynamicBodies.length - 1; i++) 314 { 315 auto body1 = dynamicBodies[i]; 316 foreach(shape1; body1.shapes) 317 { 318 for (int j = i + 1; j < dynamicBodies.length; j++) 319 { 320 auto body2 = dynamicBodies[j]; 321 foreach(shape2; body2.shapes) 322 if (shape1.boundingBox.intersectsAABB(shape2.boundingBox)) 323 { 324 Contact c; 325 c.body1 = body1; 326 c.body2 = body2; 327 checkCollisionPair(shape1, shape2, c); 328 } 329 } 330 } 331 } 332 } 333 334 void findStaticCollisionsBruteForce() 335 { 336 foreach(body1; dynamicBodies) 337 { 338 foreach(shape1; body1.shapes) 339 { 340 foreach(body2; staticBodies) 341 { 342 foreach(shape2; body2.shapes) 343 { 344 Contact c; 345 c.body1 = body1; 346 c.body2 = body2; 347 c.shape2pos = shape2.position; 348 checkCollisionPair(shape1, shape2, c); 349 } 350 } 351 } 352 } 353 354 // Find collisions between dynamic bodies 355 // and the BVH world (static triangle mesh) 356 if (bvhRoot !is null) 357 foreach(rb; dynamicBodies) 358 foreach(shape; rb.shapes) 359 { 360 // There may be more than one contact at a time 361 static Contact[5] contacts; 362 static Triangle[5] contactTris; 363 uint numContacts = 0; 364 365 Contact c; 366 c.body1 = rb; 367 c.body2 = proxyTri; 368 c.fact = false; 369 370 Sphere sphere = shape.boundingSphere; 371 372 bvhRoot.traverseBySphere(sphere, (ref Triangle tri) 373 { 374 // Update temporary triangle to check collision 375 proxyTriShape.transformation = translationMatrix(tri.barycenter); 376 proxyTriGeom.v[0] = tri.v[0] - tri.barycenter; 377 proxyTriGeom.v[1] = tri.v[1] - tri.barycenter; 378 proxyTriGeom.v[2] = tri.v[2] - tri.barycenter; 379 380 bool collided = checkCollision(shape, proxyTriShape, c); 381 382 if (collided) 383 { 384 if (numContacts < contacts.length) 385 { 386 c.shape1RelPoint = c.point - shape.position; 387 c.shape2RelPoint = c.point - tri.barycenter; 388 c.body1RelPoint = c.point - c.body1.worldCenterOfMass; 389 c.body2RelPoint = c.point - tri.barycenter; 390 c.shape1 = shape; 391 c.shape2 = proxyTriShape; 392 c.shape2pos = tri.barycenter; 393 contacts[numContacts] = c; 394 contactTris[numContacts] = tri; 395 numContacts++; 396 } 397 } 398 }); 399 400 /* 401 * NOTE: 402 * There is a problem when rolling bodies over a triangle mesh. Instead of rolling 403 * straight it will get influenced when hitting triangle edges. 404 * Current solution is to solve only the contact with deepest penetration and 405 * throw out all others. Other possible approach is to merge all contacts that 406 * are within epsilon of each other. When merging the contacts, average and 407 * re-normalize the normals, and average the penetration depth value. 408 */ 409 410 int deepestContactIdx = -1; 411 float maxPen = 0.0f; 412 float bestGroundness = -1.0f; 413 foreach(i; 0..numContacts) 414 { 415 if (contacts[i].penetration > maxPen) 416 { 417 deepestContactIdx = i; 418 maxPen = contacts[i].penetration; 419 } 420 421 //Vector3f dirToContact = (contacts[i].point - rb.position).normalized; 422 //float groundness = dot(gravity.normalized, dirToContact); 423 //if (groundness > 0.7f) 424 // rb.onGround = true; 425 } 426 427 if (deepestContactIdx >= 0) 428 { 429 auto co = &contacts[deepestContactIdx]; 430 co.calcFDir(); 431 432 auto m = manifolds.get(shape.id, proxyTriShape.id); 433 if (m is null) 434 { 435 PersistentContactManifold m1; 436 m1.addContact(*co); 437 manifolds.set(shape.id, proxyTriShape.id, m1); 438 439 c.body1.contactEvent(c); 440 c.body2.contactEvent(c); 441 } 442 else 443 { 444 m.addContact(*co); 445 } 446 } 447 else 448 manifolds.remove(shape.id, proxyTriShape.id); 449 } 450 } 451 452 void checkCollisionPair(ShapeComponent shape1, ShapeComponent shape2, ref Contact c) 453 { 454 if (checkCollision(shape1, shape2, c)) 455 { 456 c.body1RelPoint = c.point - c.body1.worldCenterOfMass; 457 c.body2RelPoint = c.point - c.body2.worldCenterOfMass; 458 c.shape1RelPoint = c.point - shape1.position; 459 c.shape2RelPoint = c.point - shape2.position; 460 c.shape1 = shape1; 461 c.shape2 = shape2; 462 c.calcFDir(); 463 464 auto m = manifolds.get(shape1.id, shape2.id); 465 if (m is null) 466 { 467 PersistentContactManifold m1; 468 m1.addContact(c); 469 manifolds.set(shape1.id, shape2.id, m1); 470 471 c.body1.contactEvent(c); 472 c.body2.contactEvent(c); 473 } 474 else 475 { 476 m.addContact(c); 477 } 478 } 479 else 480 { 481 manifolds.remove(shape1.id, shape2.id); 482 } 483 } 484 485 void solveConstraints(double dt) 486 { 487 foreach(ref m; manifolds) 488 foreach(i; 0..m.numContacts) 489 { 490 auto c = &m.contacts[i]; 491 prepareContact(c); 492 } 493 494 foreach(c; constraints) 495 { 496 c.prepare(dt); 497 } 498 499 foreach(iteration; 0..constraintIterations) 500 { 501 foreach(c; constraints) 502 c.step(); 503 504 foreach(ref m; manifolds) 505 foreach(i; 0..m.numContacts) 506 { 507 auto c = &m.contacts[i]; 508 solveContact(c, dt); 509 } 510 } 511 } 512 } 513