imstkSphModel.h

/*
** This file is part of the Interactive Medical Simulation Toolkit (iMSTK)
** iMSTK is distributed under the Apache License, Version 2.0.
** See accompanying NOTICE for details.
*/

#pragma once

#include "imstkDynamicalModel.h"
#include "imstkSphState.h"
#include "imstkSPHKernels.h"
#include "imstkNeighborSearch.h"
#include "imstkSphBoundaryConditions.h"

namespace imstk
{
class PointSet;

class SphModelConfig
{
private:
    void initialize();

public:
    SphModelConfig(const double particleRadius);
    SphModelConfig(const double particleRadius, const double speedOfSound, const double restDensity);

    double m_minTimestep = 1.0e-6;
    double m_maxTimestep = 1.0e-3;
    double m_cflFactor   = 1.0;

    // particle parameters
    double m_particleRadius    = 0.0;
    double m_particleRadiusSqr = 0.0; 

    // material parameters
    double m_restDensity       = 1000.0;
    double m_restDensitySqr    = 1000000.0;    
    double m_restDensityInv    = 1.0 / 1000.0; 
    double m_particleMass      = 1.0;
    double m_particleMassScale = 1.0;          
    double m_eta = 0.5;                        

    bool m_bNormalizeDensity    = false;
    bool m_bDensityWithBoundary = false;

    // pressure
    double m_pressureStiffness = 50000.0;

    // viscosity and surface tension/cohesion
    double m_dynamicViscosityCoeff   = 1.0e-2;
    double m_viscosityBoundary       = 1.0e-5;
    double m_surfaceTensionStiffness = 1.0;
    double m_frictionBoundary = 0.1;

    // kernel properties
    double m_kernelOverParticleRadiusRatio = 4.0;
    double m_kernelRadius;    
    double m_kernelRadiusSqr; 

    // gravity
    Vec3d m_gravity = Vec3d(0.0, -9.81, 0.0);

    // sound speed
    double m_speedOfSound = 18.7;

    // neighbor search
    NeighborSearch::Method m_neighborSearchMethod = NeighborSearch::Method::UniformGridBasedSearch;
};

class SphModel : public DynamicalModel<SphState>
{
public:
    SphModel();
    ~SphModel() override = default;

    void configure(const std::shared_ptr<SphModelConfig>& params) { m_modelParameters = params; }

    bool initialize() override;

    void resetToInitialState() override { this->m_currentState->setState(this->m_initialState); }

    const std::shared_ptr<SphModelConfig>& getParameters() const
    {
        assert(m_modelParameters);
        return m_modelParameters;
    }

    void setTimeStep(const double timeStep) override { setDefaultTimeStep(timeStep); }

    void setDefaultTimeStep(const double timeStep) { m_defaultDt = timeStep; }

    double getTimeStep() const override { return m_dt; }

    void setInitialVelocities(const size_t numParticles, const Vec3d& initialVelocities);

    double getParticlePressure(const double density);

    void findNearestParticleToVertex(const VecDataArray<double, 3>& points, const std::vector<std::vector<size_t>>& indices);

    void setBoundaryConditions(std::shared_ptr<SphBoundaryConditions> sphBoundaryConditions) { m_sphBoundaryConditions = sphBoundaryConditions; }
    std::shared_ptr<SphBoundaryConditions> getBoundaryConditions() { return m_sphBoundaryConditions; }

    void setRestDensity(const double restDensity) { m_modelParameters->m_restDensity = restDensity; }

    std::shared_ptr<TaskNode> getFindParticleNeighborsNode() const { return m_findParticleNeighborsNode; }
    std::shared_ptr<TaskNode> getComputeDensityNode() const { return m_computeDensityNode; }
    std::shared_ptr<TaskNode> getComputePressureNode() const { return m_computePressureAccelNode; }
    std::shared_ptr<TaskNode> getComputeSurfaceTensionNode() const { return m_computeSurfaceTensionNode; }
    std::shared_ptr<TaskNode> getComputeTimeStepSizeNode() const { return m_computeTimeStepSizeNode; }
    std::shared_ptr<TaskNode> getSumAccelsNode() const { return m_sumAccelsNode; }
    std::shared_ptr<TaskNode> getIntegrateNode() const { return m_integrateNode; }
    std::shared_ptr<TaskNode> getComputeViscosityNode() const { return m_computeViscosityNode; }
    std::shared_ptr<TaskNode> getUpdateVelocityNode() const { return m_updateVelocityNode; }
    std::shared_ptr<TaskNode> getMoveParticlesNode() const { return m_moveParticlesNode; }

protected:
    void initGraphEdges(std::shared_ptr<TaskNode> source, std::shared_ptr<TaskNode> sink) override;

private:
    void computeTimeStepSize();

    double computeCFLTimeStepSize();

    void findParticleNeighbors();

    void computeNeighborRelativePositions();

    void collectNeighborDensity();

    void computeDensity();

    void normalizeDensity();

    void computePressureAcceleration();

    void sumAccels();

    void updateVelocity(const double timestep);

    void computeViscosity();

    void computeSurfaceTension();

    void moveParticles(const double timestep);

//void computePressureOutlet();

protected:
    std::shared_ptr<TaskNode> m_findParticleNeighborsNode = nullptr;
    std::shared_ptr<TaskNode> m_computeDensityNode        = nullptr;
    std::shared_ptr<TaskNode> m_computePressureAccelNode  = nullptr;
    std::shared_ptr<TaskNode> m_computeSurfaceTensionNode = nullptr;
    std::shared_ptr<TaskNode> m_computeTimeStepSizeNode   = nullptr;
    std::shared_ptr<TaskNode> m_sumAccelsNode              = nullptr;
    std::shared_ptr<TaskNode> m_integrateNode              = nullptr;
    std::shared_ptr<TaskNode> m_updateVelocityNode         = nullptr;
    std::shared_ptr<TaskNode> m_computeViscosityNode       = nullptr;
    std::shared_ptr<TaskNode> m_moveParticlesNode          = nullptr;
    std::shared_ptr<TaskNode> m_normalizeDensityNode       = nullptr;
    std::shared_ptr<TaskNode> m_collectNeighborDensityNode = nullptr;

private:
    std::shared_ptr<PointSet> m_pointSetGeometry;

    double m_dt = 0.0;                                  
    double m_defaultDt;                                 

    SphSimulationKernels m_kernels;                     
    std::shared_ptr<SphModelConfig> m_modelParameters;  
    std::shared_ptr<NeighborSearch> m_neighborSearcher; 

    std::shared_ptr<VecDataArray<double, 3>> m_pressureAccels       = nullptr;
    std::shared_ptr<VecDataArray<double, 3>> m_surfaceTensionAccels = nullptr;
    std::shared_ptr<VecDataArray<double, 3>> m_viscousAccels    = nullptr;
    std::shared_ptr<VecDataArray<double, 3>> m_neighborVelContr = nullptr;
    std::shared_ptr<VecDataArray<double, 3>> m_particleShift    = nullptr;

    std::shared_ptr<VecDataArray<double, 3>> m_initialVelocities = nullptr;
    std::shared_ptr<DataArray<double>>       m_initialDensities  = nullptr;

    int m_timeStepCount = 0;

    std::shared_ptr<SphBoundaryConditions> m_sphBoundaryConditions = nullptr;

    std::vector<size_t> m_minIndices;
};
} // namespace imstk