179 |
|
|
180 |
|
readyCheck(); |
181 |
|
|
182 |
+ |
// remove center of mass drift velocity (in case we passed in a configuration |
183 |
+ |
// that was drifting |
184 |
+ |
tStats->removeCOMdrift(); |
185 |
+ |
|
186 |
|
// initialize the forces before the first step |
187 |
|
|
188 |
|
calcForce(1, 1); |
185 |
– |
|
186 |
– |
//temp test |
187 |
– |
tStats->getPotential(); |
189 |
|
|
190 |
|
if (nConstrained){ |
191 |
|
preMove(); |
214 |
|
MPIcheckPoint(); |
215 |
|
#endif // is_mpi |
216 |
|
|
217 |
< |
while (info->getTime() < runTime){ |
217 |
> |
while (info->getTime() < runTime && !stopIntegrator()){ |
218 |
|
if ((info->getTime() + dt) >= currStatus){ |
219 |
|
calcPot = 1; |
220 |
|
calcStress = 1; |
346 |
|
integrableObjects[i]->getVel(vel); |
347 |
|
integrableObjects[i]->getPos(pos); |
348 |
|
integrableObjects[i]->getFrc(frc); |
349 |
+ |
|
350 |
+ |
std::cerr << "i =\t" << i << "\t" << frc[0] << "\t" << frc[1]<< "\t" << frc[2] << "\n"; |
351 |
|
|
352 |
|
mass = integrableObjects[i]->getMass(); |
353 |
|
|
693 |
|
|
694 |
|
double angle; |
695 |
|
double A[3][3], I[3][3]; |
696 |
+ |
int i, j, k; |
697 |
|
|
698 |
|
// use the angular velocities to propagate the rotation matrix a |
699 |
|
// full time step |
701 |
|
sd->getA(A); |
702 |
|
sd->getI(I); |
703 |
|
|
704 |
< |
// rotate about the x-axis |
705 |
< |
angle = dt2 * ji[0] / I[0][0]; |
706 |
< |
this->rotate( 1, 2, angle, ji, A ); |
704 |
> |
if (sd->isLinear()) { |
705 |
> |
i = sd->linearAxis(); |
706 |
> |
j = (i+1)%3; |
707 |
> |
k = (i+2)%3; |
708 |
> |
|
709 |
> |
angle = dt2 * ji[j] / I[j][j]; |
710 |
> |
this->rotate( k, i, angle, ji, A ); |
711 |
|
|
712 |
< |
// rotate about the y-axis |
713 |
< |
angle = dt2 * ji[1] / I[1][1]; |
706 |
< |
this->rotate( 2, 0, angle, ji, A ); |
712 |
> |
angle = dt * ji[k] / I[k][k]; |
713 |
> |
this->rotate( i, j, angle, ji, A); |
714 |
|
|
715 |
< |
// rotate about the z-axis |
716 |
< |
angle = dt * ji[2] / I[2][2]; |
710 |
< |
this->rotate( 0, 1, angle, ji, A); |
715 |
> |
angle = dt2 * ji[j] / I[j][j]; |
716 |
> |
this->rotate( k, i, angle, ji, A ); |
717 |
|
|
718 |
< |
// rotate about the y-axis |
719 |
< |
angle = dt2 * ji[1] / I[1][1]; |
720 |
< |
this->rotate( 2, 0, angle, ji, A ); |
721 |
< |
|
722 |
< |
// rotate about the x-axis |
723 |
< |
angle = dt2 * ji[0] / I[0][0]; |
724 |
< |
this->rotate( 1, 2, angle, ji, A ); |
725 |
< |
|
718 |
> |
} else { |
719 |
> |
// rotate about the x-axis |
720 |
> |
angle = dt2 * ji[0] / I[0][0]; |
721 |
> |
this->rotate( 1, 2, angle, ji, A ); |
722 |
> |
|
723 |
> |
// rotate about the y-axis |
724 |
> |
angle = dt2 * ji[1] / I[1][1]; |
725 |
> |
this->rotate( 2, 0, angle, ji, A ); |
726 |
> |
|
727 |
> |
// rotate about the z-axis |
728 |
> |
angle = dt * ji[2] / I[2][2]; |
729 |
> |
this->rotate( 0, 1, angle, ji, A); |
730 |
> |
|
731 |
> |
// rotate about the y-axis |
732 |
> |
angle = dt2 * ji[1] / I[1][1]; |
733 |
> |
this->rotate( 2, 0, angle, ji, A ); |
734 |
> |
|
735 |
> |
// rotate about the x-axis |
736 |
> |
angle = dt2 * ji[0] / I[0][0]; |
737 |
> |
this->rotate( 1, 2, angle, ji, A ); |
738 |
> |
|
739 |
> |
} |
740 |
|
sd->setA( A ); |
741 |
|
} |
742 |
|
|