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