Skip to content

Commit 5b92480

Browse files
Q-MinhQ-Minh
committed
Implement CPU version of VBD integrator
1 parent 6564535 commit 5b92480

21 files changed

+621
-2
lines changed

‎source/pbat/CMakeLists.txt‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -35,4 +35,5 @@ add_subdirectory(geometry)
3535
add_subdirectory(gpu)
3636
add_subdirectory(math)
3737
add_subdirectory(physics)
38-
add_subdirectory(profiling)
38+
add_subdirectory(profiling)
39+
add_subdirectory(sim)

‎source/pbat/gpu/vbd/VbdImplKernels.cu‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -71,7 +71,7 @@ __global__ void MinimizeBackwardEuler(BackwardEulerMinimization BDF)
7171
if (abs(Determinant(Hi)) <= BDF.detHZero) // Skip nearly rank-deficient hessian
7272
return;
7373

74-
xi = xi - (Inverse(Hi) * gi);
74+
xi -= (Inverse(Hi) * gi);
7575

7676
// 4. Commit vertex descent step
7777
ToBuffers(xi, BDF.x, i);

‎source/pbat/sim/CMakeLists.txt‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
target_sources(PhysicsBasedAnimationToolkit_PhysicsBasedAnimationToolkit
2+
PUBLIC
3+
FILE_SET api
4+
FILES
5+
"Sim.h"
6+
)
7+
8+
add_subdirectory(vbd)
9+
add_subdirectory(xpbd)

‎source/pbat/sim/Sim.h‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
#ifndef PBAT_SIM_SIM_H
2+
#define PBAT_SIM_SIM_H
3+
4+
namespace pbat {
5+
namespace sim {
6+
7+
#include "Vbd.h"
8+
#include "Xpbd.h"
9+
10+
} // namespace sim
11+
} // namespace pbat
12+
13+
#endif // PBAT_SIM_SIM_H

‎source/pbat/sim/vbd/CMakeLists.txt‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
target_sources(PhysicsBasedAnimationToolkit_PhysicsBasedAnimationToolkit
2+
PUBLIC
3+
FILE_SET api
4+
FILES
5+
"Vbd.h"
6+
"Data.h"
7+
"Enums.h"
8+
"Kernels.h"
9+
"Integrator.h"
10+
)
11+
target_sources(PhysicsBasedAnimationToolkit_PhysicsBasedAnimationToolkit
12+
PRIVATE
13+
"Data.cpp"
14+
"Kernels.cpp"
15+
"Integrator.cpp"
16+
)

‎source/pbat/sim/vbd/Data.cpp‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
#include "Data.h"
2+
3+
namespace pbat {
4+
namespace sim {
5+
namespace vbd {
6+
7+
} // namespace vbd
8+
} // namespace sim
9+
} // namespace pbat

‎source/pbat/sim/vbd/Data.h‎

Lines changed: 84 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,84 @@
1+
#ifndef PBAT_SIM_VBD_DATA_H
2+
#define PBAT_SIM_VBD_DATA_H
3+
4+
#include "Enums.h"
5+
#include "PhysicsBasedAnimationToolkitExport.h"
6+
#include "pbat/Aliases.h"
7+
8+
namespace pbat {
9+
namespace sim {
10+
namespace vbd {
11+
12+
PBAT_API struct Data
13+
{
14+
public:
15+
Data& VolumeMesh(Eigen::Ref<MatrixX const> const& V, Eigen::Ref<IndexMatrixX const> const& E);
16+
Data&
17+
SurfaceMesh(Eigen::Ref<IndexVectorX const> const& V, Eigen::Ref<IndexMatrixX const> const& F);
18+
Data& Velocity(Eigen::Ref<MatrixX const> const& v);
19+
Data& Acceleration(Eigen::Ref<MatrixX const> const& aext);
20+
Data& Mass(Eigen::Ref<VectorX const> const& m);
21+
Data& Quadrature(
22+
Eigen::Ref<VectorX const> const& wg,
23+
Eigen::Ref<MatrixX const> const& GP,
24+
Eigen::Ref<MatrixX const> const& lame);
25+
Data& VertexAdjacency(
26+
Eigen::Ref<IndexVectorX const> const& GVGp,
27+
Eigen::Ref<IndexVectorX const> const& GVGn,
28+
Eigen::Ref<IndexVectorX const> const& GVGe,
29+
Eigen::Ref<IndexVectorX const> const& GVGilocal);
30+
Data& Partitions(std::vector<std::vector<Index>> const& partitions);
31+
Data& DirichletConstrainedVertices(IndexVectorX const& dbc);
32+
Data& InitializationStrategy(EInitializationStrategy strategy);
33+
Data& RayleighDamping(Scalar kD);
34+
Data& CollisionPenalty(Scalar kC);
35+
Data& NumericalZero(Scalar zero);
36+
37+
public:
38+
IndexVectorX V; ///< Collision vertices
39+
IndexMatrixX F; ///< Collision triangles (on the boundary of T)
40+
IndexMatrixX T; ///< Tetrahedra
41+
42+
MatrixX x; ///< Vertex positions
43+
MatrixX v; ///< Vertex velocities
44+
MatrixX aext; ///< Vertex external accelerations
45+
VectorX m; ///< Vertex masses
46+
47+
MatrixX xt; ///< Previous vertex positions
48+
MatrixX xtilde; ///< Inertial target for vertex positions
49+
MatrixX xchebm2; ///< x^{k-2} used in Chebyshev semi-iterative method
50+
MatrixX xchebm1; ///< x^{k-1} used in Chebyshev semi-iterative method
51+
MatrixX vt; ///< Previous vertex velocities
52+
53+
VectorX wg; ///< |#quad.pts.| quadrature weights
54+
MatrixX GP; ///< |#elem.nodes|x|#dims*#quad.pts.| shape function gradients at quad. pts.
55+
MatrixX lame; ///< 2x|#quad.pts.| Lame coefficients
56+
57+
IndexVectorX GVGp; ///< |#verts+1| prefixes into GVGn
58+
IndexVectorX
59+
GVGn; ///< |# of vertex-quad.pt. edges| neighbours s.t.
60+
///< GVGn[k] for GVGp[i] <= k < GVGp[i+1] gives the quad.pts. involving vertex i
61+
IndexVectorX GVGe; ///< |# of vertex-quad.pt. edges| element indices s.t.
62+
///< GVGe[k] for GVGp[i] <= k < GVGp[i+1] gives the element index of
63+
///< adjacent to vertex i for the neighbouring quad.pt.
64+
IndexVectorX GVGilocal; ///< |# of vertex-quad.pt. edges| local vertex indices s.t.
65+
///< GVGilocal[k] for GVGp[i] <= k < GVGp[i+1] gives the local index of
66+
///< vertex i for the neighbouring quad.pts.
67+
68+
IndexVectorX dbc; ///< Dirichlet constrained vertices
69+
70+
std::vector<std::vector<Index>>
71+
partitions; ///< partitions[c] gives the c^{th} group of vertices which can all be
72+
///< integrated independently in parallel
73+
74+
EInitializationStrategy initializationStrategy; ///< BCD optimization initialization strategy
75+
Scalar kD; ///< Uniform damping coefficient
76+
Scalar kC; ///< Uniform collision penalty
77+
Scalar detHZero; ///< Numerical zero for hessian pseudo-singularity check
78+
};
79+
80+
} // namespace vbd
81+
} // namespace sim
82+
} // namespace pbat
83+
84+
#endif // PBAT_SIM_XPBD_DATA_H

‎source/pbat/sim/vbd/Enums.h‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
#ifndef PBAT_SIM_VBD_ENUMS_H
2+
#define PBAT_SIM_VBD_ENUMS_H
3+
4+
namespace pbat {
5+
namespace sim {
6+
namespace vbd {
7+
8+
enum class EInitializationStrategy {
9+
Position,
10+
Inertia,
11+
KineticEnergyMinimum,
12+
AdaptiveVbd,
13+
AdaptivePbat
14+
};
15+
16+
} // namespace vbd
17+
} // namespace sim
18+
} // namespace pbat
19+
20+
#endif // PBAT_SIM_VBD_ENUMS_H

‎source/pbat/sim/vbd/Integrator.cpp‎

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
#include "Integrator.h"
2+
3+
#include "Kernels.h"
4+
#include "pbat/math/linalg/mini/Mini.h"
5+
#include "pbat/physics/StableNeoHookeanEnergy.h"
6+
7+
#include <tbb/parallel_for.h>
8+
#include <type_traits>
9+
10+
namespace pbat {
11+
namespace sim {
12+
namespace vbd {
13+
14+
Integrator::Integrator(Data&& dataIn) : data(dataIn) {}
15+
16+
void Integrator::Step(Scalar dt, Index iterations, Index substeps, Scalar rho)
17+
{
18+
Scalar sdt = dt / (static_cast<Scalar>(substeps));
19+
Scalar sdt2 = sdt * sdt;
20+
auto const nVertices = data.x.cols();
21+
using IndexType = std::remove_const_t<decltype(nVertices)>;
22+
bool const bUseChebyshevAcceleration = rho > Scalar(0) and rho < Scalar(1);
23+
using namespace math::linalg;
24+
using mini::FromEigen;
25+
using mini::ToEigen;
26+
for (auto s = 0; s < substeps; ++s)
27+
{
28+
// Store previous positions
29+
data.xt = data.x;
30+
// Compute inertial target positions
31+
tbb::parallel_for(IndexType(0), nVertices, [&](auto i) {
32+
auto xtilde = kernels::InertialTarget(
33+
FromEigen(data.xt.col(i).head<3>()),
34+
FromEigen(data.vt.col(i).head<3>()),
35+
FromEigen(data.aext.col(i).head<3>()),
36+
sdt,
37+
sdt2);
38+
data.xtilde.col(i) = ToEigen(xtilde);
39+
});
40+
// Initialize block coordinate descent's, i.e. BCD's, solution
41+
tbb::parallel_for(IndexType(0), nVertices, [&](auto i) {
42+
auto x = kernels::InitialPositionsForSolve(
43+
FromEigen(data.xt.col(i).head<3>()),
44+
FromEigen(data.vt.col(i).head<3>()),
45+
FromEigen(data.v.col(i).head<3>()),
46+
FromEigen(data.aext.col(i).head<3>()),
47+
sdt,
48+
sdt2,
49+
data.initializationStrategy);
50+
data.x.col(i) = ToEigen(x);
51+
});
52+
// Initialize Chebyshev semi-iterative method
53+
Scalar omega{};
54+
Scalar rho2 = rho * rho;
55+
// Minimize Backward Euler, i.e. BDF1, objective
56+
for (auto k = 0; k < iterations; ++k)
57+
{
58+
if (bUseChebyshevAcceleration)
59+
omega = kernels::ChebyshevOmega(k, rho2, omega);
60+
61+
for (auto const& partition : data.partitions)
62+
{
63+
auto const nVerticesInPartition = static_cast<std::size_t>(partition.size());
64+
tbb::parallel_for(std::size_t(0), nVerticesInPartition, [&](auto v) {
65+
auto i = partition[v];
66+
auto begin = data.GVGp[i];
67+
auto end = data.GVGp[i + 1];
68+
// Compute vertex elastic hessian
69+
mini::SMatrix<Scalar, 3, 3> Hi = mini::Zeros<Scalar, 3, 3>();
70+
mini::SVector<Scalar, 3> gi = mini::Zeros<Scalar, 3, 1>();
71+
for (auto n = begin; n < end; ++n)
72+
{
73+
auto ilocal = data.GVGilocal[n];
74+
auto e = data.GVGe[n];
75+
auto g = data.GVGn[n];
76+
auto lamee = data.lame.col(g);
77+
auto wg = data.wg[g];
78+
auto Te = data.T.col(e);
79+
mini::SMatrix<Scalar, 4, 3> GPe = FromEigen(data.GP.block<4, 3>(0, e * 3));
80+
mini::SMatrix<Scalar, 3, 4> xe =
81+
FromEigen(data.x(Eigen::all, Te).block<3, 4>(0, 0));
82+
mini::SMatrix<Scalar, 3, 3> Fe = xe * GPe;
83+
physics::StableNeoHookeanEnergy<3> Psi{};
84+
mini::SVector<Scalar, 9> gF;
85+
mini::SMatrix<Scalar, 9, 9> HF;
86+
Psi.gradAndHessian(Fe, lamee(0), lamee(1), gF, HF);
87+
kernels::AccumulateElasticHessian(ilocal, wg, GPe, HF, Hi);
88+
kernels::AccumulateElasticGradient(ilocal, wg, GPe, gF, gi);
89+
}
90+
// Update vertex position
91+
Scalar m = data.m[i];
92+
mini::SVector<Scalar, 3> xti = FromEigen(data.xt.col(i).head<3>());
93+
mini::SVector<Scalar, 3> xtildei = FromEigen(data.xtilde.col(i).head<3>());
94+
mini::SVector<Scalar, 3> xi = FromEigen(data.x.col(i).head<3>());
95+
kernels::AddDamping(sdt, xti, xi, data.kD, gi, Hi);
96+
kernels::AddInertiaDerivatives(sdt2, m, xtildei, xi, gi, Hi);
97+
kernels::IntegratePositions(gi, Hi, xi, data.detHZero);
98+
data.x.col(i) = ToEigen(xi);
99+
});
100+
}
101+
102+
if (bUseChebyshevAcceleration)
103+
{
104+
tbb::parallel_for(IndexType(0), nVertices, [&](auto i) {
105+
auto xkm2eig = data.xchebm2.col(i).head<3>();
106+
auto xkm1eig = data.xchebm1.col(i).head<3>();
107+
auto xkeig = data.x.col(i).head<3>();
108+
auto xkm2 = FromEigen(xkm2eig);
109+
auto xkm1 = FromEigen(xkm1eig);
110+
auto xk = FromEigen(xkeig);
111+
kernels::ChebyshevUpdate(k, omega, xkm2, xkm1, xk);
112+
});
113+
}
114+
}
115+
// Update velocity
116+
tbb::parallel_for(IndexType(0), nVertices, [&](auto i) {
117+
kernels::IntegrateVelocity(
118+
FromEigen(data.xt.col(i).head<3>()),
119+
FromEigen(data.x.col(i).head<3>()),
120+
sdt);
121+
});
122+
}
123+
}
124+
125+
} // namespace vbd
126+
} // namespace sim
127+
} // namespace pbat

‎source/pbat/sim/vbd/Integrator.h‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
#ifndef PBAT_SIM_VBD_INTEGRATOR_H
2+
#define PBAT_SIM_VBD_INTEGRATOR_H
3+
4+
#include "Data.h"
5+
#include "PhysicsBasedAnimationToolkitExport.h"
6+
#include "pbat/Aliases.h"
7+
8+
namespace pbat {
9+
namespace sim {
10+
namespace vbd {
11+
12+
class Integrator
13+
{
14+
public:
15+
PBAT_API Integrator(Data&& data);
16+
17+
PBAT_API void
18+
Step(Scalar dt, Index iterations, Index substeps = Index{1}, Scalar rho = Scalar{1});
19+
20+
PBAT_API Data data;
21+
};
22+
23+
} // namespace vbd
24+
} // namespace sim
25+
} // namespace pbat
26+
27+
#endif // PBAT_SIM_VBD_INTEGRATOR_H

0 commit comments

Comments
 (0)