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.cm;
30 
31 import std.math;
32 
33 import dlib.math.vector;
34 import dlib.math.matrix;
35 import dlib.math.quaternion;
36 import dlib.math.utils;
37 import dlib.geometry.plane;
38 import dlib.geometry.triangle;
39 
40 import dmech.contact;
41 import dmech.rigidbody;
42 import dmech.geometry;
43 import dmech.shape;
44 import dmech.mpr;
45 import dmech.clipping;
46 
47 // TODO: move projectPointOnPlane and planeBasis to dlib.geometry.plane
48 /*
49 Vector3f projectPointOnPlane(Vector3f point, Vector3f planeOrigin, Vector3f planeNormal)
50 {
51     Vector3f v = point - planeOrigin;
52     float dist = dot(v, planeNormal);
53     Vector3f projectedPoint = point - dist * planeNormal;
54     return projectedPoint;
55 }
56 
57 // Z direction in this basis is plane normal
58 Matrix4x4f planeBasis(Vector3f planeOrigin, Vector3f planeNormal)
59 {
60     Matrix4x4f basis = directionToMatrix(planeNormal);
61     basis.translation = planeOrigin;
62     return basis;
63 }
64 */
65 
66 Vector2f projectPoint(
67     Vector3f point, 
68     Vector3f origin, 
69     Vector3f right,
70     Vector3f up)
71 {
72     Vector2f res;
73     res.x = dot(right, point - origin);
74     res.y = dot(up, point - origin);
75     return res;
76 }
77 
78 Vector3f unprojectPoint(
79     Vector2f point, 
80     Vector3f origin, 
81     Vector3f right,
82     Vector3f up)
83 {
84     Vector3f res;
85     res = origin + up * point.y + right * point.x;
86     return res;
87 }
88 
89 /*
90  * Contact manifold.
91  * Stores information about two bodies' collision.
92  * Contacts are generated at once using axis rotation method
93  * (for incremental solution see pcm.d).
94  */
95 struct ContactManifold
96 {
97     Contact[8] contacts;
98     uint numContacts = 0;
99     Contact mainContact;
100     Feature f1;
101     Feature f2;
102     Vector2f[8] pts;
103 
104     // Find all contact points at once
105     void computeContacts(Contact c)
106     {
107         mainContact = c;
108 
109         uint numPts = 0;
110 
111         // Calculate tangent space for contact normal
112         Vector3f n0, n1;
113         if (dot(c.normal, Vector3f(1,0,0)) < 0.5f)
114             n0 = cross(c.normal, Vector3f(1,0,0)); 
115         else
116             n0 = cross(c.normal, Vector3f(0,0,1));
117         n1 = cross(n0, c.normal);
118         n0.normalize();
119         n1.normalize();
120 
121         // If colliding with a sphere, there's only one contact
122         if (c.shape1.geometry.type == GeomType.Sphere ||
123             c.shape2.geometry.type == GeomType.Sphere)
124         {
125             contacts[0] = c;
126             contacts[0].fdir1 = n0;
127             contacts[0].fdir2 = n1;
128             numContacts = 1;
129             return;
130         }
131 
132         if (c.shape1.geometry.type == GeomType.Ellipsoid ||
133             c.shape2.geometry.type == GeomType.Ellipsoid)
134         {
135             contacts[0] = c;
136             contacts[0].fdir1 = n0;
137             contacts[0].fdir2 = n1;
138             numContacts = 1;
139             return;
140         }
141         
142         Vector3f right = n0;
143         Vector3f up = n1;
144         
145         // Calculate basis for contact plane
146         Plane contactPlane;
147         contactPlane.fromPointAndNormal(c.point, c.normal);
148 
149         // Scan contact features by rotating axis
150         f1.numVertices = 0;
151         f2.numVertices = 0;
152         
153         float eps = 0.05f;
154         
155         // If colliding with a cylinder, use smaller contact area
156         if (c.shape1.geometry.type == GeomType.Cylinder ||
157             c.shape2.geometry.type == GeomType.Cylinder)
158         {
159             eps = 0.01f;
160         }
161 
162         if (c.shape1.geometry.type == GeomType.Cone ||
163             c.shape2.geometry.type == GeomType.Cone)
164         {
165             eps = 0.01f;
166         }
167 
168         float startAng = PI / 4;
169         float stepAng = PI / 2; // 90 degrees
170 
171         uint numAxes1 = 4;
172 
173         if (c.shape1.geometry.type == GeomType.Triangle)
174         {
175             numAxes1 = 3;
176             startAng = 0;
177             stepAng = radtodeg(120.0f); // 2*PI/3 = 120 degrees
178         }
179 
180         for(uint i = 0; i < numAxes1; i++)
181         {
182             float ang = startAng + stepAng * i;
183 
184             Vector3f axis = (c.normal + n0 * cos(ang) * eps + n1 * sin(ang) * eps).normalized;
185             
186             Vector3f p;
187 
188             supportTransformed(c.shape1, -axis, p);
189             
190             if (contactPlane.distance(p) < 0.0f)
191             {
192                 Vector2f planePoint = projectPoint(p, c.point, right, up);
193                 f1.vertices[f1.numVertices] = planePoint;
194                 f1.numVertices++;
195             }
196         }
197 
198         uint numAxes2 = 4;
199 
200         if (c.shape2.geometry.type == GeomType.Triangle)
201         {
202             numAxes2 = 3;
203             startAng = 0;
204             stepAng = radtodeg(120.0f); // 2*PI/3 = 120 degrees
205         }
206         
207         for(uint i = 0; i < numAxes2; i++)
208         {
209             float ang = startAng + stepAng * i;
210 
211             Vector3f axis = (c.normal + n0 * cos(ang) * eps + n1 * sin(ang) * eps).normalized;
212             
213             Vector3f p;
214 
215             supportTransformed(c.shape2, axis, p);
216             if (contactPlane.distance(p) > 0.0f)
217             {
218                 Vector2f planePoint = projectPoint(p, c.point, right, up);
219                 f2.vertices[f2.numVertices] = planePoint;
220                 f2.numVertices++;
221             }
222         }
223 
224         // Clip features in 2D space: find their overlapping polygon
225         clip(f1, f2, pts, numPts);
226 
227         // Transform the resulting points back into 3D space       
228         Contact[8] newManifold;
229         Vector3f centroid = Vector3f(0.0f, 0.0f, 0.0f);
230         for(uint i = 0; i < numPts; i++)
231         {
232             Contact newc;
233             newc.fact = true;
234             newc.body1 = c.body1;
235             newc.body2 = c.body2;
236             newc.shape1 = c.shape1;
237             newc.shape2 = c.shape2;
238             newc.point = unprojectPoint(pts[i], c.point, right, up);
239             newc.normal = c.normal;
240             newc.penetration = c.penetration;
241             if (numPts > 1)
242                 centroid += newc.point;
243             else
244             {
245                 newc.fdir1 = n0;
246                 newc.fdir2 = n1;
247             }
248             newManifold[i] = newc;
249         }
250 
251         if (numPts > 1)
252         {
253             centroid /= numPts;
254 
255             for(uint i = 0; i < numPts; i++)
256             {
257                 newManifold[i].fdir1 = (newManifold[i].point - centroid).normalized;
258                 newManifold[i].fdir2 = cross(newManifold[i].fdir1, c.normal);
259             }
260         }
261         
262         // Update the existing manifold
263         updateManifold(newManifold, numPts);
264     }
265     
266     void updateManifold(ref Contact[8] newManifold, uint numNewContacts)
267     {
268         numContacts = 0;
269     
270         Contact[8] res;
271         bool[8] used;
272         uint resNum = 0;
273         
274         foreach(i; 0..numContacts)
275         {        
276             Vector3f p1 = contacts[i].point;
277             
278             foreach(j; 0..numNewContacts)
279             {
280                 Vector3f p2 = newManifold[j].point;
281                 
282                 if (distance(p1, p2) < 0.1f)
283                 {
284                     res[resNum] = contacts[i];
285                     res[resNum].shape1 = newManifold[j].shape1;
286                     res[resNum].shape2 = newManifold[j].shape2;
287                     res[resNum].body1 = newManifold[j].body1;
288                     res[resNum].body2 = newManifold[j].body2;
289                     res[resNum].point = p2; //(p1 + p2) * 0.5f;
290                     res[resNum].normal = newManifold[j].normal;
291                     res[resNum].penetration = newManifold[j].penetration;
292                     res[resNum].fdir1 = newManifold[j].fdir1;
293                     res[resNum].fdir2 = newManifold[j].fdir2;
294                     resNum++;
295                     used[j] = true;
296                 }
297             }
298         }
299         
300         foreach(i; 0..numNewContacts)
301         {
302             if (!used[i])
303             {
304                 res[resNum] = newManifold[i];
305                 resNum++;
306             }
307         }
308         
309         foreach(i; 0..resNum)
310         {
311             contacts[i] = res[i];
312             numContacts++;
313         }
314     }
315 }
316