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