PeriDyno 1.0.0
Loading...
Searching...
No Matches
TJConstraintSolver.cu
Go to the documentation of this file.
1#include "TJConstraintSolver.h"
2#include "SharedFuncsForRigidBody.h"
3
4namespace dyno
5{
6 IMPLEMENT_TCLASS(TJConstraintSolver, TDataType)
7
8 template<typename TDataType>
9 TJConstraintSolver<TDataType>::TJConstraintSolver()
10 :ConstraintModule()
11 {
12 this->inContacts()->tagOptional(true);
13 }
14
15 template<typename TDataType>
16 TJConstraintSolver<TDataType>::~TJConstraintSolver()
17 {
18 }
19
20 template<typename TDataType>
21 void TJConstraintSolver<TDataType>::initializeJacobian(Real dt)
22 {
23 int constraint_size = 0;
24 int contact_size = this->inContacts()->size();
25
26 auto topo = this->inDiscreteElements()->constDataPtr();
27
28 int ballAndSocketJoint_size = topo->ballAndSocketJoints().size();
29 int sliderJoint_size = topo->sliderJoints().size();
30 int hingeJoint_size = topo->hingeJoints().size();
31 int fixedJoint_size = topo->fixedJoints().size();
32 int pointJoint_size = topo->pointJoints().size();
33
34 if (this->varFrictionEnabled()->getData())
35 {
36 constraint_size += 3 * contact_size;
37 }
38 else
39 {
40 constraint_size = contact_size;
41 }
42
43 if (ballAndSocketJoint_size != 0)
44 {
45 constraint_size += 3 * ballAndSocketJoint_size;
46 }
47
48 if (sliderJoint_size != 0)
49 {
50 constraint_size += 8 * sliderJoint_size;
51 }
52
53 if (hingeJoint_size != 0)
54 {
55 constraint_size += 8 * hingeJoint_size;
56 }
57
58 if (fixedJoint_size != 0)
59 {
60 constraint_size += 6 * fixedJoint_size;
61 }
62
63 if (pointJoint_size != 0)
64 {
65 constraint_size += 3 * pointJoint_size;
66 }
67
68 if (constraint_size == 0)
69 {
70 return;
71 }
72
73 mVelocityConstraints.resize(constraint_size);
74
75 if (contact_size != 0)
76 {
77 auto& contacts = this->inContacts()->getData();
78 setUpContactAndFrictionConstraints(
79 mVelocityConstraints,
80 mContactsInLocalFrame,
81 this->inCenter()->getData(),
82 this->inRotationMatrix()->getData(),
83 this->varFrictionEnabled()->getData()
84 );
85 }
86
87 if (ballAndSocketJoint_size != 0)
88 {
89 auto& joints = topo->ballAndSocketJoints();
90 int begin_index = contact_size;
91
92 if (this->varFrictionEnabled()->getData())
93 {
94 begin_index += 2 * contact_size;
95 }
96
97 setUpBallAndSocketJointConstraints(
98 mVelocityConstraints,
99 joints,
100 this->inCenter()->getData(),
101 this->inRotationMatrix()->getData(),
102 begin_index
103 );
104 }
105
106 if (sliderJoint_size != 0)
107 {
108 auto& joints = topo->sliderJoints();
109 int begin_index = contact_size;
110
111 if (this->varFrictionEnabled()->getData())
112 {
113 begin_index += 2 * contact_size;
114 }
115 begin_index += 3 * ballAndSocketJoint_size;
116 setUpSliderJointConstraints(
117 mVelocityConstraints,
118 joints,
119 this->inCenter()->getData(),
120 this->inRotationMatrix()->getData(),
121 this->inQuaternion()->getData(),
122 begin_index
123 );
124 }
125
126 if (hingeJoint_size != 0)
127 {
128 auto& joints = topo->hingeJoints();
129 int begin_index = contact_size + 3 * ballAndSocketJoint_size + 8 * sliderJoint_size;
130 if (this->varFrictionEnabled()->getData())
131 {
132 begin_index += 2 * contact_size;
133 }
134 setUpHingeJointConstraints(
135 mVelocityConstraints,
136 joints,
137 this->inCenter()->getData(),
138 this->inRotationMatrix()->getData(),
139 this->inQuaternion()->getData(),
140 begin_index
141 );
142 }
143
144 if (fixedJoint_size != 0)
145 {
146 auto& joints = topo->fixedJoints();
147 int begin_index = contact_size + 3 * ballAndSocketJoint_size + 8 * sliderJoint_size + 8 * hingeJoint_size;
148 if (this->varFrictionEnabled()->getData())
149 {
150 begin_index += 2 * contact_size;
151 }
152 setUpFixedJointConstraints(
153 mVelocityConstraints,
154 joints,
155 this->inRotationMatrix()->getData(),
156 this->inQuaternion()->getData(),
157 begin_index
158 );;
159 }
160
161 if (pointJoint_size != 0)
162 {
163 auto& joints = topo->pointJoints();
164 int begin_index = contact_size + 3 * ballAndSocketJoint_size + 8 * sliderJoint_size + 8 * hingeJoint_size + 6 * fixedJoint_size;
165 if (this->varFrictionEnabled()->getData())
166 {
167 begin_index += 2 * contact_size;
168 }
169 setUpPointJointConstraints(
170 mVelocityConstraints,
171 joints,
172 this->inCenter()->getData(),
173 begin_index
174 );
175 }
176
177 auto sizeOfRigids = this->inCenter()->size();
178 mContactNumber.resize(sizeOfRigids);
179
180 mJ.resize(4 * constraint_size);
181 mB.resize(4 * constraint_size);
182 mK_1.resize(constraint_size);
183 mK_2.resize(constraint_size);
184 mK_3.resize(constraint_size);
185 mEta.resize(constraint_size);
186 mLambda.resize(constraint_size);
187
188 mJ.reset();
189 mB.reset();
190 mK_1.reset();
191 mK_2.reset();
192 mK_3.reset();
193 mEta.reset();
194 mLambda.reset();
195
196 mContactNumber.reset();
197
198 calculateJacobianMatrix(
199 mJ,
200 mB,
201 this->inCenter()->getData(),
202 this->inInertia()->getData(),
203 this->inMass()->getData(),
204 this->inRotationMatrix()->getData(),
205 mVelocityConstraints
206 );
207
208 calculateK(
209 mVelocityConstraints,
210 mJ,
211 mB,
212 this->inCenter()->getData(),
213 this->inInertia()->getData(),
214 this->inMass()->getData(),
215 mK_1,
216 mK_2,
217 mK_3
218 );
219
220 mErrors.resize(constraint_size);
221 mErrors.reset();
222
223
224 calculateEtaVectorForPJSBaumgarte(
225 mEta,
226 mJ,
227 this->inVelocity()->getData(),
228 this->inAngularVelocity()->getData(),
229 this->inCenter()->getData(),
230 this->inQuaternion()->getData(),
231 mVelocityConstraints,
232 mErrors,
233 this->varSlop()->getValue(),
234 this->varBaumgarteBias()->getValue(),
235 this->varSubStepping()->getValue(),
236 dt
237 );
238
239 if (contact_size != 0)
240 {
241 calculateContactPoints(
242 this->inContacts()->getData(),
243 mContactNumber);
244 }
245 }
246
247 template<typename TDataType>
248 void TJConstraintSolver<TDataType>::constrain()
249 {
250 uint bodyNum = this->inCenter()->size();
251
252 auto topo = this->inDiscreteElements()->constDataPtr();
253
254 mImpulseC.resize(bodyNum * 2);
255 mImpulseExt.resize(bodyNum * 2);
256 mImpulseC.reset();
257 mImpulseExt.reset();
258
259 Real dt = this->inTimeStep()->getData();
260
261
262 if (!this->inContacts()->isEmpty() || topo->totalJointSize() > 0)
263 {
264 if (mContactsInLocalFrame.size() != this->inContacts()->size()) {
265 mContactsInLocalFrame.resize(this->inContacts()->size());
266 }
267
268 setUpContactsInLocalFrame(
269 mContactsInLocalFrame,
270 this->inContacts()->getData(),
271 this->inCenter()->getData(),
272 this->inRotationMatrix()->getData()
273 );
274
275 Real dh = dt / this->varSubStepping()->getValue();
276
277 for (int i = 0; i < this->varSubStepping()->getValue(); i++)
278 {
279 if (this->varGravityEnabled()->getValue())
280 {
281 setUpGravity(
282 mImpulseExt,
283 this->varGravityValue()->getValue(),
284 dh
285 );
286 }
287
288
289 updateVelocity(
290 this->inAttribute()->getData(),
291 this->inVelocity()->getData(),
292 this->inAngularVelocity()->getData(),
293 mImpulseExt,
294 this->varLinearDamping()->getValue(),
295 this->varAngularDamping()->getValue(),
296 dh
297 );
298
299 mImpulseC.reset();
300 initializeJacobian(dh);
301 for (int j = 0; j < this->varIterationNumberForVelocitySolver()->getValue(); j++)
302 {
303 JacobiIteration(
304 mLambda,
305 mImpulseC,
306 mJ,
307 mB,
308 mEta,
309 mVelocityConstraints,
310 mContactNumber,
311 mK_1,
312 mK_2,
313 mK_3,
314 this->inMass()->getData(),
315 this->inFrictionCoefficients()->getData(),
316 this->varFrictionCoefficient()->getData(),
317 this->varGravityValue()->getData(),
318 dh
319 );
320 }
321
322 updateVelocity(
323 this->inAttribute()->getData(),
324 this->inVelocity()->getData(),
325 this->inAngularVelocity()->getData(),
326 mImpulseC,
327 this->varLinearDamping()->getValue(),
328 this->varAngularDamping()->getValue(),
329 dh
330 );
331
332 updateGesture(
333 this->inAttribute()->getData(),
334 this->inCenter()->getData(),
335 this->inQuaternion()->getData(),
336 this->inRotationMatrix()->getData(),
337 this->inInertia()->getData(),
338 this->inVelocity()->getData(),
339 this->inAngularVelocity()->getData(),
340 this->inInitialInertia()->getData(),
341 dh
342 );
343 }
344 }
345
346 else
347 {
348 if (this->varGravityEnabled()->getValue())
349 {
350 setUpGravity(
351 mImpulseExt,
352 this->varGravityValue()->getValue(),
353 dt
354 );
355 }
356
357
358 updateVelocity(
359 this->inAttribute()->getData(),
360 this->inVelocity()->getData(),
361 this->inAngularVelocity()->getData(),
362 mImpulseExt,
363 this->varLinearDamping()->getValue(),
364 this->varAngularDamping()->getValue(),
365 dt
366 );
367
368 updateGesture(
369 this->inAttribute()->getData(),
370 this->inCenter()->getData(),
371 this->inQuaternion()->getData(),
372 this->inRotationMatrix()->getData(),
373 this->inInertia()->getData(),
374 this->inVelocity()->getData(),
375 this->inAngularVelocity()->getData(),
376 this->inInitialInertia()->getData(),
377 dt
378 );
379 }
380
381 }
382
383 DEFINE_CLASS(TJConstraintSolver);
384}