8dc666d588f54524671cd8cd7a15547496f0fded
[blender.git] / extern / bullet / BulletDynamics / ConstraintSolver / OdeConstraintSolver.cpp
1 /*
2  * Copyright (c) 2005 Erwin Coumans http://www.erwincoumans.com
3  *
4  * Permission to use, copy, modify, distribute and sell this software
5  * and its documentation for any purpose is hereby granted without fee,
6  * provided that the above copyright notice appear in all copies.
7  * Erwin Coumans makes no representations about the suitability 
8  * of this software for any purpose.  
9  * It is provided "as is" without express or implied warranty.
10 */
11
12 #include "OdeConstraintSolver.h"
13
14 #include "NarrowPhaseCollision/PersistentManifold.h"
15 #include "Dynamics/RigidBody.h"
16 #include "ContactConstraint.h"
17 #include "Solve2LinearConstraint.h"
18 #include "ContactSolverInfo.h"
19 #include "Dynamics/BU_Joint.h"
20 #include "Dynamics/ContactJoint.h"
21
22 #define USE_SOR_SOLVER
23
24 #include "SorLcp.h"
25
26 #include <math.h>
27 #include <float.h>//FLT_MAX
28 #ifdef WIN32
29 #include <memory.h>
30 #endif
31 #include <string.h>
32 #include <stdio.h>
33
34 #ifdef WIN32
35 #include <malloc.h>
36 #else
37 #include <alloca.h>
38 #endif
39
40 class BU_Joint;
41
42 //see below
43
44 int gCurBody = 0;
45 int gCurJoint= 0;
46
47 int ConvertBody(RigidBody* body,RigidBody** bodies,int& numBodies);
48 void ConvertConstraint(PersistentManifold* manifold,BU_Joint** joints,int& numJoints,
49                                            RigidBody** bodies,int bodyId0,int bodyId1);
50
51
52
53
54
55 //iterative lcp and penalty method
56 float OdeConstraintSolver::SolveGroup(PersistentManifold** manifoldPtr, int numManifolds,const ContactSolverInfo& infoGlobal)
57 {
58         gCurBody = 0;
59         gCurJoint = 0;
60
61         float cfm = 1e-5f;
62         float erp = 0.2f;
63
64         RigidBody* bodies [128];
65
66         int numBodies = 0;
67         BU_Joint* joints [128*5];
68         int numJoints = 0;
69
70         for (int j=0;j<numManifolds;j++)
71         {
72
73                 int body0=-1,body1=-1;
74
75                 PersistentManifold* manifold = manifoldPtr[j];
76                 if (manifold->GetNumContacts() > 0)
77                 {
78                         body0 = ConvertBody((RigidBody*)manifold->GetBody0(),bodies,numBodies);
79                         body1 = ConvertBody((RigidBody*)manifold->GetBody1(),bodies,numBodies);
80                         ConvertConstraint(manifold,joints,numJoints,bodies,body0,body1);
81                 }
82         }
83
84         SolveInternal1(cfm,erp,bodies,numBodies,joints,numJoints,infoGlobal);
85
86         return 0.f;
87
88 }
89
90 /////////////////////////////////////////////////////////////////////////////////
91
92
93 typedef SimdScalar dQuaternion[4];
94 #define _R(i,j) R[(i)*4+(j)]
95
96 void dRfromQ1 (dMatrix3 R, const dQuaternion q)
97 {
98   // q = (s,vx,vy,vz)
99   SimdScalar qq1 = 2*q[1]*q[1];
100   SimdScalar qq2 = 2*q[2]*q[2];
101   SimdScalar qq3 = 2*q[3]*q[3];
102   _R(0,0) = 1 - qq2 - qq3;
103   _R(0,1) = 2*(q[1]*q[2] - q[0]*q[3]);
104   _R(0,2) = 2*(q[1]*q[3] + q[0]*q[2]);
105   _R(1,0) = 2*(q[1]*q[2] + q[0]*q[3]);
106   _R(1,1) = 1 - qq1 - qq3;
107   _R(1,2) = 2*(q[2]*q[3] - q[0]*q[1]);
108   _R(2,0) = 2*(q[1]*q[3] - q[0]*q[2]);
109   _R(2,1) = 2*(q[2]*q[3] + q[0]*q[1]);
110   _R(2,2) = 1 - qq1 - qq2;
111 }
112
113
114
115 int ConvertBody(RigidBody* body,RigidBody** bodies,int& numBodies)
116 {
117         if (!body || (body->getInvMass() == 0.f) )
118         {
119                 return -1;
120         }
121         //first try to find
122         int i,j;
123         for (i=0;i<numBodies;i++)
124         {
125                 if (bodies[i] == body)
126                         return i;
127         }
128         //if not found, create a new body
129         bodies[numBodies++] = body;
130         //convert data
131
132
133         body->m_facc.setValue(0,0,0);
134         body->m_tacc.setValue(0,0,0);
135         
136         //are the indices the same ?
137         for (i=0;i<4;i++)
138         {
139                 for ( j=0;j<3;j++)
140                 {
141                         body->m_invI[i+4*j] = 0.f;
142                         body->m_I[i+4*j] = 0.f;
143                 }
144         }
145         body->m_invI[0+4*0] =   body->getInvInertiaDiagLocal()[0];
146         body->m_invI[1+4*1] =   body->getInvInertiaDiagLocal()[1];
147         body->m_invI[2+4*2] =   body->getInvInertiaDiagLocal()[2];
148
149         body->m_I[0+0*4] = body->getInvInertiaDiagLocal()[0];
150         body->m_I[1+1*4] = body->getInvInertiaDiagLocal()[1];
151         body->m_I[2+2*4] = body->getInvInertiaDiagLocal()[2];
152         
153
154         /*
155         
156         SimdMatrix3x3 invI;
157         invI.setIdentity();
158         invI[0][0] = body->getInvInertiaDiagLocal()[0];
159         invI[1][1] = body->getInvInertiaDiagLocal()[1];
160         invI[2][2] = body->getInvInertiaDiagLocal()[2];
161         SimdMatrix3x3 inertia = invI.inverse();
162
163         for (i=0;i<3;i++)
164         {
165                 for (j=0;j<3;j++)
166                 {
167                         body->m_I[i+4*j] = inertia[i][j];
168                 }
169         }
170         */
171 //      body->m_I[3+0*4] = 0.f;
172 //      body->m_I[3+1*4] = 0.f;
173 //      body->m_I[3+2*4] = 0.f;
174 //      body->m_I[3+3*4] = 0.f;
175         
176         
177         dQuaternion q;
178
179         q[1] = body->getOrientation()[0];
180         q[2] = body->getOrientation()[1];
181         q[3] = body->getOrientation()[2];
182         q[0] = body->getOrientation()[3];
183         
184         dRfromQ1(body->m_R,q);
185         
186         return numBodies-1;
187 }
188
189
190
191
192         
193 #define MAX_JOINTS_1 8192
194
195 static ContactJoint gJointArray[MAX_JOINTS_1];
196
197
198 void ConvertConstraint(PersistentManifold* manifold,BU_Joint** joints,int& numJoints,
199                                            RigidBody** bodies,int _bodyId0,int _bodyId1)
200 {
201
202
203         manifold->RefreshContactPoints(((RigidBody*)manifold->GetBody0())->getCenterOfMassTransform(),
204                 ((RigidBody*)manifold->GetBody1())->getCenterOfMassTransform());
205
206         int bodyId0 = _bodyId0,bodyId1 = _bodyId1;
207
208         int i,numContacts = manifold->GetNumContacts();
209         
210         bool swapBodies = (bodyId0 < 0);
211
212         
213         RigidBody* body0,*body1;
214
215         if (swapBodies)
216         {
217                 bodyId0 = _bodyId1;
218                 bodyId1 = _bodyId0;
219
220                 body0 = (RigidBody*)manifold->GetBody1();
221                 body1 = (RigidBody*)manifold->GetBody0();
222
223         } else
224         {
225                 body0 = (RigidBody*)manifold->GetBody0();
226                 body1 = (RigidBody*)manifold->GetBody1();
227         }
228
229         assert(bodyId0 >= 0);
230
231         for (i=0;i<numContacts;i++)
232         {
233                 
234                 assert (gCurJoint < MAX_JOINTS_1);
235
236                 if (manifold->GetContactPoint(i).GetDistance() < 0.f)
237                 {
238                         ContactJoint* cont = new (&gJointArray[gCurJoint++]) ContactJoint( manifold ,i, swapBodies,body0,body1);
239
240                         cont->node[0].joint = cont;
241                         cont->node[0].body = bodyId0 >= 0 ? bodies[bodyId0] : 0;
242                         
243                         cont->node[1].joint = cont;
244                         cont->node[1].body = bodyId1 >= 0 ? bodies[bodyId1] : 0;
245                         
246                         joints[numJoints++] = cont;
247                         for (int i=0;i<6;i++)
248                                 cont->lambda[i] = 0.f;
249
250                         cont->flags = 0;
251                 }
252         }
253
254         //create a new contact constraint
255 };
256