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.rigidbody;
30 
31 import std.math;
32 
33 import dlib.core.memory;
34 import dlib.container.array;
35 import dlib.math.vector;
36 import dlib.math.matrix;
37 import dlib.math.quaternion;
38 import dlib.math.transformation;
39 import dlib.math.utils;
40 
41 import dmech.shape;
42 import dmech.contact;
43 
44 interface CollisionDispatcher
45 {
46     void onNewContact(RigidBody rb, Contact c);
47 }
48 
49 class RigidBody: Freeable
50 {
51     Vector3f position;
52     Quaternionf orientation;
53 
54     Vector3f linearVelocity;
55     Vector3f angularVelocity;
56 
57     Vector3f pseudoLinearVelocity;
58     Vector3f pseudoAngularVelocity;
59 
60     Vector3f linearAcceleration;
61     Vector3f angularAcceleration;
62 
63     float mass;
64     float invMass;
65 
66     Matrix3x3f inertiaTensor;
67     Matrix3x3f invInertiaTensor; //TODO: rename to invInertiaWorld
68 
69     Vector3f centerOfMass;
70 
71     Vector3f forceAccumulator;
72     Vector3f torqueAccumulator;
73 
74     float bounce;
75     float friction;
76 
77     bool active = true;
78     bool useGravity = true;
79     bool enableRotation = true;
80 
81     DynamicArray!ShapeComponent shapes;
82 
83     bool dynamic;
84 
85     float damping = 0.5f;
86     float stopThreshold = 0.15f; //0.15f
87     float stopThresholdPV = 0.0f; //0.01f
88 
89     bool useOwnGravity = false;
90     Vector3f gravity = Vector3f(0, 0, 0);
91 
92     DynamicArray!CollisionDispatcher collisionDispatchers;
93 
94     void contactEvent(Contact c)
95     {
96         foreach(d; collisionDispatchers.data)
97         {
98             d.onNewContact(this, c);
99         }
100     }
101 
102     bool raycastable = true;
103     uint numContacts = 0;
104     bool useFriction = true;
105     
106     float maxSpeed = float.max;
107 
108     this()
109     {
110         position = Vector3f(0.0f, 0.0f, 0.0f);
111         orientation = Quaternionf(0.0f, 0.0f, 0.0f, 1.0f);
112 
113         linearVelocity = Vector3f(0.0f, 0.0f, 0.0f);
114         angularVelocity = Vector3f(0.0f, 0.0f, 0.0f);
115 
116         pseudoLinearVelocity = Vector3f(0.0f, 0.0f, 0.0f);
117         pseudoAngularVelocity = Vector3f(0.0f, 0.0f, 0.0f);
118 
119         linearAcceleration = Vector3f(0.0f, 0.0f, 0.0f);
120         angularAcceleration = Vector3f(0.0f, 0.0f, 0.0f);
121 
122         mass = 1.0f;
123         invMass = 1.0f;
124 
125         inertiaTensor = matrixf(
126             mass, 0, 0,
127             0, mass, 0,
128             0, 0, mass
129         );
130         invInertiaTensor = matrixf(
131             0, 0, 0,
132             0, 0, 0,
133             0, 0, 0
134         );
135 
136         centerOfMass = Vector3f(0.0f, 0.0f, 0.0f);
137 
138         forceAccumulator = Vector3f(0.0f, 0.0f, 0.0f);
139         torqueAccumulator = Vector3f(0.0f, 0.0f, 0.0f);
140 
141         bounce = 0.0f;
142         friction = 0.5f;
143 
144         dynamic = true;
145     }
146 
147     Matrix4x4f transformation()
148     {
149         Matrix4x4f t;
150         t = translationMatrix(position);
151         t *= orientation.toMatrix4x4();
152         return t;
153     }
154 
155     void addShapeComponent(ShapeComponent shape)
156     {
157         shape.transformation = transformation() * translationMatrix(shape.centroid);
158         shapes.append(shape);
159 
160         if (!dynamic)
161             return;
162 
163         centerOfMass = Vector3f(0, 0, 0);
164         float m = 0.0f;
165 
166         foreach (sh; shapes.data)
167         {
168             m += sh.mass;
169             centerOfMass += sh.mass * sh.centroid;
170         }
171 
172         float invm = 1.0f / m;
173         centerOfMass *= invm;
174 
175         mass += shape.mass;
176         invMass = 1.0f / mass;
177 
178         // Compute inertia tensor using Huygens-Steiner theorem
179         inertiaTensor = Matrix3x3f.zero;
180         foreach (sh; shapes.data)
181         {
182             Vector3f r = centerOfMass - sh.centroid;
183             inertiaTensor +=
184                 sh.geometry.inertiaTensor(sh.mass) +
185                 (Matrix3x3f.identity * dot(r, r) - tensorProduct(r, r)) * sh.mass;
186         }
187 
188         invInertiaTensor = inertiaTensor.inverse;
189     }
190 
191     void integrateForces(float dt)
192     {
193         if (!active || !dynamic)
194             return;
195 
196         linearAcceleration = forceAccumulator * invMass;
197         if (enableRotation)
198             angularAcceleration = torqueAccumulator * invInertiaTensor;
199 
200         linearVelocity += linearAcceleration * dt;
201         if (enableRotation)
202             angularVelocity += angularAcceleration * dt;
203     }
204 
205     void integrateVelocities(float dt)
206     {
207         if (!active || !dynamic)
208             return;
209 
210         float d = clamp(1.0f - dt * damping, 0.0f, 1.0f);
211         linearVelocity *= d;
212         angularVelocity *= d;
213         
214         float speed = linearVelocity.length;
215         
216         if (speed > maxSpeed)
217             linearVelocity = linearVelocity.normalized * maxSpeed;
218 
219         if (speed > stopThreshold /* || numContacts < 3 */)
220         {
221             position += linearVelocity * dt;
222         }
223 
224         if (enableRotation)
225         if (angularVelocity.length > 0.2f /* || numContacts < 3 */) //stopThreshold
226         {
227             orientation += 0.5f * Quaternionf(angularVelocity, 0.0f) * orientation * dt;
228             orientation.normalize();
229         }
230     }
231 
232     void integratePseudoVelocities(float dt)
233     {
234         if (!active || !dynamic)
235             return;
236 
237         float d = clamp(1.0f - dt * damping, 0.0f, 1.0f);
238 
239         pseudoLinearVelocity *= d;
240         pseudoAngularVelocity *= d;
241 
242         if (pseudoLinearVelocity.length > stopThresholdPV)
243         {
244             position += pseudoLinearVelocity;
245         }
246 
247         if (enableRotation)
248         if (pseudoAngularVelocity.length > stopThresholdPV)
249         {
250             orientation += 0.5f * Quaternionf(pseudoAngularVelocity, 0.0f) * orientation;
251             orientation.normalize();
252         }
253 
254         pseudoLinearVelocity = Vector3f(0.0f, 0.0f, 0.0f);
255         pseudoAngularVelocity = Vector3f(0.0f, 0.0f, 0.0f);
256     }
257 
258     @property Vector3f worldCenterOfMass()
259     {
260         return position + orientation.rotate(centerOfMass);
261     }
262 
263     void updateInertia()
264     {
265         if (active && dynamic)
266         {
267             auto rot = orientation.toMatrix3x3;
268             //invInertiaTensor = (rot * inertiaTensor * rot.transposed).inverse;
269             invInertiaTensor = (rot * inertiaTensor.inverse * rot.transposed);
270         }
271     }
272 
273     void updateShapeComponents()
274     {
275         foreach (sh; shapes.data)
276         {
277             sh.transformation = transformation() * translationMatrix(sh.centroid);
278         }
279     }
280 
281     void resetForces()
282     {
283         forceAccumulator = Vector3f(0, 0, 0);
284         torqueAccumulator = Vector3f(0, 0, 0);
285     }
286 
287     void applyForce(Vector3f force)
288     {
289         if (!active || !dynamic)
290             return;
291 
292         forceAccumulator += force;
293     }
294 
295     void applyTorque(Vector3f torque)
296     {
297         if (!active || !dynamic)
298             return;
299 
300         torqueAccumulator += torque;
301     }
302 
303     void applyImpulse(Vector3f impulse, Vector3f point)
304     {
305         if (!active || !dynamic)
306             return;
307 
308         linearVelocity += impulse * invMass;
309         Vector3f angularImpulse = cross(point - worldCenterOfMass, impulse);
310         angularVelocity += angularImpulse * invInertiaTensor;
311     }
312 
313     void applyForceAtPos(Vector3f force, Vector3f pos)
314     {
315         if (!active || !dynamic)
316             return;
317 
318         forceAccumulator += force;
319         torqueAccumulator += cross(pos - worldCenterOfMass, force);
320     }
321 
322     void applyPseudoImpulse(Vector3f impulse, Vector3f point)
323     {
324         if (!active || !dynamic)
325             return;
326 
327         pseudoLinearVelocity += impulse * invMass;
328         Vector3f angularImpulse = cross(point - worldCenterOfMass, impulse);
329         pseudoAngularVelocity += angularImpulse * invInertiaTensor;
330     }
331 
332     @property float linearKineticEnergy()
333     {
334         if (!active || !dynamic)
335             return float.infinity;
336 
337         // 0.5 * m * v^2
338         float v = linearVelocity.length;
339         return 0.5f * mass * v * v;
340     }
341 
342     @property float angularKineticEnergy()
343     {
344         if (!active || !dynamic)
345             return float.infinity;
346 
347         // 0.5 * w * I * w
348         Vector3f w = angularVelocity;
349         return 0.5f * dot(w * inertiaTensor, w);
350     }
351 
352     ~this()
353     {
354         shapes.free();
355         collisionDispatchers.free();
356     }
357 
358     void free()
359     {
360         Delete(this);
361     }
362 }