Integration of [Kee et al. 2023]: "An Optimization-Based SPH Solver for Simulation of Hyperelastic Solids"

SPlisHSPlasH/CMakeLists.txt

@@ -128,12 +128,14 @@ set(ELASTICITY_HEADER_FILES
 	Elasticity/Elasticity_Becker2009.h	
 	Elasticity/Elasticity_Peer2018.h
 	Elasticity/Elasticity_Kugelstadt2021.h
+	Elasticity/Elasticity_Kee2023.h
 	)
 
 set(ELASTICITY_SOURCE_FILES
 	Elasticity/Elasticity_Becker2009.cpp
 	Elasticity/Elasticity_Peer2018.cpp
 	Elasticity/Elasticity_Kugelstadt2021.cpp
+	Elasticity/Elasticity_Kee2023.cpp
 	)
 
 set(UTILS_HEADER_FILES

SPlisHSPlasH/Elasticity/Elasticity_Kee2023.h

@@ -0,0 +1,235 @@
+#ifndef __Elasticity_Kee2023_h__
+#define __Elasticity_Kee2023_h__
+
+#include "SPlisHSPlasH/Common.h"
+#include "SPlisHSPlasH/FluidModel.h"
+#include "SPlisHSPlasH/NonPressureForceBase.h"
+#if USE_AVX
+#include "SPlisHSPlasH/Utilities/AVX_math.h"
+#include "SPlisHSPlasH/Utilities/CholeskyAVXSolver.h"
+#endif
+
+
+namespace SPH
+{
+	/** \brief This class implements the elasticity solver
+	* by Kee et al. [K23].
+	*
+	* References:
+	* - [K23] Kee et al..
+	* Kee 2023.
+	*
+	*/
+	class Elasticity_Kee2023 : public NonPressureForceBase
+	{
+	protected:
+
+		struct Factorization
+		{
+			Real m_dt;
+			Real m_mu;
+			Eigen::SparseMatrix<Real, Eigen::RowMajor> m_DT_K;
+			Eigen::SparseMatrix<Real, Eigen::RowMajor> m_D;
+			Eigen::SparseMatrix<Real, Eigen::ColMajor> m_matHTH;
+
+#ifdef USE_AVX
+			CholeskyAVXSolver *m_cholesky;
+			Factorization() { m_cholesky = nullptr; }
+			~Factorization() { delete m_cholesky; }
+#else
+			Factorization() {}
+			~Factorization() {}
+			// L-BFGS prefactored Cholesky (N×N, constant proxy)
+			Eigen::SparseMatrix<Real, Eigen::ColMajor> m_matL;
+			Eigen::SparseMatrix<Real, Eigen::ColMajor> m_matLT;
+			Eigen::VectorXi m_permInd;
+			Eigen::VectorXi m_permInvInd;
+#endif
+		};
+
+		struct ElasticObject
+		{
+			std::string m_md5;
+			std::vector<unsigned int> m_particleIndices;
+			unsigned int m_nFixed;
+
+			std::shared_ptr<Factorization> m_factorization;
+#ifdef USE_AVX
+			VectorXr m_rhs;
+			VectorXr m_sol;
+			std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_dx;
+			std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_f_avx;
+			std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_sol_avx;
+			std::vector<Quaternion8f, AlignmentAllocator<Quaternion8f, 32>> m_quats_avx;
+#else
+			std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_f;
+			std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_xk;
+			std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_xTilde;
+			std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_dx;
+			std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_dx_perm;
+			std::vector<Quaternionr, Eigen::aligned_allocator<Quaternionr>> m_quats;
+			std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_gradient;  // gradient at current xk (for line search)
+
+			// Newton: per-particle 9×9 Hessian (K_i = d²ψ/dvec(F)²)
+			std::vector<Eigen::Matrix<Real, 9, 9>> m_hessian9x9;
+
+			// Newton PCG workspace
+			std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_pcg_r;   // residual
+			std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_pcg_p;   // search direction
+			std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_pcg_Ap;  // A * p
+			std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_pcg_z;   // preconditioned residual
+			std::vector<Matrix3r, Eigen::aligned_allocator<Matrix3r>> m_pcg_precond;  // block-diagonal preconditioner (inverted 3x3 blocks)
+
+			// L-BFGS secant history (circular queue)
+			std::vector<std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>>> m_lbfgs_s;  // position differences s_k = x_{k+1} - x_k
+			std::vector<std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>>> m_lbfgs_y;  // gradient differences y_k = g_{k+1} - g_k
+			std::vector<Real> m_lbfgs_rho;     // 1 / (y_k^T s_k)
+			std::vector<Real> m_lbfgs_alpha;   // temporary for two-loop recursion
+			std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_lbfgs_last_sol;   // previous sol for s_k computation
+			std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_lbfgs_q;          // temporary for two-loop recursion
+			int m_lbfgs_count = 0;             // number of stored secant pairs
+#endif
+
+			ElasticObject() { m_factorization = nullptr; }
+			~ElasticObject() { m_factorization = nullptr; }
+		};
+
+		Real m_youngsModulus;
+		Real m_poissonRatio;
+		Vector3r m_fixedBoxMin;
+		Vector3r m_fixedBoxMax;
+		Vector3r m_fixedBox2Min;
+		Vector3r m_fixedBox2Max;
+
+		// initial particle indices, used to access their original positions
+		std::vector<unsigned int> m_current_to_initial_index;
+		std::vector<unsigned int> m_initial_to_current_index;
+		// initial particle neighborhood
+		std::vector<std::vector<unsigned int>> m_initialNeighbors;
+		// volumes in rest configuration
+		std::vector<Real> m_restVolumes;
+		std::vector<Matrix3r> m_rotations;
+		std::vector<Real> m_stress;
+		std::vector<int> m_fixedGroupId;		// 0: free, 1: box1, 2: box2
+		std::vector<Matrix3r> m_L;
+		std::vector<Matrix3r> m_F;
+		std::vector<Matrix3r> m_PL;
+		Real m_alpha;
+		int m_maxNeighbors;
+		int m_solverType;			// 0: Newton, 1: LBFGS
+		int m_lbfgsWindowSize;
+		int m_materialType;			// 0: Stable Neo-Hookean, 1: Co-rotated
+		int m_maxIter;
+		Real m_maxError;
+		int m_maxIterCG;      // max CG iterations for Newton linear solve
+		Real m_tolCG;         // CG convergence tolerance
+		int m_maxLSIter;
+		Real m_lsArmijoParam;
+		Real m_lsBeta;
+		bool m_useLineSearch;
+		unsigned int m_totalNeighbors;
+		std::vector<ElasticObject*> m_objects;
+		Real m_lambda;
+		Real m_mu;
+
+#ifdef USE_AVX
+		typedef Eigen::SimplicialLLT<Eigen::SparseMatrix<double>, Eigen::Lower, Eigen::AMDOrdering<int>> SolverLLT;
+#else
+		typedef Eigen::SimplicialLLT<Eigen::SparseMatrix<double>, Eigen::Lower, Eigen::AMDOrdering<int>> SolverLLT;
+#endif
+
+		void determineFixedParticles();
+		std::string computeMD5(const unsigned int objIndex);
+		void initValues();
+		void initSystem();
+		void initFactorization(std::shared_ptr<Factorization> factorization, std::vector<unsigned int> &particleIndices, const unsigned int nFixed, const Real dt, const Real mu);
+		void findObjects();
+		void computeMatrixL();
+
+		void stepElasticitySolver();
+
+#ifndef USE_AVX
+		void computeXTilde(ElasticObject* obj);
+		void updateVelocity(ElasticObject* obj, const std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>>& xk, Real fdt);
+		Real computeEnergy(ElasticObject* obj);
+		Real computePsi(const Matrix3r& F, const Matrix3r& R) const;
+		Real computeEnergyAndGradient(ElasticObject* obj);
+		void computeHessian(ElasticObject* obj);
+		void computeCorotatedHessian9x9(ElasticObject* obj);
+		void computeStableNeoHookeanHessian9x9(ElasticObject* obj);
+		void computeNewtonPreconditioner(ElasticObject* obj);
+		void newtonMatvec(ElasticObject* obj, const std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>>& x,
+			std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>>& Ax);
+		int matFreePCG(ElasticObject* obj);
+		void prefactorizedLLTSolve(ElasticObject* obj);
+		Real newtonSolve(ElasticObject* obj, int& cgIter);
+		Real lbfgsSolve(ElasticObject* obj);
+		Real lineSearch(ElasticObject* obj, Real energy, int& lsIter);
+#endif
+
+		Matrix3r computeP(const Matrix3r& F, const Matrix3r& R) const;
+
+		virtual void initParameters();
+		/** 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.
+		*/
+		virtual void deferredInit();
+
+		//////////////////////////////////////////////////////////////////////////
+		// multiplication of symmetric matrix, represented by a 6D vector, and a
+		// 3D vector
+		//////////////////////////////////////////////////////////////////////////
+		FORCE_INLINE void symMatTimesVec(const Vector6r & M, const Vector3r & v, Vector3r &res)
+		{
+			res[0] = M[0] * v[0] + M[3] * v[1] + M[4] * v[2];
+			res[1] = M[3] * v[0] + M[1] * v[1] + M[5] * v[2];
+			res[2] = M[4] * v[0] + M[5] * v[1] + M[2] * v[2];
+		}
+
+		void rotationMatricesToAVXQuaternions();
+
+	public:
+		static std::string METHOD_NAME;
+		static int YOUNGS_MODULUS;
+		static int POISSON_RATIO;
+		static int FIXED_BOX_MIN;
+		static int FIXED_BOX_MAX;
+		static int FIXED_BOX2_MIN;
+		static int FIXED_BOX2_MAX;
+		static int ALPHA;
+		static int MAX_NEIGHBORS;
+		static int SOLVER_TYPE;
+		static int LBFGS_WINDOW_SIZE;
+		static int MATERIAL_TYPE;
+		static int MAX_ITER;
+		static int MAX_ERROR;
+		static int MAX_ITER_CG;
+		static int TOL_CG;
+		static int MAX_LS_ITER;
+		static int LS_ARMIJO_PARAM;
+		static int LS_BETA;
+		static int USE_LINE_SEARCH;
+
+		static int ENUM_SOLVER_NEWTON;
+		static int ENUM_SOLVER_LBFGS;
+		static int ENUM_MATERIAL_STABLE_NEOHOOKEAN;
+		static int ENUM_MATERIAL_COROTATED;
+
+		Elasticity_Kee2023(FluidModel *model);
+		virtual ~Elasticity_Kee2023(void);
+
+		static NonPressureForceBase* creator(FluidModel* model) { return new Elasticity_Kee2023(model); }
+		virtual std::string getMethodName() { return METHOD_NAME; }
+
+		virtual void step();
+		virtual void reset();
+		virtual void performNeighborhoodSearchSort();
+
+		virtual void saveState(BinaryFileWriter &binWriter);
+		virtual void loadState(BinaryFileReader &binReader);
+	};
+}
+
+#endif

SPlisHSPlasH/Elasticity/Elasticity_Kugelstadt2021.cpp

@@ -308,6 +308,29 @@ void Elasticity_Kugelstadt2021::initValues()
 		}
 	}
 
+	// Check neighbor list symmetry
+	{
+		unsigned int maxActualNeighbors = 0;
+		int asymmetricPairs = 0;
+		for (unsigned int i = 0; i < numParticles; i++)
+		{
+			if (m_initialNeighbors[i].size() > maxActualNeighbors)
+				maxActualNeighbors = (unsigned int)m_initialNeighbors[i].size();
+			for (size_t jn = 0; jn < m_initialNeighbors[i].size(); jn++)
+			{
+				unsigned int j = m_initialNeighbors[i][jn];
+				bool found = false;
+				for (size_t kn = 0; kn < m_initialNeighbors[j].size(); kn++)
+				{
+					if (m_initialNeighbors[j][kn] == i) { found = true; break; }
+				}
+				if (!found) asymmetricPairs++;
+			}
+		}
+		LOG_INFO << "Neighbor stats: maxNeighbors=" << maxActualNeighbors
+			<< ", asymmetricPairs=" << asymmetricPairs;
+	}
+
 	// mark all particles in the bounding box as fixed
 	determineFixedParticles();
 

SPlisHSPlasH/NonPressureForceRegistration.cpp

@@ -17,6 +17,7 @@
 #include "Elasticity/Elasticity_Becker2009.h"
 #include "Elasticity/Elasticity_Peer2018.h"
 #include "Elasticity/Elasticity_Kugelstadt2021.h"
+#include "Elasticity/Elasticity_Kee2023.h"
 
 #include "SurfaceTension/SurfaceTension_Becker2007.h"
 #include "SurfaceTension/SurfaceTension_Akinci2013.h"
@@ -39,6 +40,7 @@ void Simulation::registerNonpressureForces()
 	addElasticityMethod(Elasticity_Becker2009::METHOD_NAME, Elasticity_Becker2009::creator);
 	addElasticityMethod(Elasticity_Peer2018::METHOD_NAME, Elasticity_Peer2018::creator);
 	addElasticityMethod(Elasticity_Kugelstadt2021::METHOD_NAME, Elasticity_Kugelstadt2021::creator);
+	addElasticityMethod(Elasticity_Kee2023::METHOD_NAME, Elasticity_Kee2023::creator);
 
 	addSurfaceTensionMethod("None", [](FluidModel*) -> NonPressureForceBase* { return nullptr; });
 	addSurfaceTensionMethod(SurfaceTension_Becker2007::METHOD_NAME, SurfaceTension_Becker2007::creator);