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 }