@@ -21,7 +21,7 @@ PositionBasedElasticRodsTSC::~PositionBasedElasticRodsTSC()
21
21
void PositionBasedElasticRodsTSC::step (SimulationModel &model)
22
22
{
23
23
TimeManager *tm = TimeManager::getCurrent ();
24
- const Real h = tm ->getTimeStepSize ();
24
+ const Real hOld = tm ->getTimeStepSize ();
25
25
PositionBasedElasticRodsModel &ermodel = (PositionBasedElasticRodsModel&)model;
26
26
27
27
// ////////////////////////////////////////////////////////////////////////
@@ -33,86 +33,102 @@ void PositionBasedElasticRodsTSC::step(SimulationModel &model)
33
33
ParticleData &pg = ermodel.getGhostParticles ();
34
34
35
35
const int numBodies = (int )rb.size ();
36
- #pragma omp parallel if(numBodies > MIN_PARALLEL_SIZE) default(shared)
37
- {
38
- #pragma omp for schedule(static) nowait
39
- for (int i = 0 ; i < numBodies; i++)
40
- {
41
- rb[i]->getLastPosition () = rb[i]->getOldPosition ();
42
- rb[i]->getOldPosition () = rb[i]->getPosition ();
43
- TimeIntegration::semiImplicitEuler (h, rb[i]->getMass (), rb[i]->getPosition (), rb[i]->getVelocity (), rb[i]->getAcceleration ());
44
- rb[i]->getLastRotation () = rb[i]->getOldRotation ();
45
- rb[i]->getOldRotation () = rb[i]->getRotation ();
46
- TimeIntegration::semiImplicitEulerRotation (h, rb[i]->getMass (), rb[i]->getInertiaTensorInverseW (), rb[i]->getRotation (), rb[i]->getAngularVelocity (), rb[i]->getTorque ());
47
- rb[i]->rotationUpdated ();
48
- }
49
36
50
- // ////////////////////////////////////////////////////////////////////////
51
- // particle model
52
- // ////////////////////////////////////////////////////////////////////////
53
- # pragma omp for schedule(static)
54
- for ( int i = 0 ; i < ( int ) pd. size (); i++ )
37
+ Real h = hOld / (Real)m_subSteps;
38
+ tm -> setTimeStepSize (h);
39
+ for ( unsigned int step = 0 ; step < m_subSteps; step++)
40
+ {
41
+ # pragma omp parallel if(numBodies > MIN_PARALLEL_SIZE) default(shared )
55
42
{
56
- pd.getLastPosition (i) = pd.getOldPosition (i);
57
- pd.getOldPosition (i) = pd.getPosition (i);
58
- pd.getVelocity (i) *= (1 .0f - m_damping);
43
+ #pragma omp for schedule(static) nowait
44
+ for (int i = 0 ; i < numBodies; i++)
45
+ {
46
+ rb[i]->getLastPosition () = rb[i]->getOldPosition ();
47
+ rb[i]->getOldPosition () = rb[i]->getPosition ();
48
+ TimeIntegration::semiImplicitEuler (h, rb[i]->getMass (), rb[i]->getPosition (), rb[i]->getVelocity (), rb[i]->getAcceleration ());
49
+ rb[i]->getLastRotation () = rb[i]->getOldRotation ();
50
+ rb[i]->getOldRotation () = rb[i]->getRotation ();
51
+ TimeIntegration::semiImplicitEulerRotation (h, rb[i]->getMass (), rb[i]->getInertiaTensorInverseW (), rb[i]->getRotation (), rb[i]->getAngularVelocity (), rb[i]->getTorque ());
52
+ rb[i]->rotationUpdated ();
53
+ }
54
+
55
+ // ////////////////////////////////////////////////////////////////////////
56
+ // particle model
57
+ // ////////////////////////////////////////////////////////////////////////
58
+ #pragma omp for schedule(static)
59
+ for (int i = 0 ; i < (int ) pd.size (); i++)
60
+ {
61
+ pd.getLastPosition (i) = pd.getOldPosition (i);
62
+ pd.getOldPosition (i) = pd.getPosition (i);
63
+ pd.getVelocity (i) *= (1 .0f - m_damping);
59
64
60
- TimeIntegration::semiImplicitEuler (h, pd.getMass (i), pd.getPosition (i), pd.getVelocity (i), pd.getAcceleration (i));
61
- }
65
+ TimeIntegration::semiImplicitEuler (h, pd.getMass (i), pd.getPosition (i), pd.getVelocity (i), pd.getAcceleration (i));
66
+ }
62
67
63
- for (int i = 0 ; i < (int )pg.size (); i++)
64
- {
65
- pg.getLastPosition (i) = pg.getOldPosition (i);
66
- pg.getOldPosition (i) = pg.getPosition (i);
68
+ for (int i = 0 ; i < (int )pg.size (); i++)
69
+ {
70
+ pg.getLastPosition (i) = pg.getOldPosition (i);
71
+ pg.getOldPosition (i) = pg.getPosition (i);
67
72
68
- pg.getVelocity (i) *= (1 .0f - m_damping);
73
+ pg.getVelocity (i) *= (1 .0f - m_damping);
69
74
70
- TimeIntegration::semiImplicitEuler (h, pg.getMass (i), pg.getPosition (i), pg.getVelocity (i), pg.getAcceleration (i));
75
+ TimeIntegration::semiImplicitEuler (h, pg.getMass (i), pg.getPosition (i), pg.getVelocity (i), pg.getAcceleration (i));
76
+ }
71
77
}
72
- }
73
78
74
- positionConstraintProjection (model);
79
+ positionConstraintProjection (model);
75
80
76
- #pragma omp parallel if(numBodies > MIN_PARALLEL_SIZE) default(shared)
77
- {
78
- // Update velocities
79
- #pragma omp for schedule(static) nowait
80
- for (int i = 0 ; i < numBodies; i++)
81
- {
82
- if (m_velocityUpdateMethod == 0 )
81
+ #pragma omp parallel if(numBodies > MIN_PARALLEL_SIZE) default(shared)
82
+ {
83
+ // Update velocities
84
+ #pragma omp for schedule(static) nowait
85
+ for (int i = 0 ; i < numBodies; i++)
86
+ {
87
+ if (m_velocityUpdateMethod == 0 )
88
+ {
89
+ TimeIntegration::velocityUpdateFirstOrder (h, rb[i]->getMass (), rb[i]->getPosition (), rb[i]->getOldPosition (), rb[i]->getVelocity ());
90
+ TimeIntegration::angularVelocityUpdateFirstOrder (h, rb[i]->getMass (), rb[i]->getRotation (), rb[i]->getOldRotation (), rb[i]->getAngularVelocity ());
91
+ }
92
+ else
93
+ {
94
+ TimeIntegration::velocityUpdateSecondOrder (h, rb[i]->getMass (), rb[i]->getPosition (), rb[i]->getOldPosition (), rb[i]->getLastPosition (), rb[i]->getVelocity ());
95
+ TimeIntegration::angularVelocityUpdateSecondOrder (h, rb[i]->getMass (), rb[i]->getRotation (), rb[i]->getOldRotation (), rb[i]->getLastRotation (), rb[i]->getAngularVelocity ());
96
+ }
97
+ }
98
+
99
+ // Update velocities
100
+ #pragma omp for schedule(static)
101
+ for (int i = 0 ; i < (int ) pd.size (); i++)
83
102
{
84
- TimeIntegration::velocityUpdateFirstOrder (h, rb[i]->getMass (), rb[i]->getPosition (), rb[i]->getOldPosition (), rb[i]->getVelocity ());
85
- TimeIntegration::angularVelocityUpdateFirstOrder (h, rb[i]->getMass (), rb[i]->getRotation (), rb[i]->getOldRotation (), rb[i]->getAngularVelocity ());
103
+ if (m_velocityUpdateMethod == 0 )
104
+ TimeIntegration::velocityUpdateFirstOrder (h, pd.getMass (i), pd.getPosition (i), pd.getOldPosition (i), pd.getVelocity (i));
105
+ else
106
+ TimeIntegration::velocityUpdateSecondOrder (h, pd.getMass (i), pd.getPosition (i), pd.getOldPosition (i), pd.getLastPosition (i), pd.getVelocity (i));
86
107
}
87
- else
108
+
109
+ #pragma omp for schedule(static)
110
+ for (int i = 0 ; i < (int )pg.size (); i++)
88
111
{
89
- TimeIntegration::velocityUpdateSecondOrder (h, rb[i]->getMass (), rb[i]->getPosition (), rb[i]->getOldPosition (), rb[i]->getLastPosition (), rb[i]->getVelocity ());
90
- TimeIntegration::angularVelocityUpdateSecondOrder (h, rb[i]->getMass (), rb[i]->getRotation (), rb[i]->getOldRotation (), rb[i]->getLastRotation (), rb[i]->getAngularVelocity ());
112
+ if (m_velocityUpdateMethod == 0 )
113
+ TimeIntegration::velocityUpdateFirstOrder (h, pg.getMass (i), pg.getPosition (i), pg.getOldPosition (i), pg.getVelocity (i));
114
+ else
115
+ TimeIntegration::velocityUpdateSecondOrder (h, pg.getMass (i), pg.getPosition (i), pg.getOldPosition (i), pg.getLastPosition (i), pg.getVelocity (i));
91
116
}
92
- // update geometry
93
- rb[i]->getGeometry ().updateMeshTransformation (rb[i]->getPosition (), rb[i]->getRotationMatrix ());
94
- }
117
+ }
118
+ }
119
+ h = hOld;
120
+ tm ->setTimeStepSize (hOld);
95
121
122
+ #pragma omp parallel if(numBodies > MIN_PARALLEL_SIZE) default(shared)
123
+ {
96
124
// Update velocities
97
- #pragma omp for schedule(static)
98
- for (int i = 0 ; i < (int ) pd.size (); i++)
99
- {
100
- if (m_velocityUpdateMethod == 0 )
101
- TimeIntegration::velocityUpdateFirstOrder (h, pd.getMass (i), pd.getPosition (i), pd.getOldPosition (i), pd.getVelocity (i));
102
- else
103
- TimeIntegration::velocityUpdateSecondOrder (h, pd.getMass (i), pd.getPosition (i), pd.getOldPosition (i), pd.getLastPosition (i), pd.getVelocity (i));
104
- }
105
-
106
- #pragma omp for schedule(static)
107
- for (int i = 0 ; i < (int )pg.size (); i++)
125
+ #pragma omp for schedule(static) nowait
126
+ for (int i = 0 ; i < numBodies; i++)
108
127
{
109
- if (m_velocityUpdateMethod == 0 )
110
- TimeIntegration::velocityUpdateFirstOrder (h, pg.getMass (i), pg.getPosition (i), pg.getOldPosition (i), pg.getVelocity (i));
111
- else
112
- TimeIntegration::velocityUpdateSecondOrder (h, pg.getMass (i), pg.getPosition (i), pg.getOldPosition (i), pg.getLastPosition (i), pg.getVelocity (i));
128
+ if (rb[i]->getMass () != 0.0 )
129
+ rb[i]->getGeometry ().updateMeshTransformation (rb[i]->getPosition (), rb[i]->getRotationMatrix ());
113
130
}
114
-
115
- }
131
+ }
116
132
117
133
if (m_collisionDetection)
118
134
m_collisionDetection->collisionDetection (model);
0 commit comments