Program Listing for File Elasticity_Peer2018.h

↰ Return to documentation for file (SPlisHSPlasH/Elasticity/Elasticity_Peer2018.h)

#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