Class Elasticity_Kee2023
Defined in File Elasticity_Kee2023.h
Nested Relationships
Inheritance Relationships
Base Type
public SPH::NonPressureForceBase(Class NonPressureForceBase)
Class Documentation
-
class Elasticity_Kee2023 : public SPH::NonPressureForceBase
This class implements the elasticity solver by Kee et al. [K23].
References:
[K23] Min Hyung Kee, Kiwon Um, HyunMo Kang, JungHyun Han. An Optimization-based SPH Solver for Simulation of Hyperelastic Solids. Computer Graphics Forum, 42: 225-233, 2023. https://doi.org/10.1111/cgf.14756
Public Functions
-
Elasticity_Kee2023(FluidModel *model)
-
virtual ~Elasticity_Kee2023(void)
-
inline virtual std::string getMethodName()
returns the name of the method
-
virtual void step()
Perform a step of the elasticity solver.
-
virtual void reset()
-
virtual void performNeighborhoodSearchSort()
-
virtual void saveState(BinaryFileWriter &binWriter)
-
virtual void loadState(BinaryFileReader &binReader)
Public Static Functions
-
static inline NonPressureForceBase *creator(FluidModel *model)
Public Static Attributes
-
static std::string METHOD_NAME = "Kee et al. 2023"
-
static int YOUNGS_MODULUS = -1
-
static int POISSON_RATIO = -1
-
static int FIXED_BOX_MIN = -1
-
static int FIXED_BOX_MAX = -1
-
static int FIXED_BOX2_MIN = -1
-
static int FIXED_BOX2_MAX = -1
-
static int ALPHA = -1
-
static int MAX_NEIGHBORS = -1
-
static int SOLVER_TYPE = -1
-
static int LBFGS_WINDOW_SIZE = -1
-
static int MATERIAL_TYPE = -1
-
static int MAX_ITER = -1
-
static int MAX_ERROR = -1
-
static int MAX_ITER_CG = -1
-
static int TOL_CG = -1
-
static int MAX_LS_ITER = -1
-
static int LS_ARMIJO_PARAM = -1
-
static int LS_BETA = -1
-
static int USE_LINE_SEARCH = -1
-
static int ENUM_SOLVER_NEWTON = -1
-
static int ENUM_SOLVER_LBFGS = -1
-
static int ENUM_MATERIAL_STABLE_NEOHOOKEAN = -1
-
static int ENUM_MATERIAL_COROTATED = -1
Protected Types
-
typedef Eigen::SimplicialLLT<Eigen::SparseMatrix<double>, Eigen::Lower, Eigen::AMDOrdering<int>> SolverLLT
Protected Functions
-
void determineFixedParticles()
Mark all particles in the bounding box as fixed.
-
std::string computeMD5(const unsigned int objIndex)
Compute an MD4 check sum using the neighborhood structure in order to recognize known particle models (cache).
-
void initValues()
Initialize the particle neighborhoods in the reference configuration. Fix particles which lie in the user-defined bounding box. Find out if there are multiple separate objects in the phase. Finally, compute kernel gradient correction matrices and factorization.
-
void initSystem()
Initialize the solver for the linear system by either computing a factorization or loading a factorization from the cache.
Compute the factorization of the linear system matrix. This is only done once at the beginning of the simulation.
-
void findObjects()
Find separate objects by object id.
-
void computeMatrixL()
Compute kernel gradient correction matrices (Eq. 8).
-
void precomputeValues()
Precompute V_j * gradW(xi0 - xj0) for all neighbor pairs. Depends only on rest positions and volumes — constant for the simulation. Called once at init (and after any neighbor sort). AVX version packs 8 neighbors per Vector3f8 entry (Kugelstadt pattern).
-
void stepElasticitySolver()
Solve the optimization problem for elastic forces.
min_x (1/2)(x - x_tilde)^T M (x - x_tilde) + dt^2 * Psi(x)
Newton: A * dx = -gradient, where A = M + dt^2 * D^T * K * D (prefactored Cholesky) L-BFGS: quasi-Newton with prefactored A as initial inverse Hessian H_0, secant pairs (s_k, y_k) in a circular queue of size m (window size)
-
void computeXTilde(ElasticObject *obj)
Compute predicted position: xTilde = x + dt * v
-
void updateVelocity(ElasticObject *obj, Real fdt)
Update velocities from solved positions: v = (xk - x) / dt
-
Real computeEnergy(ElasticObject *obj)
Evaluate the objective energy at the current xk (no gradient computation). Used for backtracking line search.
-
Real computePsi(const Matrix3r &F, const Matrix3r &R) const
Compute energy density Psi based on material type.
Stable Neo-Hookean (Smith et al. 2018 Eq. 14): Psi = (μ/2)(I_C - 3) + (λ/2)(J - α)² - (μ/2)log(I_C + 1)
Co-rotated: Psi = μ||F - R||² + (λ/2)(tr(R^T F) - 3)²
-
Real computeEnergyAndGradient(ElasticObject *obj)
Compute the energy and gradient of the objective function.
Energy: E(x) = (1/2)(x - x_tilde)^T M (x - x_tilde) + dt^2 * sum_i V_i * Psi(F_i) + (1/2) x^T H^T K2 H x
Energy density Psi depends on the selected material type:
Stable Neo-Hookean (Smith et al. 2018): Psi = (mu/2)(I_C - 3) + (lambda/2)(J - alpha)^2 - (mu/2)log(I_C + 1)
Co-rotated: Psi = mu * ||F - R||^2_F + (lambda/2) * (J - 1)^2
Gradient (stored in m_gradient) of the objective: g_i = m_i * (sol_i - x_tilde_i)
dt^2 * V_i * sum_j V_j * (P_i^T * L_i + P_j^T * L_j) * gradW(x_i0 - x_j0)
(H^T K2 H xk)_i
Returns the total energy E(x).
-
void computeHessian(ElasticObject *obj)
Compute the Hessian of the elastic energy w.r.t. positions.
For Newton: should be updated every iteration. For L-BFGS: constant proxy K = (2*mu + lambda) * V_i is prefactored in initFactorization() and used as H_0.
-
void computeCorotatedHessian9x9(ElasticObject *obj)
Compute the per-particle 9x9 Hessian d²ψ/d(vecF)² for the co-rotated model.
Uses iARAP approach (Lin et al. 2022) for full decomposition:
Singular values extracted from quartic roots
R computed from Cauchy-Green invariants
V computed analytically from S = R^T F
Co-rotated energy density: ψ = μ||F-R||² + (λ/2)(tr(R^T F) - 3)²
The 9x9 Hessian decomposes into: H = 2μ * H_ARAP + λ * H_vol
ARAP part: H_ARAP = I₉ - Σₖ λₖ (tₖ tₖᵀ) λₖ = 2/(σₐ+σᵦ), clamped to 1 when σₐ+σᵦ < 2
PSD projection: eigendecompose 9x9, clamp negative eigenvalues to 0.
-
void computeStableNeoHookeanHessian9x9(ElasticObject *obj)
Compute the per-particle 9x9 Hessian for Stable Neo-Hookean (Smith et al. 2018).
Uses iARAP decomposition (Lin et al. 2022) for U, V, sigma.
Energy: ψ = μ/2 (I_C - 3) + λ/2 (J - α)² - μ/2 log(I_C + 1) where α = 1 + μ/λ - μ/(4λ)
Hessian decomposes into 4 terms (Eq. 22): A = μ_T I + μ_F f⊗f + λ g⊗g + γ_H H_vol
Direct approach: add each term separately, use volume eigendecomposition for H_vol.
-
void computeNewtonPreconditioner(ElasticObject *obj)
Assemble the 3N×3N Newton system matrix and Cholesky factorize.
A[3j+b, 3k+c] = m_j δ_{jk}δ_{bc}
dt² Σᵢ Vᵢ Σ_{a,a’} D[3i+a,j] H_i[3b+a,3c+a’] D[3i+a’,k]
HTH[j,k] δ_{bc}
where H_i is the per-particle 9×9 Hessian, D is the deformation gradient operator, and HTH is the zero-energy mode control matrix. Compute block-diagonal preconditioner for Newton PCG. Extracts and inverts 3×3 diagonal blocks of A = M + D^T·K·D + H^T·K_ze·H. Note: K already includes dt² * Vi (pre-multiplied in computeCorotatedHessian9x9)
-
void newtonMatvec(ElasticObject *obj)
Compute A*x = M*x + H^T*K_ze*H*x + dt²*D^T*K*D*x for Newton PCG.
-
int matFreePCG(ElasticObject *obj)
Solve A*dx = -gradient using Preconditioned Conjugate Gradient. Uses block-diagonal preconditioner (3×3 blocks). Writes result into obj->m_dx_avx / obj->m_dx.
-
Real newtonSolve(ElasticObject *obj, int &cgIter)
Newton solve: compute energy+gradient, assemble true Hessian, solve 3N×3N system. Returns dx in obj->m_dx, gradient in obj->m_gradient. Returns energy at current xk.
-
void prefactorizedLLTSolve(ElasticObject *obj)
Solve the linear system A * x = b using the prefactored Cholesky decomposition.
Input: dx contains b (right-hand side) Output: dx contains the solution
Forward substitution: L * y = P * b Backward substitution: L^T * z = y Inverse permutation: dx = P^T * z
-
Real lbfgsSolve(ElasticObject *obj)
L-BFGS solve: quasi-Newton with prefactored Cholesky as H_0. Uses two-loop recursion with secant pairs in a circular queue. Returns dx in obj->m_dx, gradient in obj->m_gradient. Returns energy at current xk.
Secant pairs are computed from actual position/gradient differences: s_{k-1} = xk - xk_prev (includes line search alpha) y_{k-1} = gradient_current - gradient_prev (via m_gradient)
-
Real lineSearch(ElasticObject *obj, Real energy, int &lsIter)
Backtracking line search with Armijo condition. Reads search direction dx from obj->m_dx, gradient from obj->m_gradient. Returns the step size alpha (0 if line search failed). Does NOT modify obj->m_xk — caller applies the step.
-
Matrix3r computeP(const Matrix3r &F, const Matrix3r &R) const
Compute first Piola-Kirchhoff stress P based on material type.
Stable Neo-Hookean (Smith et al. 2018 Eq. 18): P = μ(1 - 1/(I_C+1)) F + λ(J - α) cof(F)
Co-rotated: P = 2μ(F - R) + λ(tr(R^T F) - 3) R
-
virtual void initParameters()
-
virtual void deferredInit()
This function is called after the simulation scene is loaded and all parameters are initialized. While reading a scene file several parameters can change. The deferred init function should initialize all values which depend on these parameters.
- inline FORCE_INLINE void symMatTimesVec (const Vector6r &M, const Vector3r &v, Vector3r &res)
-
void rotationMatricesToAVXQuaternions()
Protected Attributes
-
std::vector<unsigned int> m_current_to_initial_index
-
std::vector<unsigned int> m_initial_to_current_index
-
std::vector<std::vector<unsigned int>> m_initialNeighbors
-
std::vector<int> m_fixedGroupId
-
int m_maxNeighbors
-
int m_solverType
-
int m_lbfgsWindowSize
-
int m_materialType
-
int m_maxIter
-
int m_maxIterCG
-
int m_maxLSIter
-
bool m_useLineSearch
-
unsigned int m_totalNeighbors
-
std::vector<ElasticObject*> m_objects
-
std::vector<unsigned int> m_precomputed_indices
-
struct ElasticObject
-
struct Factorization