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