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