Skip to content

Commit 58e600d

Browse files
minh702minh702
authored
Integration of [Kee et al. 2023]: "An Optimization-Based SPH Solver for Simulation of Hyperelastic Solids" (#372)
* WIP: Kee et al. 2023. SNH, CoRot anaytic eigen analysis for PSD fix; L-BFGS; validating line search * Kee2023: L-BFGS buffer renewal, line search with flat energy termination, SVD polar decomposition fix - L-BFGS: runtime window size change with buffer renewal; window=0 for constant Hessian - Line search: backtracking Armijo + flat energy termination (|dE| < 1e-4); convergence check before LS - SVD polar decomposition: flip U.col(2) instead of R.col(2) for correct R under inversion - Non-convergence warning at max iterations - Kugelstadt2021: added neighbor list asymmetry check - Beam.json: updated scene config for validation * [del] remove utilities/version.h * Kee2023: enable USE_AVX build with AVX Cholesky LLT solve (Phase A). * Kee2023: AVX sparse matvec for Dx and Hx in computeEnergyAndGradient (Phase B). * Kee2023: full AVX L-BFGS pipeline with Scalarf8 state, ~3x speedup (Phase C). * Kee2023: rename AVX members with _avx suffix, matching Kugelstadt convention. * Kee2023: precomputeValues for V_j*gradW, AVX Section 3 with precomputed kernel. * Kee2023: fix AVX Section 3 force loop with precomputed V_gradW8, re-guard Newton for AVX build. * Kee2023: AVX sparse matvec for Dx and Hx in computeEnergyAndGradient (Phase B). * Kee2023: full AVX Newton matrix-free PCG (matvec + CG iteration in Scalarf8). * Kee2023: fix MSVC build errors - replace M_PI with literals, use omp critical instead of reduction(max). * Kee2023: remove hardcoded twist, switch to AnimationField (instability issues, see doc). * Kee2023: update Beam validation scene (SNH, narrower fixed range, poisson 0.45). * Kee2023: complete AnimationField-driven rotation for Beam validation (one-end twist).
1 parent 6bff55a commit 58e600d

11 files changed

Lines changed: 4319 additions & 12 deletions

‎SPlisHSPlasH/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -128,12 +128,14 @@ set(ELASTICITY_HEADER_FILES
128128
Elasticity/Elasticity_Becker2009.h
129129
Elasticity/Elasticity_Peer2018.h
130130
Elasticity/Elasticity_Kugelstadt2021.h
131+
Elasticity/Elasticity_Kee2023.h
131132
)
132133

133134
set(ELASTICITY_SOURCE_FILES
134135
Elasticity/Elasticity_Becker2009.cpp
135136
Elasticity/Elasticity_Peer2018.cpp
136137
Elasticity/Elasticity_Kugelstadt2021.cpp
138+
Elasticity/Elasticity_Kee2023.cpp
137139
)
138140

139141
set(UTILS_HEADER_FILES

‎SPlisHSPlasH/Elasticity/Elasticity_Kee2023.cpp‎

Lines changed: 3429 additions & 0 deletions
Large diffs are not rendered by default.

‎SPlisHSPlasH/Elasticity/Elasticity_Kee2023.h‎

Lines changed: 260 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,260 @@
1+
#ifndef __Elasticity_Kee2023_h__
2+
#define __Elasticity_Kee2023_h__
3+
4+
#include "SPlisHSPlasH/Common.h"
5+
#include "SPlisHSPlasH/FluidModel.h"
6+
#include "SPlisHSPlasH/NonPressureForceBase.h"
7+
#if USE_AVX
8+
#include "SPlisHSPlasH/Utilities/AVX_math.h"
9+
#include "SPlisHSPlasH/Utilities/CholeskyAVXSolver.h"
10+
#endif
11+
12+
13+
namespace SPH
14+
{
15+
/** \brief This class implements the elasticity solver
16+
* by Kee et al. [K23].
17+
*
18+
* References:
19+
* - [K23] Kee et al..
20+
* Kee 2023.
21+
*
22+
*/
23+
class Elasticity_Kee2023 : public NonPressureForceBase
24+
{
25+
protected:
26+
27+
struct Factorization
28+
{
29+
Real m_dt;
30+
Real m_mu;
31+
Eigen::SparseMatrix<Real, Eigen::RowMajor> m_DT_K;
32+
Eigen::SparseMatrix<Real, Eigen::RowMajor> m_D;
33+
Eigen::SparseMatrix<Real, Eigen::ColMajor> m_matHTH;
34+
35+
#ifdef USE_AVX
36+
CholeskyAVXSolver *m_cholesky;
37+
Factorization() { m_cholesky = nullptr; }
38+
~Factorization() { delete m_cholesky; }
39+
#else
40+
Factorization() {}
41+
~Factorization() {}
42+
// L-BFGS prefactored Cholesky (N×N, constant proxy)
43+
Eigen::SparseMatrix<Real, Eigen::ColMajor> m_matL;
44+
Eigen::SparseMatrix<Real, Eigen::ColMajor> m_matLT;
45+
Eigen::VectorXi m_permInd;
46+
Eigen::VectorXi m_permInvInd;
47+
#endif
48+
};
49+
50+
struct ElasticObject
51+
{
52+
std::string m_md5;
53+
std::vector<unsigned int> m_particleIndices;
54+
unsigned int m_nFixed;
55+
std::shared_ptr<Factorization> m_factorization;
56+
57+
// Newton: per-particle 9×9 Hessian and block-diagonal preconditioner.
58+
// Same type in both builds (SVD is always scalar).
59+
std::vector<Eigen::Matrix<Real, 9, 9>> m_hessian9x9;
60+
std::vector<Matrix3r, Eigen::aligned_allocator<Matrix3r>> m_pcg_precond;
61+
62+
#ifdef USE_AVX
63+
// AVX state (Scalarf8, 3 active lanes x/y/z).
64+
std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_f_avx; // F = D·xk workspace (3*numParticles)
65+
std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_xk_avx; // current iterate
66+
std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_xTilde_avx; // inertial target
67+
std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_dx_avx; // solver step (also LLT solve RHS/result)
68+
std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_gradient_avx; // ∇E at xk
69+
70+
// Newton CG workspace (Scalarf8, coord-packed x/y/z in lanes 0/1/2)
71+
std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_pcg_r_avx;
72+
std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_pcg_p_avx;
73+
std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_pcg_Ap_avx;
74+
std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_pcg_z_avx;
75+
76+
// L-BFGS secant history
77+
std::vector<std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>>> m_lbfgs_s_avx;
78+
std::vector<std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>>> m_lbfgs_y_avx;
79+
std::vector<Real> m_lbfgs_rho;
80+
std::vector<Real> m_lbfgs_alpha;
81+
std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_lbfgs_last_sol_avx;
82+
std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_lbfgs_q_avx;
83+
int m_lbfgs_count = 0;
84+
#else
85+
std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_f; // F = D·xk workspace
86+
std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_xk; // current iterate
87+
std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_xTilde; // inertial target
88+
std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_dx; // Newton/LBFGS step (also LLT solve RHS/result)
89+
std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_gradient; // ∇E at xk
90+
91+
// Newton PCG workspace
92+
std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_pcg_r;
93+
std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_pcg_p;
94+
std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_pcg_Ap;
95+
std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_pcg_z;
96+
97+
// L-BFGS secant history
98+
std::vector<std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>>> m_lbfgs_s;
99+
std::vector<std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>>> m_lbfgs_y;
100+
std::vector<Real> m_lbfgs_rho;
101+
std::vector<Real> m_lbfgs_alpha;
102+
std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_lbfgs_last_sol;
103+
std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_lbfgs_q;
104+
int m_lbfgs_count = 0;
105+
106+
// Permutation workspace for scalar LLT (manual forward/backward sub)
107+
std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_dx_perm;
108+
#endif
109+
110+
ElasticObject() { m_factorization = nullptr; }
111+
~ElasticObject() { m_factorization = nullptr; }
112+
};
113+
114+
Real m_youngsModulus;
115+
Real m_poissonRatio;
116+
Vector3r m_fixedBoxMin;
117+
Vector3r m_fixedBoxMax;
118+
Vector3r m_fixedBox2Min;
119+
Vector3r m_fixedBox2Max;
120+
121+
// initial particle indices, used to access their original positions
122+
std::vector<unsigned int> m_current_to_initial_index;
123+
std::vector<unsigned int> m_initial_to_current_index;
124+
// initial particle neighborhood
125+
std::vector<std::vector<unsigned int>> m_initialNeighbors;
126+
// volumes in rest configuration
127+
std::vector<Real> m_restVolumes;
128+
std::vector<Matrix3r> m_rotations;
129+
std::vector<Real> m_stress;
130+
std::vector<int> m_fixedGroupId; // 0: free, 1: box1, 2: box2
131+
std::vector<Matrix3r> m_L;
132+
std::vector<Matrix3r> m_F;
133+
std::vector<Matrix3r> m_PL;
134+
Real m_alpha;
135+
int m_maxNeighbors;
136+
int m_solverType; // 0: Newton, 1: LBFGS
137+
int m_lbfgsWindowSize;
138+
int m_materialType; // 0: Stable Neo-Hookean, 1: Co-rotated
139+
int m_maxIter;
140+
Real m_maxError;
141+
int m_maxIterCG; // max CG iterations for Newton linear solve
142+
Real m_tolCG; // CG convergence tolerance
143+
int m_maxLSIter;
144+
Real m_lsArmijoParam;
145+
Real m_lsBeta;
146+
bool m_useLineSearch;
147+
unsigned int m_totalNeighbors;
148+
std::vector<ElasticObject*> m_objects;
149+
Real m_lambda;
150+
Real m_mu;
151+
152+
// Precomputed V_j * gradW(xi0 - xj0) per neighbor pair.
153+
// Scalar format (always available, used by Newton CG).
154+
std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_precomp_V_gradW;
155+
std::vector<unsigned int> m_precomputed_indices;
156+
#ifdef USE_AVX
157+
// AVX-packed format (8 neighbors per entry, used by L-BFGS force loops).
158+
std::vector<Vector3f8, Eigen::aligned_allocator<Vector3f8>> m_precomp_V_gradW8;
159+
std::vector<unsigned int> m_precomputed_indices8;
160+
#endif
161+
162+
#ifdef USE_AVX
163+
typedef Eigen::SimplicialLLT<Eigen::SparseMatrix<double>, Eigen::Lower, Eigen::AMDOrdering<int>> SolverLLT;
164+
#else
165+
typedef Eigen::SimplicialLLT<Eigen::SparseMatrix<double>, Eigen::Lower, Eigen::AMDOrdering<int>> SolverLLT;
166+
#endif
167+
168+
void determineFixedParticles();
169+
std::string computeMD5(const unsigned int objIndex);
170+
void initValues();
171+
void initSystem();
172+
void initFactorization(std::shared_ptr<Factorization> factorization, std::vector<unsigned int> &particleIndices, const unsigned int nFixed, const Real dt, const Real mu);
173+
void findObjects();
174+
void computeMatrixL();
175+
void precomputeValues();
176+
177+
void stepElasticitySolver();
178+
179+
void computeXTilde(ElasticObject* obj);
180+
void updateVelocity(ElasticObject* obj, Real fdt);
181+
Real computeEnergy(ElasticObject* obj);
182+
Real computePsi(const Matrix3r& F, const Matrix3r& R) const;
183+
Real computeEnergyAndGradient(ElasticObject* obj);
184+
void computeHessian(ElasticObject* obj);
185+
void computeCorotatedHessian9x9(ElasticObject* obj);
186+
void computeStableNeoHookeanHessian9x9(ElasticObject* obj);
187+
void computeNewtonPreconditioner(ElasticObject* obj);
188+
void newtonMatvec(ElasticObject* obj);
189+
int matFreePCG(ElasticObject* obj);
190+
Real newtonSolve(ElasticObject* obj, int& cgIter);
191+
void prefactorizedLLTSolve(ElasticObject* obj);
192+
Real lbfgsSolve(ElasticObject* obj);
193+
Real lineSearch(ElasticObject* obj, Real energy, int& lsIter);
194+
195+
Matrix3r computeP(const Matrix3r& F, const Matrix3r& R) const;
196+
197+
virtual void initParameters();
198+
/** This function is called after the simulation scene is loaded and all
199+
* parameters are initialized. While reading a scene file several parameters
200+
* can change. The deferred init function should initialize all values which
201+
* depend on these parameters.
202+
*/
203+
virtual void deferredInit();
204+
205+
//////////////////////////////////////////////////////////////////////////
206+
// multiplication of symmetric matrix, represented by a 6D vector, and a
207+
// 3D vector
208+
//////////////////////////////////////////////////////////////////////////
209+
FORCE_INLINE void symMatTimesVec(const Vector6r & M, const Vector3r & v, Vector3r &res)
210+
{
211+
res[0] = M[0] * v[0] + M[3] * v[1] + M[4] * v[2];
212+
res[1] = M[3] * v[0] + M[1] * v[1] + M[5] * v[2];
213+
res[2] = M[4] * v[0] + M[5] * v[1] + M[2] * v[2];
214+
}
215+
216+
void rotationMatricesToAVXQuaternions();
217+
218+
public:
219+
static std::string METHOD_NAME;
220+
static int YOUNGS_MODULUS;
221+
static int POISSON_RATIO;
222+
static int FIXED_BOX_MIN;
223+
static int FIXED_BOX_MAX;
224+
static int FIXED_BOX2_MIN;
225+
static int FIXED_BOX2_MAX;
226+
static int ALPHA;
227+
static int MAX_NEIGHBORS;
228+
static int SOLVER_TYPE;
229+
static int LBFGS_WINDOW_SIZE;
230+
static int MATERIAL_TYPE;
231+
static int MAX_ITER;
232+
static int MAX_ERROR;
233+
static int MAX_ITER_CG;
234+
static int TOL_CG;
235+
static int MAX_LS_ITER;
236+
static int LS_ARMIJO_PARAM;
237+
static int LS_BETA;
238+
static int USE_LINE_SEARCH;
239+
240+
static int ENUM_SOLVER_NEWTON;
241+
static int ENUM_SOLVER_LBFGS;
242+
static int ENUM_MATERIAL_STABLE_NEOHOOKEAN;
243+
static int ENUM_MATERIAL_COROTATED;
244+
245+
Elasticity_Kee2023(FluidModel *model);
246+
virtual ~Elasticity_Kee2023(void);
247+
248+
static NonPressureForceBase* creator(FluidModel* model) { return new Elasticity_Kee2023(model); }
249+
virtual std::string getMethodName() { return METHOD_NAME; }
250+
251+
virtual void step();
252+
virtual void reset();
253+
virtual void performNeighborhoodSearchSort();
254+
255+
virtual void saveState(BinaryFileWriter &binWriter);
256+
virtual void loadState(BinaryFileReader &binReader);
257+
};
258+
}
259+
260+
#endif

‎SPlisHSPlasH/Elasticity/Elasticity_Kugelstadt2021.cpp‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -308,6 +308,29 @@ void Elasticity_Kugelstadt2021::initValues()
308308
}
309309
}
310310

311+
// Check neighbor list symmetry
312+
{
313+
unsigned int maxActualNeighbors = 0;
314+
int asymmetricPairs = 0;
315+
for (unsigned int i = 0; i < numParticles; i++)
316+
{
317+
if (m_initialNeighbors[i].size() > maxActualNeighbors)
318+
maxActualNeighbors = (unsigned int)m_initialNeighbors[i].size();
319+
for (size_t jn = 0; jn < m_initialNeighbors[i].size(); jn++)
320+
{
321+
unsigned int j = m_initialNeighbors[i][jn];
322+
bool found = false;
323+
for (size_t kn = 0; kn < m_initialNeighbors[j].size(); kn++)
324+
{
325+
if (m_initialNeighbors[j][kn] == i) { found = true; break; }
326+
}
327+
if (!found) asymmetricPairs++;
328+
}
329+
}
330+
LOG_INFO << "Neighbor stats: maxNeighbors=" << maxActualNeighbors
331+
<< ", asymmetricPairs=" << asymmetricPairs;
332+
}
333+
311334
// mark all particles in the bounding box as fixed
312335
determineFixedParticles();
313336

‎SPlisHSPlasH/NonPressureForceRegistration.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
#include "Elasticity/Elasticity_Becker2009.h"
1818
#include "Elasticity/Elasticity_Peer2018.h"
1919
#include "Elasticity/Elasticity_Kugelstadt2021.h"
20+
#include "Elasticity/Elasticity_Kee2023.h"
2021

2122
#include "SurfaceTension/SurfaceTension_Becker2007.h"
2223
#include "SurfaceTension/SurfaceTension_Akinci2013.h"
@@ -39,6 +40,7 @@ void Simulation::registerNonpressureForces()
3940
addElasticityMethod(Elasticity_Becker2009::METHOD_NAME, Elasticity_Becker2009::creator);
4041
addElasticityMethod(Elasticity_Peer2018::METHOD_NAME, Elasticity_Peer2018::creator);
4142
addElasticityMethod(Elasticity_Kugelstadt2021::METHOD_NAME, Elasticity_Kugelstadt2021::creator);
43+
addElasticityMethod(Elasticity_Kee2023::METHOD_NAME, Elasticity_Kee2023::creator);
4244

4345
addSurfaceTensionMethod("None", [](FluidModel*) -> NonPressureForceBase* { return nullptr; });
4446
addSurfaceTensionMethod(SurfaceTension_Becker2007::METHOD_NAME, SurfaceTension_Becker2007::creator);

0 commit comments

Comments
 (0)