# | Line 1 | Line 1 | |
---|---|---|
1 | #include <iostream> | |
2 | < | #include <cstdlib> |
3 | < | #include <cmath> |
4 | < | |
2 | > | #include <stdlib.h> |
3 | > | #include <math.h> |
4 | > | #include "Rattle.hpp" |
5 | > | #include "Roll.hpp" |
6 | #ifdef IS_MPI | |
7 | #include "mpiSimulation.hpp" | |
8 | #include <unistd.h> | |
9 | #endif //is_mpi | |
10 | ||
11 | + | #ifdef PROFILE |
12 | + | #include "mdProfile.hpp" |
13 | + | #endif // profile |
14 | + | |
15 | #include "Integrator.hpp" | |
16 | #include "simError.h" | |
17 | ||
# | Line 25 | Line 30 | template<typename T> Integrator<T>::Integrator(SimInfo | |
30 | if (info->the_integrator != NULL){ | |
31 | delete info->the_integrator; | |
32 | } | |
28 | – | info->the_integrator = this; |
33 | ||
34 | nAtoms = info->n_atoms; | |
35 | + | integrableObjects = info->integrableObjects; |
36 | ||
37 | + | consFramework = new RollFramework(info); |
38 | + | |
39 | + | if(consFramework == NULL){ |
40 | + | sprintf(painCave.errMsg, |
41 | + | "Integrator::Intergrator() Error: Memory allocation error for RattleFramework" ); |
42 | + | painCave.isFatal = 1; |
43 | + | simError(); |
44 | + | } |
45 | + | |
46 | + | /* |
47 | // check for constraints | |
48 | ||
49 | constrainedA = NULL; | |
# | Line 41 | Line 56 | template<typename T> Integrator<T>::Integrator(SimInfo | |
56 | nConstrained = 0; | |
57 | ||
58 | checkConstraints(); | |
59 | + | */ |
60 | } | |
61 | ||
62 | template<typename T> Integrator<T>::~Integrator(){ | |
63 | + | if (consFramework != NULL) |
64 | + | delete consFramework; |
65 | + | /* |
66 | if (nConstrained){ | |
67 | delete[] constrainedA; | |
68 | delete[] constrainedB; | |
# | Line 52 | Line 71 | template<typename T> Integrator<T>::~Integrator(){ | |
71 | delete[] moved; | |
72 | delete[] oldPos; | |
73 | } | |
74 | + | */ |
75 | } | |
76 | ||
77 | + | /* |
78 | template<typename T> void Integrator<T>::checkConstraints(void){ | |
79 | isConstrained = 0; | |
80 | ||
# | Line 65 | Line 86 | template<typename T> void Integrator<T>::checkConstrai | |
86 | ||
87 | SRI** theArray; | |
88 | for (int i = 0; i < nMols; i++){ | |
89 | < | theArray = (SRI * *) molecules[i].getMyBonds(); |
89 | > | |
90 | > | theArray = (SRI * *) molecules[i].getMyBonds(); |
91 | for (int j = 0; j < molecules[i].getNBonds(); j++){ | |
92 | constrained = theArray[j]->is_constrained(); | |
93 | ||
# | Line 87 | Line 109 | template<typename T> void Integrator<T>::checkConstrai | |
109 | if (constrained){ | |
110 | dummy_plug = theArray[j]->get_constraint(); | |
111 | temp_con[nConstrained].set_a(dummy_plug->get_a()); | |
112 | < | temp_con[nConstrained].set_b(dummy_plug->get_b()); |
112 | > | temp_con[nConstrained].set_b(Dummy_plug->get_b()); |
113 | temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr()); | |
114 | ||
115 | nConstrained++; | |
# | Line 111 | Line 133 | template<typename T> void Integrator<T>::checkConstrai | |
133 | } | |
134 | } | |
135 | ||
136 | + | |
137 | if (nConstrained > 0){ | |
138 | isConstrained = 1; | |
139 | ||
# | Line 132 | Line 155 | template<typename T> void Integrator<T>::checkConstrai | |
155 | } | |
156 | ||
157 | ||
158 | < | // save oldAtoms to check for lode balanceing later on. |
158 | > | // save oldAtoms to check for lode balancing later on. |
159 | ||
160 | oldAtoms = nAtoms; | |
161 | ||
# | Line 144 | Line 167 | template<typename T> void Integrator<T>::checkConstrai | |
167 | ||
168 | delete[] temp_con; | |
169 | } | |
170 | + | */ |
171 | ||
148 | – | |
172 | template<typename T> void Integrator<T>::integrate(void){ | |
150 | – | int i, j; // loop counters |
173 | ||
174 | double runTime = info->run_time; | |
175 | double sampleTime = info->sampleTime; | |
176 | double statusTime = info->statusTime; | |
177 | double thermalTime = info->thermalTime; | |
178 | + | double resetTime = info->resetTime; |
179 | ||
180 | + | double difference; |
181 | double currSample; | |
182 | double currThermal; | |
183 | double currStatus; | |
184 | + | double currReset; |
185 | ||
186 | int calcPot, calcStress; | |
162 | – | int isError; |
187 | ||
188 | tStats = new Thermo(info); | |
189 | statOut = new StatWriter(info); | |
190 | dumpOut = new DumpWriter(info); | |
191 | ||
192 | atoms = info->atoms; | |
169 | – | DirectionalAtom* dAtom; |
193 | ||
194 | dt = info->dt; | |
195 | dt2 = 0.5 * dt; | |
196 | + | |
197 | + | readyCheck(); |
198 | + | |
199 | + | // remove center of mass drift velocity (in case we passed in a configuration |
200 | + | // that was drifting |
201 | + | tStats->removeCOMdrift(); |
202 | + | //tStats->removeAngularMomentum(); |
203 | + | |
204 | + | // initialize the retraints if necessary |
205 | + | if (info->useSolidThermInt && !info->useLiquidThermInt) { |
206 | + | myFF->initRestraints(); |
207 | + | } |
208 | ||
209 | // initialize the forces before the first step | |
210 | ||
211 | calcForce(1, 1); | |
177 | – | // myFF->doForces(1,1); |
212 | ||
213 | + | //execute constraint algorithm to make sure at the very beginning the system is constrained |
214 | + | //consFramework->doPreConstraint(); |
215 | + | //consFramework->doConstrainA(); |
216 | + | //calcForce(1, 1); |
217 | + | //consFramework->doConstrainB(); |
218 | + | |
219 | if (info->setTemp){ | |
220 | thermalize(); | |
221 | } | |
222 | ||
183 | – | calcPot = 0; |
184 | – | calcStress = 0; |
185 | – | currSample = sampleTime; |
186 | – | currThermal = thermalTime; |
187 | – | currStatus = statusTime; |
188 | – | |
223 | calcPot = 0; | |
224 | calcStress = 0; | |
225 | currSample = sampleTime + info->getTime(); | |
226 | currThermal = thermalTime+ info->getTime(); | |
227 | currStatus = statusTime + info->getTime(); | |
228 | < | >>>>>>> 1.18 |
228 | > | currReset = resetTime + info->getTime(); |
229 | ||
230 | dumpOut->writeDump(info->getTime()); | |
231 | statOut->writeStat(info->getTime()); | |
232 | ||
199 | – | readyCheck(); |
233 | ||
234 | #ifdef IS_MPI | |
235 | strcpy(checkPointMsg, "The integrator is ready to go."); | |
236 | MPIcheckPoint(); | |
237 | #endif // is_mpi | |
238 | ||
239 | < | while (info->getTime() < runTime){ |
240 | < | if ((info->getTime() + dt) >= currStatus){ |
239 | > | while (info->getTime() < runTime && !stopIntegrator()){ |
240 | > | difference = info->getTime() + dt - currStatus; |
241 | > | if (difference > 0 || fabs(difference) < 1e-4 ){ |
242 | calcPot = 1; | |
243 | calcStress = 1; | |
244 | } | |
245 | ||
246 | + | #ifdef PROFILE |
247 | + | startProfile( pro1 ); |
248 | + | #endif |
249 | + | |
250 | integrateStep(calcPot, calcStress); | |
251 | ||
252 | + | #ifdef PROFILE |
253 | + | endProfile( pro1 ); |
254 | + | |
255 | + | startProfile( pro2 ); |
256 | + | #endif // profile |
257 | + | |
258 | info->incrTime(dt); | |
259 | ||
260 | if (info->setTemp){ | |
# | Line 226 | Line 270 | template<typename T> void Integrator<T>::integrate(voi | |
270 | } | |
271 | ||
272 | if (info->getTime() >= currStatus){ | |
273 | < | statOut->writeStat(info->getTime()); |
274 | < | calcPot = 0; |
273 | > | statOut->writeStat(info->getTime()); |
274 | > | calcPot = 0; |
275 | calcStress = 0; | |
276 | currStatus += statusTime; | |
277 | < | } |
277 | > | } |
278 | ||
279 | + | if (info->resetIntegrator){ |
280 | + | if (info->getTime() >= currReset){ |
281 | + | this->resetIntegrator(); |
282 | + | currReset += resetTime; |
283 | + | } |
284 | + | } |
285 | + | |
286 | + | #ifdef PROFILE |
287 | + | endProfile( pro2 ); |
288 | + | #endif //profile |
289 | + | |
290 | #ifdef IS_MPI | |
291 | strcpy(checkPointMsg, "successfully took a time step."); | |
292 | MPIcheckPoint(); | |
293 | #endif // is_mpi | |
294 | } | |
295 | ||
296 | < | dumpOut->writeFinal(info->getTime()); |
296 | > | // dump out a file containing the omega values for the final configuration |
297 | > | if (info->useSolidThermInt && !info->useLiquidThermInt) |
298 | > | myFF->dumpzAngle(); |
299 | > | |
300 | ||
301 | delete dumpOut; | |
302 | delete statOut; | |
# | Line 247 | Line 305 | template<typename T> void Integrator<T>::integrateStep | |
305 | template<typename T> void Integrator<T>::integrateStep(int calcPot, | |
306 | int calcStress){ | |
307 | // Position full step, and velocity half step | |
250 | – | preMove(); |
308 | ||
309 | + | #ifdef PROFILE |
310 | + | startProfile(pro3); |
311 | + | #endif //profile |
312 | + | |
313 | + | //save old state (position, velocity etc) |
314 | + | consFramework->doPreConstraint(); |
315 | + | |
316 | + | #ifdef PROFILE |
317 | + | endProfile(pro3); |
318 | + | |
319 | + | startProfile(pro4); |
320 | + | #endif // profile |
321 | + | |
322 | moveA(); | |
323 | ||
324 | < | if (nConstrained){ |
325 | < | constrainA(); |
326 | < | } |
324 | > | #ifdef PROFILE |
325 | > | endProfile(pro4); |
326 | > | |
327 | > | startProfile(pro5); |
328 | > | #endif//profile |
329 | ||
330 | ||
331 | #ifdef IS_MPI | |
# | Line 261 | Line 333 | template<typename T> void Integrator<T>::integrateStep | |
333 | MPIcheckPoint(); | |
334 | #endif // is_mpi | |
335 | ||
264 | – | |
336 | // calc forces | |
266 | – | |
337 | calcForce(calcPot, calcStress); | |
338 | ||
339 | #ifdef IS_MPI | |
# | Line 271 | Line 341 | template<typename T> void Integrator<T>::integrateStep | |
341 | MPIcheckPoint(); | |
342 | #endif // is_mpi | |
343 | ||
344 | + | #ifdef PROFILE |
345 | + | endProfile( pro5 ); |
346 | ||
347 | + | startProfile( pro6 ); |
348 | + | #endif //profile |
349 | + | |
350 | + | consFramework->doPreConstraint(); |
351 | + | |
352 | // finish the velocity half step | |
353 | ||
354 | moveB(); | |
355 | ||
356 | < | if (nConstrained){ |
357 | < | constrainB(); |
358 | < | } |
356 | > | #ifdef PROFILE |
357 | > | endProfile(pro6); |
358 | > | #endif // profile |
359 | ||
360 | #ifdef IS_MPI | |
361 | strcpy(checkPointMsg, "Succesful moveB\n"); | |
# | Line 288 | Line 365 | template<typename T> void Integrator<T>::moveA(void){ | |
365 | ||
366 | ||
367 | template<typename T> void Integrator<T>::moveA(void){ | |
368 | < | int i, j; |
368 | > | size_t i, j; |
369 | DirectionalAtom* dAtom; | |
370 | double Tb[3], ji[3]; | |
294 | – | double A[3][3], I[3][3]; |
295 | – | double angle; |
371 | double vel[3], pos[3], frc[3]; | |
372 | double mass; | |
373 | + | double omega; |
374 | + | |
375 | + | for (i = 0; i < integrableObjects.size() ; i++){ |
376 | + | integrableObjects[i]->getVel(vel); |
377 | + | integrableObjects[i]->getPos(pos); |
378 | + | integrableObjects[i]->getFrc(frc); |
379 | + | |
380 | + | mass = integrableObjects[i]->getMass(); |
381 | ||
299 | – | for (i = 0; i < nAtoms; i++){ |
300 | – | atoms[i]->getVel(vel); |
301 | – | atoms[i]->getPos(pos); |
302 | – | atoms[i]->getFrc(frc); |
303 | – | |
304 | – | mass = atoms[i]->getMass(); |
305 | – | |
382 | for (j = 0; j < 3; j++){ | |
383 | // velocity half step | |
384 | vel[j] += (dt2 * frc[j] / mass) * eConvert; | |
# | Line 310 | Line 386 | template<typename T> void Integrator<T>::moveA(void){ | |
386 | pos[j] += dt * vel[j]; | |
387 | } | |
388 | ||
389 | < | atoms[i]->setVel(vel); |
390 | < | atoms[i]->setPos(pos); |
389 | > | integrableObjects[i]->setVel(vel); |
390 | > | integrableObjects[i]->setPos(pos); |
391 | ||
392 | < | if (atoms[i]->isDirectional()){ |
317 | < | dAtom = (DirectionalAtom *) atoms[i]; |
392 | > | if (integrableObjects[i]->isDirectional()){ |
393 | ||
394 | // get and convert the torque to body frame | |
395 | ||
396 | < | dAtom->getTrq(Tb); |
397 | < | dAtom->lab2Body(Tb); |
396 | > | integrableObjects[i]->getTrq(Tb); |
397 | > | integrableObjects[i]->lab2Body(Tb); |
398 | ||
399 | // get the angular momentum, and propagate a half step | |
400 | ||
401 | < | dAtom->getJ(ji); |
401 | > | integrableObjects[i]->getJ(ji); |
402 | ||
403 | for (j = 0; j < 3; j++) | |
404 | ji[j] += (dt2 * Tb[j]) * eConvert; | |
405 | ||
406 | < | // use the angular velocities to propagate the rotation matrix a |
332 | < | // full time step |
406 | > | this->rotationPropagation( integrableObjects[i], ji ); |
407 | ||
408 | < | dAtom->getA(A); |
335 | < | dAtom->getI(I); |
336 | < | |
337 | < | // rotate about the x-axis |
338 | < | angle = dt2 * ji[0] / I[0][0]; |
339 | < | this->rotate(1, 2, angle, ji, A); |
340 | < | |
341 | < | // rotate about the y-axis |
342 | < | angle = dt2 * ji[1] / I[1][1]; |
343 | < | this->rotate(2, 0, angle, ji, A); |
344 | < | |
345 | < | // rotate about the z-axis |
346 | < | angle = dt * ji[2] / I[2][2]; |
347 | < | this->rotate(0, 1, angle, ji, A); |
348 | < | |
349 | < | // rotate about the y-axis |
350 | < | angle = dt2 * ji[1] / I[1][1]; |
351 | < | this->rotate(2, 0, angle, ji, A); |
408 | > | integrableObjects[i]->setJ(ji); |
409 | ||
353 | – | // rotate about the x-axis |
354 | – | angle = dt2 * ji[0] / I[0][0]; |
355 | – | this->rotate(1, 2, angle, ji, A); |
356 | – | |
357 | – | |
358 | – | dAtom->setJ(ji); |
359 | – | dAtom->setA(A); |
410 | } | |
411 | } | |
412 | + | |
413 | + | consFramework->doConstrainA(); |
414 | } | |
415 | ||
416 | ||
417 | template<typename T> void Integrator<T>::moveB(void){ | |
418 | int i, j; | |
367 | – | DirectionalAtom* dAtom; |
419 | double Tb[3], ji[3]; | |
420 | double vel[3], frc[3]; | |
421 | double mass; | |
422 | ||
423 | < | for (i = 0; i < nAtoms; i++){ |
424 | < | atoms[i]->getVel(vel); |
425 | < | atoms[i]->getFrc(frc); |
423 | > | for (i = 0; i < integrableObjects.size(); i++){ |
424 | > | integrableObjects[i]->getVel(vel); |
425 | > | integrableObjects[i]->getFrc(frc); |
426 | ||
427 | < | mass = atoms[i]->getMass(); |
427 | > | mass = integrableObjects[i]->getMass(); |
428 | ||
429 | // velocity half step | |
430 | for (j = 0; j < 3; j++) | |
431 | vel[j] += (dt2 * frc[j] / mass) * eConvert; | |
432 | ||
433 | < | atoms[i]->setVel(vel); |
433 | > | integrableObjects[i]->setVel(vel); |
434 | ||
435 | < | if (atoms[i]->isDirectional()){ |
385 | < | dAtom = (DirectionalAtom *) atoms[i]; |
386 | < | |
387 | < | // get and convert the torque to body frame |
435 | > | if (integrableObjects[i]->isDirectional()){ |
436 | ||
437 | < | dAtom->getTrq(Tb); |
390 | < | dAtom->lab2Body(Tb); |
437 | > | // get and convert the torque to body frame |
438 | ||
439 | + | integrableObjects[i]->getTrq(Tb); |
440 | + | integrableObjects[i]->lab2Body(Tb); |
441 | + | |
442 | // get the angular momentum, and propagate a half step | |
443 | ||
444 | < | dAtom->getJ(ji); |
444 | > | integrableObjects[i]->getJ(ji); |
445 | ||
446 | for (j = 0; j < 3; j++) | |
447 | ji[j] += (dt2 * Tb[j]) * eConvert; | |
448 | ||
449 | ||
450 | < | dAtom->setJ(ji); |
450 | > | integrableObjects[i]->setJ(ji); |
451 | } | |
452 | + | |
453 | } | |
454 | + | |
455 | + | consFramework->doConstrainB(); |
456 | } | |
457 | ||
458 | + | /* |
459 | template<typename T> void Integrator<T>::preMove(void){ | |
460 | int i, j; | |
461 | double pos[3]; | |
# | Line 418 | Line 472 | template<typename T> void Integrator<T>::constrainA(){ | |
472 | } | |
473 | ||
474 | template<typename T> void Integrator<T>::constrainA(){ | |
475 | < | int i, j, k; |
475 | > | int i, j; |
476 | int done; | |
477 | double posA[3], posB[3]; | |
478 | double velA[3], velB[3]; | |
# | Line 558 | Line 612 | template<typename T> void Integrator<T>::constrainA(){ | |
612 | painCave.isFatal = 1; | |
613 | simError(); | |
614 | } | |
615 | + | |
616 | } | |
617 | ||
618 | template<typename T> void Integrator<T>::constrainB(void){ | |
619 | < | int i, j, k; |
619 | > | int i, j; |
620 | int done; | |
621 | double posA[3], posB[3]; | |
622 | double velA[3], velB[3]; | |
# | Line 570 | Line 625 | template<typename T> void Integrator<T>::constrainB(vo | |
625 | int a, b, ax, ay, az, bx, by, bz; | |
626 | double rma, rmb; | |
627 | double dx, dy, dz; | |
628 | < | double rabsq, pabsq, rvab; |
574 | < | double diffsq; |
628 | > | double rvab; |
629 | double gab; | |
630 | int iteration; | |
631 | ||
# | Line 658 | Line 712 | template<typename T> void Integrator<T>::constrainB(vo | |
712 | iteration); | |
713 | painCave.isFatal = 1; | |
714 | simError(); | |
715 | + | } |
716 | + | } |
717 | + | */ |
718 | + | template<typename T> void Integrator<T>::rotationPropagation |
719 | + | ( StuntDouble* sd, double ji[3] ){ |
720 | + | |
721 | + | double angle; |
722 | + | double A[3][3], I[3][3]; |
723 | + | int i, j, k; |
724 | + | |
725 | + | // use the angular velocities to propagate the rotation matrix a |
726 | + | // full time step |
727 | + | |
728 | + | sd->getA(A); |
729 | + | sd->getI(I); |
730 | + | |
731 | + | if (sd->isLinear()) { |
732 | + | i = sd->linearAxis(); |
733 | + | j = (i+1)%3; |
734 | + | k = (i+2)%3; |
735 | + | |
736 | + | angle = dt2 * ji[j] / I[j][j]; |
737 | + | this->rotate( k, i, angle, ji, A ); |
738 | + | |
739 | + | angle = dt * ji[k] / I[k][k]; |
740 | + | this->rotate( i, j, angle, ji, A); |
741 | + | |
742 | + | angle = dt2 * ji[j] / I[j][j]; |
743 | + | this->rotate( k, i, angle, ji, A ); |
744 | + | |
745 | + | } else { |
746 | + | // rotate about the x-axis |
747 | + | angle = dt2 * ji[0] / I[0][0]; |
748 | + | this->rotate( 1, 2, angle, ji, A ); |
749 | + | |
750 | + | // rotate about the y-axis |
751 | + | angle = dt2 * ji[1] / I[1][1]; |
752 | + | this->rotate( 2, 0, angle, ji, A ); |
753 | + | |
754 | + | // rotate about the z-axis |
755 | + | angle = dt * ji[2] / I[2][2]; |
756 | + | sd->addZangle(angle); |
757 | + | this->rotate( 0, 1, angle, ji, A); |
758 | + | |
759 | + | // rotate about the y-axis |
760 | + | angle = dt2 * ji[1] / I[1][1]; |
761 | + | this->rotate( 2, 0, angle, ji, A ); |
762 | + | |
763 | + | // rotate about the x-axis |
764 | + | angle = dt2 * ji[0] / I[0][0]; |
765 | + | this->rotate( 1, 2, angle, ji, A ); |
766 | + | |
767 | } | |
768 | + | sd->setA( A ); |
769 | } | |
770 | ||
771 | template<typename T> void Integrator<T>::rotate(int axes1, int axes2, | |
# | Line 726 | Line 833 | template<typename T> void Integrator<T>::rotate(int ax | |
833 | } | |
834 | } | |
835 | ||
836 | < | // rotate the Rotation matrix acording to: |
836 | > | // rotate the Rotation matrix acording to: |
837 | // A[][] = A[][] * transpose(rot[][]) | |
838 | ||
839 | ||
# | Line 751 | Line 858 | template<typename T> void Integrator<T>::thermalize(){ | |
858 | template<typename T> void Integrator<T>::thermalize(){ | |
859 | tStats->velocitize(); | |
860 | } | |
861 | + | |
862 | + | template<typename T> double Integrator<T>::getConservedQuantity(void){ |
863 | + | return tStats->getTotalE(); |
864 | + | } |
865 | + | template<typename T> string Integrator<T>::getAdditionalParameters(void){ |
866 | + | //By default, return a null string |
867 | + | //The reason we use string instead of char* is that if we use char*, we will |
868 | + | //return a pointer point to local variable which might cause problem |
869 | + | return string(); |
870 | + | } |
871 | + | |
872 | + | |
873 | + | template<typename T> void Integrator<T>::printQuaternion(StuntDouble* sd){ |
874 | + | Mat4x4d S; |
875 | + | double I[3][3]; |
876 | + | Vector4d j4; |
877 | + | Vector3d j; |
878 | + | Vector3d tempJ; |
879 | + | Vector4d qdot; |
880 | + | Vector4d omega4; |
881 | + | Mat4x4d I4; |
882 | + | Quaternion q; |
883 | + | double I0; |
884 | + | Vector4d p_qua; |
885 | + | |
886 | + | if (sd->isDirectional()){ |
887 | + | sd->getQ(q.vec); |
888 | + | sd->getI(I); |
889 | + | sd->getJ(j.vec); |
890 | + | |
891 | + | //omega4[0] = 0.0; |
892 | + | //omega4[1] = j[0]/I[0][0]; |
893 | + | //omega4[2] = j[1]/I[1][1]; |
894 | + | //omega4[3] = j[2]/I[2][2]; |
895 | + | |
896 | + | //S = getS(q); |
897 | + | //qdot = 0.5 * S * omega4; |
898 | + | |
899 | + | //I0 = (qdot[1] * q[1] * I[0][0] + qdot[2] * q[2] * I[1][1] + qdot[3] * q[3] * I[2][2])/(qdot[1] * q[1]+ qdot[2] * q[2] + qdot[3] * q[3]); |
900 | + | |
901 | + | //I4.element[0][0] = I0; |
902 | + | //I4.element[1][1] = I[0][0]; |
903 | + | //I4.element[2][2] = I[1][1]; |
904 | + | //I4.element[3][3] = I[2][2]; |
905 | + | |
906 | + | S = getS(q); |
907 | + | j4[0] = 0.0; |
908 | + | j4[1] = j[0]; |
909 | + | j4[2] = j[1]; |
910 | + | j4[3] = j[2]; |
911 | + | |
912 | + | p_qua = 2 * S * j4; |
913 | + | |
914 | + | j4 = 0.5 * S.transpose() * p_qua; |
915 | + | //cout << "q0^2 + q1^2 + q2^2 + q3^2 = " << q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3] << endl; |
916 | + | //cout << "q0*q0dot + q1*q1dot + q2 *q2dot + q3*q3dot = " <<q[0]*qdot[0] + q[1]*qdot[1] + q[2]*qdot[2] + q[3]*qdot[3] << endl; |
917 | + | //cout << "q1*q1dot* Ixx + q2*q2dot* Iyy + q3 *q3dot* Izz = " << qdot[1] * q[1] * I[0][0] + qdot[2] * q[2] * I[1][1] + qdot[3] * q[3] * I[2][2] << endl; |
918 | + | //cout << "q1*q1dot + q2 *q2dot + q3*q3dot = " << qdot[1] * q[1]+ qdot[2] * q[2] + qdot[3] * q[3] << endl; |
919 | + | //cout << "I0 = " << I0 << endl; |
920 | + | cout << "p_qua[0] = " << p_qua[0] << endl; |
921 | + | } |
922 | + | } |
923 | + | |
924 | + | template<typename T> Mat4x4d Integrator<T>::getS(const Quaternion& q){ |
925 | + | Mat4x4d result; |
926 | + | |
927 | + | result.element[0][0] = q.x; |
928 | + | result.element[0][1] = -q.y; |
929 | + | result.element[0][2] = -q.z; |
930 | + | result.element[0][3] = -q.w; |
931 | + | |
932 | + | result.element[1][0] = q.y; |
933 | + | result.element[1][1] = q.x; |
934 | + | result.element[1][2] = -q.w; |
935 | + | result.element[1][3] = q.z; |
936 | + | |
937 | + | result.element[2][0] = q.z; |
938 | + | result.element[2][1] = q.w; |
939 | + | result.element[2][2] = q.x; |
940 | + | result.element[2][3] = -q.y; |
941 | + | |
942 | + | result.element[3][0] = q.w; |
943 | + | result.element[3][1] = -q.z; |
944 | + | result.element[3][2] = q.y; |
945 | + | result.element[3][3] = q.x; |
946 | + | |
947 | + | return result; |
948 | + | } |
949 | + |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |