.. _program_listing_file_SPlisHSPlasH_Elasticity_Elasticity_Kee2023.h:

Program Listing for File Elasticity_Kee2023.h
=============================================

|exhale_lsh| :ref:`Return to documentation for file <file_SPlisHSPlasH_Elasticity_Elasticity_Kee2023.h>` (``SPlisHSPlasH/Elasticity/Elasticity_Kee2023.h``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #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
   {
       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;
   
               // Newton: per-particle 9×9 Hessian and block-diagonal preconditioner.
               // Same type in both builds (SVD is always scalar).
               std::vector<Eigen::Matrix<Real, 9, 9>> m_hessian9x9;
               std::vector<Matrix3r, Eigen::aligned_allocator<Matrix3r>> m_pcg_precond;
   
   #ifdef USE_AVX
               // AVX state (Scalarf8, 3 active lanes x/y/z).
               std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_f_avx;          // F = D·xk workspace (3*numParticles)
               std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_xk_avx;         // current iterate
               std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_xTilde_avx;     // inertial target
               std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_dx_avx;         // solver step (also LLT solve RHS/result)
               std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_gradient_avx;   // ∇E at xk
   
               // Newton CG workspace (Scalarf8, coord-packed x/y/z in lanes 0/1/2)
               std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_pcg_r_avx;
               std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_pcg_p_avx;
               std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_pcg_Ap_avx;
               std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_pcg_z_avx;
   
               // L-BFGS secant history
               std::vector<std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>>> m_lbfgs_s_avx;
               std::vector<std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>>> m_lbfgs_y_avx;
               std::vector<Real> m_lbfgs_rho;
               std::vector<Real> m_lbfgs_alpha;
               std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_lbfgs_last_sol_avx;
               std::vector<Scalarf8, AlignmentAllocator<Scalarf8, 32>> m_lbfgs_q_avx;
               int m_lbfgs_count = 0;
   #else
               std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_f;        // F = D·xk workspace
               std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_xk;       // current iterate
               std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_xTilde;   // inertial target
               std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_dx;       // Newton/LBFGS step (also LLT solve RHS/result)
               std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_gradient; // ∇E at xk
   
               // Newton PCG workspace
               std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_pcg_r;
               std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_pcg_p;
               std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_pcg_Ap;
               std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_pcg_z;
   
               // L-BFGS secant history
               std::vector<std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>>> m_lbfgs_s;
               std::vector<std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>>> m_lbfgs_y;
               std::vector<Real> m_lbfgs_rho;
               std::vector<Real> m_lbfgs_alpha;
               std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_lbfgs_last_sol;
               std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_lbfgs_q;
               int m_lbfgs_count = 0;
   
               // Permutation workspace for scalar LLT (manual forward/backward sub)
               std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_dx_perm;
   #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;
   
           // Precomputed V_j * gradW(xi0 - xj0) per neighbor pair.
           // Scalar format (always available, used by Newton CG).
           std::vector<Vector3r, Eigen::aligned_allocator<Vector3r>> m_precomp_V_gradW;
           std::vector<unsigned int> m_precomputed_indices;
   #ifdef USE_AVX
           // AVX-packed format (8 neighbors per entry, used by L-BFGS force loops).
           std::vector<Vector3f8, Eigen::aligned_allocator<Vector3f8>> m_precomp_V_gradW8;
           std::vector<unsigned int> m_precomputed_indices8;
   #endif
   
   #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 precomputeValues();
   
           void stepElasticitySolver();
   
           void computeXTilde(ElasticObject* obj);
           void updateVelocity(ElasticObject* obj, 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);
           int matFreePCG(ElasticObject* obj);
           Real newtonSolve(ElasticObject* obj, int& cgIter);
           void prefactorizedLLTSolve(ElasticObject* obj);
           Real lbfgsSolve(ElasticObject* obj);
           Real lineSearch(ElasticObject* obj, Real energy, int& lsIter);
   
           Matrix3r computeP(const Matrix3r& F, const Matrix3r& R) const;
   
           virtual void initParameters();
           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