Skip to content

Commit 4fad902

Browse files
author
bender
committed
- added second order velocity update
1 parent 1b0a564 commit 4fad902

File tree

12 files changed

+198
-74
lines changed

12 files changed

+198
-74
lines changed

‎Demos/BarDemo/TetModel.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -60,8 +60,9 @@ void TetModel::reset()
6060
for (unsigned int i = 0; i < nPoints; i++)
6161
{
6262
const Eigen::Vector3f& x0 = pd.getPosition0(i);
63-
pd.getPosition(i) = x0;
63+
pd.getPosition(i) = x0;
6464
pd.getLastPosition(i) = pd.getPosition(i);
65+
pd.getOldPosition(i) = pd.getPosition(i);
6566
pd.getVelocity(i).setZero();
6667
pd.getAcceleration(i).setZero();
6768
}

‎Demos/BarDemo/TimeStepTetModel.cpp‎

Lines changed: 20 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
#include "TimeStepTetModel.h"
22
#include "Demos/Utils/TimeManager.h"
33
#include "PositionBasedDynamics/PositionBasedDynamics.h"
4+
#include "PositionBasedDynamics/TimeIntegration.h"
45

56
using namespace PBD;
67
using namespace std;
@@ -18,17 +19,25 @@ void TimeStepTetModel::step(TetModel &model)
1819
{
1920
TimeManager *tm = TimeManager::getCurrent ();
2021
const float h = tm->getTimeStepSize();
22+
ParticleData &pd = model.getParticleMesh().getVertexData();
2123

2224
clearAccelerations(model);
23-
semiImplicitEulerStep(model, h);
25+
for (unsigned int i = 0; i < pd.size(); i++)
26+
{
27+
pd.getLastPosition(i) = pd.getOldPosition(i);
28+
pd.getOldPosition(i) = pd.getPosition(i);
29+
TimeIntegration::semiImplicitEuler(h, pd.getMass(i), pd.getPosition(i), pd.getVelocity(i), pd.getAcceleration(i));
30+
}
2431

2532
constraintProjection(model);
2633

27-
// Update velocities
28-
ParticleData &pd = model.getParticleMesh().getVertexData();
34+
// Update velocities
2935
for (unsigned int i = 0; i < pd.size(); i++)
3036
{
31-
pd.getVelocity(i) = 1.0f / h * (pd.getPosition(i) - pd.getLastPosition(i));
37+
if (m_velocityUpdateMethod == 0)
38+
TimeIntegration::velocityUpdateFirstOrder(h, pd.getMass(i), pd.getPosition(i), pd.getOldPosition(i), pd.getVelocity(i));
39+
else
40+
TimeIntegration::velocityUpdateSecondOrder(h, pd.getMass(i), pd.getPosition(i), pd.getOldPosition(i), pd.getLastPosition(i), pd.getVelocity(i));
3241
}
3342

3443
// compute new time
@@ -53,27 +62,6 @@ void TimeStepTetModel::clearAccelerations(TetModel &model)
5362
}
5463
}
5564

56-
void TimeStepTetModel::semiImplicitEulerStep(TetModel &model, const float h)
57-
{
58-
ParticleData &pd = model.getParticleMesh().getVertexData();
59-
60-
// h * x'
61-
for (unsigned int i=0; i < pd.size(); i++)
62-
{
63-
if (pd.getMass(i) != 0.0)
64-
{
65-
Eigen::Vector3f &pos = pd.getPosition(i);
66-
Eigen::Vector3f &vel = pd.getVelocity(i);
67-
const Eigen::Vector3f &accel = pd.getAcceleration(i);
68-
Eigen::Vector3f &lastPos = pd.getLastPosition(i);
69-
lastPos = pos;
70-
71-
vel += accel * h;
72-
pos += vel * h;
73-
}
74-
}
75-
}
76-
7765
void TimeStepTetModel::reset(TetModel &model)
7866
{
7967

@@ -114,13 +102,13 @@ void TimeStepTetModel::constraintProjection(TetModel &model)
114102
float restLength = (x2_0 - x1_0).norm();
115103

116104
Eigen::Vector3f corr1, corr2;
117-
const bool res = PositionBasedDynamics::solveDistanceConstraint(x1, invMass1, x2, invMass2, restLength, model.getStiffness(), model.getStiffness(), corr1, corr2);
105+
const bool res = PositionBasedDynamics::solveDistanceConstraint(x1, invMass1, x2, invMass2, restLength, 0.0f, model.getStiffness(), corr1, corr2);
118106

119107
if (res)
120108
{
121-
if (invMass1 != 0)
109+
if (invMass1 != 0.0f)
122110
x1 += corr1;
123-
if (invMass2 != 0)
111+
if (invMass2 != 0.0f)
124112
x2 += corr2;
125113
}
126114
}
@@ -219,13 +207,13 @@ void TimeStepTetModel::constraintProjection(TetModel &model)
219207

220208
if (res)
221209
{
222-
if (invMass1 != 0)
210+
if (invMass1 != 0.0f)
223211
x1 += corr1;
224-
if (invMass2 != 0)
212+
if (invMass2 != 0.0f)
225213
x2 += corr2;
226-
if (invMass3 != 0)
214+
if (invMass3 != 0.0f)
227215
x3 += corr3;
228-
if (invMass4 != 0)
216+
if (invMass4 != 0.0f)
229217
x4 += corr4;
230218
}
231219
}

‎Demos/BarDemo/TimeStepTetModel.h‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,10 +9,10 @@ namespace PBD
99
class TimeStepTetModel
1010
{
1111
protected:
12+
unsigned int m_velocityUpdateMethod;
1213
unsigned int m_simulationMethod;
1314

1415
void clearAccelerations(TetModel &model);
15-
void semiImplicitEulerStep(TetModel &model, const float h);
1616
void constraintProjection(TetModel &model);
1717

1818
public:
@@ -22,6 +22,8 @@ namespace PBD
2222
void step(TetModel &model);
2323
void reset(TetModel &model);
2424

25+
unsigned int getVelocityUpdateMethod() const { return m_velocityUpdateMethod; }
26+
void setVelocityUpdateMethod(unsigned int val) { m_velocityUpdateMethod = val; }
2527
unsigned int getSimulationMethod() const { return m_simulationMethod; }
2628
void setSimulationMethod(unsigned int val) { m_simulationMethod = val; }
2729
};

‎Demos/BarDemo/main.cpp‎

Lines changed: 17 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,8 @@ void TW_CALL setNormalizeShear(const void *value, void *clientData);
3636
void TW_CALL getNormalizeShear(void *value, void *clientData);
3737
void TW_CALL setSimulationMethod(const void *value, void *clientData);
3838
void TW_CALL getSimulationMethod(void *value, void *clientData);
39+
void TW_CALL setVelocityUpdateMethod(const void *value, void *clientData);
40+
void TW_CALL getVelocityUpdateMethod(void *value, void *clientData);
3941

4042
TetModel model;
4143
TimeStepTetModel simulation;
@@ -67,8 +69,10 @@ int main( int argc, char **argv )
6769

6870
TwAddVarRW(MiniGL::getTweakBar(), "Pause", TW_TYPE_BOOLCPP, &pause, " label='Pause' group=Simulation key=SPACE ");
6971
TwAddVarCB(MiniGL::getTweakBar(), "TimeStepSize", TW_TYPE_FLOAT, setTimeStep, getTimeStep, &model, " label='Time step size' min=0.0 max = 0.1 step=0.001 precision=4 group=Simulation ");
70-
TwType enumType = TwDefineEnum("SimulationMethodType", NULL, 0);
71-
TwAddVarCB(MiniGL::getTweakBar(), "SimulationMethod", enumType, setSimulationMethod, getSimulationMethod, &simulation,
72+
TwType enumType = TwDefineEnum("VelocityUpdateMethodType", NULL, 0);
73+
TwAddVarCB(MiniGL::getTweakBar(), "VelocityUpdateMethod", enumType, setVelocityUpdateMethod, getVelocityUpdateMethod, &simulation, " label='Velocity update method' enum='0 {First Order Update}, 1 {Second Order Update}' group=Simulation");
74+
TwType enumType2 = TwDefineEnum("SimulationMethodType", NULL, 0);
75+
TwAddVarCB(MiniGL::getTweakBar(), "SimulationMethod", enumType2, setSimulationMethod, getSimulationMethod, &simulation,
7276
" label='Simulation method' enum='0 {None}, 1 {Volume constraints}, 2 {FEM based PBD}, 3 {Strain based dynamics (no inversion handling)}, 4 {Shape matching (no inversion handling)}' group=Simulation");
7377
TwAddVarCB(MiniGL::getTweakBar(), "Stiffness", TW_TYPE_FLOAT, setStiffness, getStiffness, &model, " label='Stiffness' min=0.0 step=0.1 precision=4 group='Simulation' ");
7478
TwAddVarCB(MiniGL::getTweakBar(), "PoissonRatio", TW_TYPE_FLOAT, setPoissonRatio, getPoissonRatio, &model, " label='Poisson ratio XY' min=0.0 step=0.1 precision=4 group='Simulation' ");
@@ -336,3 +340,14 @@ void TW_CALL getSimulationMethod(void *value, void *clientData)
336340
{
337341
*(short *)(value) = (short)((TimeStepTetModel*)clientData)->getSimulationMethod();
338342
}
343+
344+
void TW_CALL setVelocityUpdateMethod(const void *value, void *clientData)
345+
{
346+
const short val = *(const short *)(value);
347+
((TimeStepTetModel*)clientData)->setVelocityUpdateMethod((unsigned int)val);
348+
}
349+
350+
void TW_CALL getVelocityUpdateMethod(void *value, void *clientData)
351+
{
352+
*(short *)(value) = (short)((TimeStepTetModel*)clientData)->getVelocityUpdateMethod();
353+
}

‎Demos/ClothDemo/TimeStepTriangleModel.cpp‎

Lines changed: 21 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
#include "TimeStepTriangleModel.h"
22
#include "Demos/Utils/TimeManager.h"
33
#include "PositionBasedDynamics/PositionBasedDynamics.h"
4+
#include "PositionBasedDynamics/TimeIntegration.h"
45

56
using namespace PBD;
67
using namespace std;
@@ -19,17 +20,25 @@ void TimeStepTriangleModel::step(TriangleModel &model)
1920
{
2021
TimeManager *tm = TimeManager::getCurrent ();
2122
const float h = tm->getTimeStepSize();
23+
ParticleData &pd = model.getParticleMesh().getVertexData();
2224

2325
clearAccelerations(model);
24-
semiImplicitEulerStep(model, h);
26+
for (unsigned int i = 0; i < pd.size(); i++)
27+
{
28+
pd.getLastPosition(i) = pd.getOldPosition(i);
29+
pd.getOldPosition(i) = pd.getPosition(i);
30+
TimeIntegration::semiImplicitEuler(h, pd.getMass(i), pd.getPosition(i), pd.getVelocity(i), pd.getAcceleration(i));
31+
}
2532

2633
constraintProjection(model);
2734

28-
// Update velocities
29-
ParticleData &pd = model.getParticleMesh().getVertexData();
35+
// Update velocities
3036
for (unsigned int i = 0; i < pd.size(); i++)
3137
{
32-
pd.getVelocity(i) = 1.0f / h * (pd.getPosition(i) - pd.getLastPosition(i));
38+
if (m_velocityUpdateMethod == 0)
39+
TimeIntegration::velocityUpdateFirstOrder(h, pd.getMass(i), pd.getPosition(i), pd.getOldPosition(i), pd.getVelocity(i));
40+
else
41+
TimeIntegration::velocityUpdateSecondOrder(h, pd.getMass(i), pd.getPosition(i), pd.getOldPosition(i), pd.getLastPosition(i), pd.getVelocity(i));
3342
}
3443

3544
// compute new time
@@ -54,27 +63,6 @@ void TimeStepTriangleModel::clearAccelerations(TriangleModel &model)
5463
}
5564
}
5665

57-
void TimeStepTriangleModel::semiImplicitEulerStep(TriangleModel &model, const float h)
58-
{
59-
ParticleData &pd = model.getParticleMesh().getVertexData();
60-
61-
// h * x'
62-
for (unsigned int i=0; i < pd.size(); i++)
63-
{
64-
if (pd.getMass(i) != 0.0)
65-
{
66-
Eigen::Vector3f &pos = pd.getPosition(i);
67-
Eigen::Vector3f &vel = pd.getVelocity(i);
68-
const Eigen::Vector3f &accel = pd.getAcceleration(i);
69-
Eigen::Vector3f &lastPos = pd.getLastPosition(i);
70-
lastPos = pos;
71-
72-
vel += accel * h;
73-
pos += vel * h;
74-
}
75-
}
76-
}
77-
7866
void TimeStepTriangleModel::reset(TriangleModel &model)
7967
{
8068

@@ -115,9 +103,9 @@ void TimeStepTriangleModel::constraintProjection(TriangleModel &model)
115103

116104
if (res)
117105
{
118-
if (invMass1 != 0)
106+
if (invMass1 != 0.0f)
119107
x1 += corr1;
120-
if (invMass2 != 0)
108+
if (invMass2 != 0.0f)
121109
x2 += corr2;
122110
}
123111
}
@@ -153,11 +141,11 @@ void TimeStepTriangleModel::constraintProjection(TriangleModel &model)
153141

154142
if (res)
155143
{
156-
if (invMass1 != 0)
144+
if (invMass1 != 0.0f)
157145
x1 += corr1;
158-
if (invMass2 != 0)
146+
if (invMass2 != 0.0f)
159147
x2 += corr2;
160-
if (invMass3 != 0)
148+
if (invMass3 != 0.0f)
161149
x3 += corr3;
162150
}
163151
}
@@ -192,11 +180,11 @@ void TimeStepTriangleModel::constraintProjection(TriangleModel &model)
192180

193181
if (res)
194182
{
195-
if (invMass1 != 0)
183+
if (invMass1 != 0.0f)
196184
x1 += corr1;
197-
if (invMass2 != 0)
185+
if (invMass2 != 0.0f)
198186
x2 += corr2;
199-
if (invMass3 != 0)
187+
if (invMass3 != 0.0f)
200188
x3 += corr3;
201189
}
202190
}

‎Demos/ClothDemo/TimeStepTriangleModel.h‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,11 +9,11 @@ namespace PBD
99
class TimeStepTriangleModel
1010
{
1111
protected:
12+
unsigned int m_velocityUpdateMethod;
1213
unsigned int m_simulationMethod;
1314
unsigned int m_bendingMethod;
1415

1516
void clearAccelerations(TriangleModel &model);
16-
void semiImplicitEulerStep(TriangleModel &model, const float h);
1717
void constraintProjection(TriangleModel &model);
1818

1919
public:
@@ -23,6 +23,8 @@ namespace PBD
2323
void step(TriangleModel &model);
2424
void reset(TriangleModel &model);
2525

26+
unsigned int getVelocityUpdateMethod() const { return m_velocityUpdateMethod; }
27+
void setVelocityUpdateMethod(unsigned int val) { m_velocityUpdateMethod = val; }
2628
unsigned int getSimulationMethod() const { return m_simulationMethod; }
2729
void setSimulationMethod(unsigned int val) { m_simulationMethod = val; }
2830
unsigned int getBendingMethod() const { return m_bendingMethod; }

‎Demos/ClothDemo/TriangleModel.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,6 +72,7 @@ void TriangleModel::reset()
7272
const Eigen::Vector3f& x0 = pd.getPosition0(i);
7373
pd.getPosition(i) = x0;
7474
pd.getLastPosition(i) = pd.getPosition(i);
75+
pd.getOldPosition(i) = pd.getPosition(i);
7576
pd.getVelocity(i).setZero();
7677
pd.getAcceleration(i).setZero();
7778
}

‎Demos/ClothDemo/main.cpp‎

Lines changed: 21 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,9 @@ void TW_CALL setBendingMethod(const void *value, void *clientData);
4848
void TW_CALL getBendingMethod(void *value, void *clientData);
4949
void TW_CALL setSimulationMethod(const void *value, void *clientData);
5050
void TW_CALL getSimulationMethod(void *value, void *clientData);
51+
void TW_CALL setVelocityUpdateMethod(const void *value, void *clientData);
52+
void TW_CALL getVelocityUpdateMethod(void *value, void *clientData);
53+
5154

5255
TriangleModel model;
5356
TimeStepTriangleModel simulation;
@@ -80,8 +83,10 @@ int main( int argc, char **argv )
8083

8184
TwAddVarRW(MiniGL::getTweakBar(), "Pause", TW_TYPE_BOOLCPP, &pause, " label='Pause' group=Simulation key=SPACE ");
8285
TwAddVarCB(MiniGL::getTweakBar(), "TimeStepSize", TW_TYPE_FLOAT, setTimeStep, getTimeStep, &model, " label='Time step size' min=0.0 max = 0.1 step=0.001 precision=4 group=Simulation ");
83-
TwType enumType = TwDefineEnum("SimulationMethodType", NULL, 0);
84-
TwAddVarCB(MiniGL::getTweakBar(), "SimulationMethod", enumType, setSimulationMethod, getSimulationMethod, &simulation, " label='Simulation method' enum='0 {None}, 1 {Distance constraints}, 2 {FEM based PBD}, 3 {Strain based dynamics}' group=Simulation");
86+
TwType enumType = TwDefineEnum("VelocityUpdateMethodType", NULL, 0);
87+
TwAddVarCB(MiniGL::getTweakBar(), "VelocityUpdateMethod", enumType, setVelocityUpdateMethod, getVelocityUpdateMethod, &simulation, " label='Velocity update method' enum='0 {First Order Update}, 1 {Second Order Update}' group=Simulation");
88+
TwType enumType2 = TwDefineEnum("SimulationMethodType", NULL, 0);
89+
TwAddVarCB(MiniGL::getTweakBar(), "SimulationMethod", enumType2, setSimulationMethod, getSimulationMethod, &simulation, " label='Simulation method' enum='0 {None}, 1 {Distance constraints}, 2 {FEM based PBD}, 3 {Strain based dynamics}' group=Simulation");
8590
TwAddVarCB(MiniGL::getTweakBar(), "Stiffness", TW_TYPE_FLOAT, setStiffness, getStiffness, &model, " label='Stiffness' min=0.0 step=0.1 precision=4 group='Distance constraints' ");
8691
TwAddVarCB(MiniGL::getTweakBar(), "XXStiffness", TW_TYPE_FLOAT, setXXStiffness, getXXStiffness, &model, " label='Stiffness XX' min=0.0 step=0.1 precision=4 group='Strain based dynamics' ");
8792
TwAddVarCB(MiniGL::getTweakBar(), "YYStiffness", TW_TYPE_FLOAT, setYYStiffness, getYYStiffness, &model, " label='Stiffness YY' min=0.0 step=0.1 precision=4 group='Strain based dynamics' ");
@@ -93,8 +98,8 @@ int main( int argc, char **argv )
9398
TwAddVarCB(MiniGL::getTweakBar(), "YXPoissonRatioFEM", TW_TYPE_FLOAT, setYXPoissonRatio, getYXPoissonRatio, &model, " label='Poisson ratio YX' min=0.0 step=0.1 precision=4 group='FEM based PBD' ");
9499
TwAddVarCB(MiniGL::getTweakBar(), "NormalizeStretch", TW_TYPE_BOOL32, setNormalizeStretch, getNormalizeStretch, &model, " label='Normalize stretch' group='Strain based dynamics' ");
95100
TwAddVarCB(MiniGL::getTweakBar(), "NormalizeShear", TW_TYPE_BOOL32, setNormalizeShear, getNormalizeShear, &model, " label='Normalize shear' group='Strain based dynamics' ");
96-
TwType enumType2 = TwDefineEnum("BendingMethodType", NULL, 0);
97-
TwAddVarCB(MiniGL::getTweakBar(), "BendingMethod", enumType2, setBendingMethod, getBendingMethod, &simulation, " label='Bending method' enum='0 {None}, 1 {Dihedral angle}, 2 {Isometric bending}' group=Bending");
101+
TwType enumType3 = TwDefineEnum("BendingMethodType", NULL, 0);
102+
TwAddVarCB(MiniGL::getTweakBar(), "BendingMethod", enumType3, setBendingMethod, getBendingMethod, &simulation, " label='Bending method' enum='0 {None}, 1 {Dihedral angle}, 2 {Isometric bending}' group=Bending");
98103
TwAddVarCB(MiniGL::getTweakBar(), "BendingStiffness", TW_TYPE_FLOAT, setBendingStiffness, getBendingStiffness, &model, " label='Bending stiffness' min=0.0 step=0.01 precision=4 group=Bending ");
99104

100105
glutMainLoop ();
@@ -421,3 +426,15 @@ void TW_CALL getSimulationMethod(void *value, void *clientData)
421426
{
422427
*(short *)(value) = (short)((TimeStepTriangleModel*)clientData)->getSimulationMethod();
423428
}
429+
430+
void TW_CALL setVelocityUpdateMethod(const void *value, void *clientData)
431+
{
432+
const short val = *(const short *)(value);
433+
((TimeStepTriangleModel*)clientData)->setVelocityUpdateMethod((unsigned int)val);
434+
}
435+
436+
void TW_CALL getVelocityUpdateMethod(void *value, void *clientData)
437+
{
438+
*(short *)(value) = (short)((TimeStepTriangleModel*)clientData)->getVelocityUpdateMethod();
439+
}
440+

0 commit comments

Comments
 (0)