.. _program_listing_file_SPlisHSPlasH_Elasticity_Elasticity_Peer2018.h:

Program Listing for File Elasticity_Peer2018.h
==============================================

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

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

.. code-block:: cpp

   #ifndef __Elasticity_Peer2018_h__
   #define __Elasticity_Peer2018_h__
   
   #include "SPlisHSPlasH/Common.h"
   #include "SPlisHSPlasH/FluidModel.h"
   #include "SPlisHSPlasH/NonPressureForceBase.h"
   #include "SPlisHSPlasH/Utilities/MatrixFreeSolver.h"
   
   namespace SPH
   {
       class Elasticity_Peer2018 : public NonPressureForceBase
       {
       protected:
           typedef Eigen::ConjugateGradient<MatrixReplacement, Eigen::Lower | Eigen::Upper, Eigen::IdentityPreconditioner> Solver;
   
           Real m_youngsModulus;
           Real m_poissonRatio;
           Vector3r m_fixedBoxMin;
           Vector3r m_fixedBoxMax;
   
           // 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<Matrix3r> m_stress;
           std::vector<Matrix3r> m_L;
           std::vector<Matrix3r> m_RL;
           std::vector<Matrix3r> m_F;
           unsigned int m_iterations;
           unsigned int m_maxIter;
           Real m_maxError;
           Real m_alpha;
           int m_maxNeighbors;
           Solver m_solver;
   
           void determineFixedParticles();
           void initValues();
           void computeMatrixL();
           void computeRotations();
           void computeRHS(VectorXr & rhs);    
   
           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];
           }
   
           FORCE_INLINE void generateIndices(const unsigned int* map, const unsigned int* idx, std::array<unsigned int, 8>& indices, const unsigned char count = 8u)
           {
               switch (count)
               {
               case 1u:
                   indices = { map[idx[0]], 0, 0, 0, 0, 0, 0, 0 }; break;
               case 2u:
                   indices = { map[idx[0]], map[idx[1]], 0, 0, 0, 0, 0, 0 }; break;
               case 3u:
                   indices = { map[idx[0]], map[idx[1]], map[idx[2]], 0, 0, 0, 0, 0 }; break;
               case 4u:
                   indices = { map[idx[0]], map[idx[1]], map[idx[2]], map[idx[3]], 0, 0, 0, 0 }; break;
               case 5u:
                   indices = { map[idx[0]], map[idx[1]], map[idx[2]], map[idx[3]], map[idx[4]], 0, 0, 0 }; break;
               case 6u:
                   indices = { map[idx[0]], map[idx[1]], map[idx[2]], map[idx[3]], map[idx[4]], map[idx[5]], 0, 0 }; break;
               case 7u:
                   indices = { map[idx[0]], map[idx[1]], map[idx[2]], map[idx[3]], map[idx[4]], map[idx[5]], map[idx[6]], 0 }; break;
               case 8u:
                   indices = { map[idx[0]], map[idx[1]], map[idx[2]], map[idx[3]], map[idx[4]], map[idx[5]], map[idx[6]], map[idx[7]] }; break;
               }
           }
   
   
       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 ITERATIONS;
           static int MAX_ITERATIONS;
           static int MAX_ERROR;
           static int ALPHA;
           static int MAX_NEIGHBORS;
   
           Elasticity_Peer2018(FluidModel *model);
           virtual ~Elasticity_Peer2018(void);
   
           static NonPressureForceBase* creator(FluidModel* model) { return new Elasticity_Peer2018(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);
   
           static void matrixVecProd(const Real* vec, Real *result, void *userData);
       };
   }
   
   #endif