Commit 1e61004e authored by Nghia Truong's avatar Nghia Truong

ENH: Rewrite all SPH related classes, allowing to connect discrete SPH objects...

ENH: Rewrite all SPH related classes, allowing to connect discrete SPH objects with different SPH models together in the same SPH solver
parent ca1d4fb1
......@@ -21,6 +21,7 @@
#include "imstkSPHCollisionHandling.h"
#include "imstkParallelUtils.h"
#include "imstkSPHModel.h"
#include <g3log/g3log.hpp>
......@@ -37,18 +38,18 @@ SPHCollisionHandling::SPHCollisionHandling(const CollisionHandling::Side&
void
SPHCollisionHandling::processCollisionData()
{
const auto& SPHModel = m_SPHObject->getSPHModel();
const auto& sphModel = m_SPHObject->getSPHModel();
#if defined(DEBUG) || defined(_DEBUG) || !defined(NDEBUG)
LOG_IF(FATAL, (!SPHModel)) << "SPH model was not initialized";
LOG_IF(FATAL, (!sphModel)) << "SPH model was not initialized";
#endif
const auto boundaryFriction = m_SPHObject->getSPHModel()->getParameters()->m_frictionBoundary;
const auto boundaryFriction = sphModel->getParameters()->m_frictionBoundaryConstant;
#if defined(DEBUG) || defined(_DEBUG) || !defined(NDEBUG)
LOG_IF(FATAL, (boundaryFriction<0.0 || boundaryFriction>1.0))
<< "Invalid boundary friction coefficient (value must be in [0, 1])";
#endif
auto& state = SPHModel->getState();
auto& state = sphModel->getState();
ParallelUtils::parallelFor(m_colData->MAColData.getSize(),
[&](const size_t idx)
{
......
......@@ -31,16 +31,10 @@ namespace imstk
{
namespace SPH
{
template<int N>
class Poly6Kernel
{
using VecXr = Eigen::Matrix<Real, N, 1>;
public:
Poly6Kernel()
{
static_assert(N == 2 || N == 3, "Invalid kernel dimension");
}
Poly6Kernel() = default;
///
/// \brief Set the kernel radius
......@@ -49,19 +43,10 @@ public:
{
m_radius = radius;
m_radius2 = m_radius * m_radius;
if (N == 2)
{
m_k = Real(4.0) / (PI * std::pow(m_radius, 8));
m_l = -Real(24.0) / (PI * std::pow(m_radius, 8));
}
else
{
m_k = Real(315.0) / (Real(64.0) * PI * std::pow(m_radius, 9));
m_l = -Real(945.0) / (Real(32.0) * PI * std::pow(m_radius, 9));
}
m_m = m_l;
m_W0 = W(VecXr::Zero());
m_k = Real(315.0) / (Real(64.0) * PI * std::pow(m_radius, 9));
m_l = -Real(945.0) / (Real(32.0) * PI * std::pow(m_radius, 9));
m_m = m_l;
m_W0 = W(Vec3r::Zero());
}
///
......@@ -78,9 +63,9 @@ public:
/// \brief Compute weight value
/// W(r,h) = (315/(64 PI h^9))(h^2-|r|^2)^3
///
Real W(const VecXr& r) const
Real W(const Vec3r& r) const
{
const auto r2 = r.squaredNorm();
const auto r2 = r[0] * r[0] + r[1] * r[1] + r[2] * r[2];
return (r2 <= m_radius2) ? std::pow(m_radius2 - r2, 3) * m_k : Real(0);
}
......@@ -93,10 +78,10 @@ public:
/// \brief Compute weight gradient
/// grad(W(r,h)) = r(-945/(32 PI h^9))(h^2-|r|^2)^2
///
VecXr gradW(const VecXr& r) const
Vec3r gradW(const Vec3r& r) const
{
VecXr res = VecXr::Zero();
const auto r2 = r.squaredNorm();
Vec3r res = Vec3r::Zero();
const auto r2 = r[0] * r[0] + r[1] * r[1] + r[2] * r[2];
if (r2 <= m_radius2 && r2 > Real(1e-12))
{
Real tmp = m_radius2 - r2;
......@@ -110,10 +95,10 @@ public:
/// \brief Compute laplacian
/// laplacian(W(r,h)) = (-945/(32 PI h^9))(h^2-|r|^2)(-7|r|^2+3h^2)
///
Real laplacian(const VecXr& r) const
Real laplacian(const Vec3r& r) const
{
Real res = 0.;
const auto r2 = r.squaredNorm();
const auto r2 = r[0] * r[0] + r[1] * r[1] + r[2] * r[2];
if (r2 <= m_radius2)
{
Real tmp = m_radius2 - r2;
......@@ -133,16 +118,10 @@ protected:
Real m_W0; ///> Precomputed W(0)
};
template<int N>
class SpikyKernel
{
using VecXr = Eigen::Matrix<Real, N, 1>;
public:
SpikyKernel()
{
static_assert(N == 2 || N == 3, "Invalid kernel dimension");
}
SpikyKernel() = default;
///
/// \brief Set the kernel radius
......@@ -152,19 +131,10 @@ public:
m_radius = radius;
m_radius2 = m_radius * m_radius;
if (N == 2)
{
const auto radius5 = std::pow(m_radius, 5);
m_k = Real(10.0) / (PI * radius5);
m_l = -Real(30.0) / (PI * radius5);
}
else
{
const auto radius6 = std::pow(m_radius, 6);
m_k = Real(15.0) / (PI * radius6);
m_l = -Real(45.0) / (PI * radius6);
}
m_W0 = W(VecXr::Zero());
const auto radius6 = std::pow(m_radius, 6);
m_k = Real(15.0) / (PI * radius6);
m_l = -Real(45.0) / (PI * radius6);
m_W0 = W(Vec3r::Zero());
}
///
......@@ -177,9 +147,9 @@ public:
/// \brief Compute weight value
/// W(r,h) = 15/(PI*h^6) * (h-r)^3
///
Real W(const VecXr& r) const
Real W(const Vec3r& r) const
{
const auto r2 = r.squaredNorm();
const auto r2 = r[0] * r[0] + r[1] * r[1] + r[2] * r[2];
return (r2 <= m_radius2) ? std::pow(m_radius - std::sqrt(r2), 3) * m_k : Real(0);
}
......@@ -192,10 +162,10 @@ public:
/// \brief Compute weight gradient
/// grad(W(r,h)) = -r(45/(PI*h^6) * (h-r)^2)
///
VecXr gradW(const VecXr& r) const
Vec3r gradW(const Vec3r& r) const
{
VecXr res = VecXr::Zero();
const auto r2 = r.squaredNorm();
Vec3r res = Vec3r::Zero();
const auto r2 = r[0] * r[0] + r[1] * r[1] + r[2] * r[2];
if (r2 <= m_radius2 && r2 > Real(1e-12))
{
const auto rl = std::sqrt(r2);
......@@ -215,16 +185,10 @@ protected:
Real m_W0; ///> Precomputed W(0)
};
template<int N>
class CohesionKernel
{
using VecXr = Eigen::Matrix<Real, N, 1>;
public:
CohesionKernel()
{
static_assert(N == 3, "Invalid kernel dimension");
}
CohesionKernel() = default;
///
/// \brief Set the kernel radius
......@@ -234,16 +198,9 @@ public:
m_radius = radius;
m_radius2 = m_radius * m_radius;
if (N == 2)
{
LOG(FATAL) << "Unimplemented function";
}
else
{
m_k = Real(32.0) / (PI * std::pow(m_radius, 9));
m_c = std::pow(m_radius, 6) / Real(64.0);
}
m_W0 = W(VecXr::Zero());
m_k = Real(32.0) / (PI * std::pow(m_radius, 9));
m_c = std::pow(m_radius, 6) / Real(64.0);
m_W0 = W(Vec3r::Zero());
}
///
......@@ -274,10 +231,10 @@ public:
/// \brief Compute weight value
/// W(r,h) = (32/(PI h^9))(h-r)^3*r^3 if h/2 < r <= h,
/// (32/(PI h^9))(2*(h-r)^3*r^3 - h^6/64 if 0 < r <= h/2
Real W(const VecXr& r) const
Real W(const Vec3r& r) const
{
Real res = 0.;
const auto r2 = r.squaredNorm();
const auto r2 = r[0] * r[0] + r[1] * r[1] + r[2] * r[2];
if (r2 <= m_radius2)
{
const auto r1 = std::sqrt(r2);
......@@ -307,16 +264,10 @@ protected:
Real m_W0; ///> Precomputed W(0)
};
template<int N>
class AdhesionKernel
{
using VecXr = Eigen::Matrix<Real, N, 1>;
public:
AdhesionKernel()
{
static_assert(N == 3, "Invalid kernel dimension");
}
AdhesionKernel() = default;
///
/// \brief Set the kernel radius
......@@ -325,16 +276,8 @@ public:
{
m_radius = radius;
m_radius2 = m_radius * m_radius;
if (N == 2)
{
LOG(FATAL) << "Unimplemented function";
}
else
{
m_k = Real(0.007 / std::pow(m_radius, 3.25));
}
m_W0 = W(VecXr::Zero());
m_k = Real(0.007 / std::pow(m_radius, 3.25));
m_W0 = W(Vec3r::Zero());
}
///
......@@ -360,10 +303,10 @@ public:
/// \brief Compute weight value
/// W(r,h) = (0.007/h^3.25)(-4r^2/h + 6r -2h)^0.25 if h/2 < r <= h
///
Real W(const VecXr& r) const
Real W(const Vec3r& r) const
{
Real res = 0.;
const auto r2 = r.squaredNorm();
const auto r2 = r[0] * r[0] + r[1] * r[1] + r[2] * r[2];
if (r2 <= m_radius2)
{
const auto r = std::sqrt(r2);
......@@ -387,16 +330,10 @@ protected:
Real m_W0; ///> Precomputed W(0)
};
template<int N>
class ViscosityKernel
{
using VecXr = Eigen::Matrix<Real, N, 1>;
public:
ViscosityKernel()
{
static_assert(N == 2 || N == 3, "Invalid kernel dimension");
}
ViscosityKernel() = default;
///
/// \brief Set the kernel radius
......@@ -412,10 +349,10 @@ public:
/// \brief Compute laplacian
/// Laplace(r) = (45/PI/r^5) * (1 - |r| / h)
///
Real laplace(const VecXr& r) const
Real laplace(const Vec3r& r) const
{
Real res = 0.;
const auto r2 = r.squaredNorm();
const auto r2 = r[0] * r[0] + r[1] * r[1] + r[2] * r[2];
if (r2 <= m_radius2)
{
const auto d = std::sqrt(r2);
......@@ -472,9 +409,9 @@ public:
Real cohesionW(const Vec3r& r) const { return m_cohesion.W(r); }
protected:
SPH::Poly6Kernel<3> m_poly6;
SPH::SpikyKernel<3> m_spiky;
SPH::ViscosityKernel<3> m_viscosity;
SPH::CohesionKernel<3> m_cohesion;
SPH::Poly6Kernel m_poly6;
SPH::SpikyKernel m_spiky;
SPH::ViscosityKernel m_viscosity;
SPH::CohesionKernel m_cohesion;
};
} // end namespace imstk
......@@ -32,60 +32,129 @@ namespace imstk
///
/// \class SPHModelConfig
/// \brief Class that holds the SPH model parameters
/// Upon instantiating an SPH model, user must specify a particle radius
/// Other parameters are optional
///
class SPHModelConfig
{
private:
void initialize();
public:
SPHModelConfig(const Real particleRadius);
SPHModelConfig(const Real particleRadius) : m_particleRadius(particleRadius) {}
///
/// \brief Compute dependent parameters, such as particle mass, kernel radius etc.
///
void initialize();
///
/// \brief Min/Max time step sizes, used to limits the time step size during simulation
/// This is necessary to avoid integration with too large time step size which can lead to explosion,
/// or integration with too small time step size, which stalls the entire system
/// \todo Move this to solver or time integrator in the future
///
Real m_minTimestep = Real(1e-6);
Real m_maxTimestep = Real(1e-3);
Real m_CFLFactor = Real(1.0);
// particle parameters
Real m_particleRadius = Real(0);
Real m_particleRadiusSqr = Real(0); ///> \note derived quantity
///
/// \brief CFL factor, used to scale the time step size
/// (The time step size is computed based on CFL condition: dt = 2 * particle_radius / max{ |v| } *CFLFactor)
///
Real m_CFLFactor = Real(1.0);
// material parameters
Real m_restDensity = Real(1000.0);
Real m_restDensitySqr = Real(1000000.0); ///> \note derived quantity
Real m_restDensityInv = Real(1.0 / 1000.0); ///> \note derived quantity
Real m_particleMass = Real(1);
Real m_particleMassScale = Real(0.95); ///> scale particle mass to a smaller value to maintain stability
///
/// \brief Particle radius, MUST be specified by user, which is half distance between two sampled points
///
Real m_particleRadius = Real(0);
bool m_bNormalizeDensity = false;
bool m_bDensityWithBoundary = false;
///
/// \brief Rest density of the fluid (1000 kg/m^3 for water, 1.225 kg/m^3 for air)
///
Real m_restDensity = Real(1000.0);
// pressure
Real m_pressureStiffness = Real(50000.0);
///
/// \brief Scale value for particle mass, making the mass of each fluid particle slightly smaller/bigger than the ideal value
/// This is necessary if the particles were not sampled by a regular grid pattern
///
Real m_particleMassScale = Real(0.95);
// viscosity and surface tension/cohesion
Real m_viscosityCoeff = Real(1e-2);
Real m_viscosityBoundary = Real(1e-5);
Real m_surfaceTensionStiffness = Real(1);
Real m_frictionBoundary = Real(0.1);
///
/// \brief Normalize the particle densities, prodcing more accurate results
///
bool m_normalizeDensity = false;
// kernel properties
///
/// \brief This flag allows to compute fluid particle densties condidering the boundary particles (if they exist)
///
bool m_densityWithBoundary = false;
///
/// \brief Stiffness value K for computing particle pressure
///
Real m_pressureStiffnessConstant = Real(50000.0);
///
/// \brief Viscosity coefficient for computing viscosity between fluid particles
///
Real m_fluidViscosityConstant = Real(1e-2);
///
/// \brief Viscosity coefficient for computing viscosity between fluid particle and boundary particles
///
Real m_boundaryViscosityConstant = Real(1e-3);
///
/// \brief Surface tension coefficient for computing surface tension
///
Real m_surfaceTensionConstant = Real(1);
///
/// \brief Friction coefficient for computing friction between fluid particles and solid boundaries
///
Real m_frictionBoundaryConstant = Real(0.1);
///
/// \brief The ratio between SPH kernel and particles radius
/// During initialization, the SPH kernel will be computed as kernel_radius = kernel_ratio * particle_radius
///
Real m_kernelOverParticleRadiusRatio = Real(4.0);
Real m_kernelRadius; ///> \note derived quantity
Real m_kernelRadiusSqr; ///> \note derived quantity
// gravity
///
/// \brief Set gravity acceleration for the particles
///
Vec3r m_gravity = Vec3r(0, -9.81, 0);
// neighbor search
NeighborSearch::Method m_NeighborSearchMethod = NeighborSearch::Method::UniformGridBasedSearch;
///
/// \brief Method to perform neighbor search
///
NeighborSearch::Method m_neighborSearchMethod = NeighborSearch::Method::UniformGridBasedSearch;
private:
// Allow SPHModel to access private members
friend class SPHModel;
/// Squared particle radius, computed from particle radius during initialization
Real m_particleRadiusSqr = Real(0);
/// Particle mass, computed from rest density during initialization
Real m_particleMass = Real(0);
/// Squared Rest density, computed from rest density during initialization
Real m_restDensitySqr = Real(1000000.0);
/// Inversed rest density, computed from rest density during initialization
Real m_restDensityInv = Real(1.0 / 1000.0);
/// Kernel radius, computed from particle radius during initialization
Real m_kernelRadius;
/// Squared of kernel radius, computed from kernel radius during initialization
Real m_kernelRadiusSqr;
};
///
/// \class SPHModel
/// \brief SPH fluid model
///
class SPHModel : public DynamicalModel<SPHKinematicState>
class SPHModel : public DynamicalModel<SPHKinematicState>, public std::enable_shared_from_this<SPHModel>
{
public:
///
......@@ -94,12 +163,7 @@ public:
SPHModel() : DynamicalModel<SPHKinematicState>(DynamicalModelType::SPH) {}
///
/// \brief Destructor
///
virtual ~SPHModel() override = default;
///
/// \brief Set simulation parameters
/// \brief Set simulation parameters, MUST be called upon instantiating SPH model
///
void configure(const std::shared_ptr<SPHModelConfig>& params) { m_modelParameters = params; }
......@@ -178,11 +242,19 @@ public:
{ return static_cast<double>(m_dt); }
///
/// \brief Do one time step simulation
/// \brief Do one time step simulation (SPH solver should execute all steps in this function)
///
void advanceTimeStep();
private:
friend class SPHSolver;
///
/// \brief Update the bounding box of particles in all SPH models of the same SPH solver, then collect particles indices
/// This is necessary to synchronize neighbor search between SPH objects
///
void updateNeighborSearchData();
///
/// \brief Compute time step size, do nothing if using a fixed time step size for integration
///
......@@ -199,7 +271,8 @@ private:
void findParticleNeighbors();
///
/// \brief Pre-compute relative positions with neighbor particles
/// \brief Pre-compute relative positions x_{ij} = x_i - x_j with all neighbor particles j
/// x_{ij} are used multiples times, thus caching them can improve performance significantly
///
void computeNeighborRelativePositions();
......@@ -220,40 +293,80 @@ private:
void normalizeDensity();
///
/// \brief Compute particle acceleration due to pressure
/// \brief Compute particle accelerations due to pressure
///
void computePressureAcceleration();
///
/// \brief Update particle velocities due to pressure
/// \brief Compute surface normals
///
void updateVelocity(const Real timestep);
void computeSurfaceNormal();
///
/// \brief Compute viscosity
/// \brief Compute surface tension using Akinci et at. 2013 model
/// (Versatile Surface Tension and Adhesion for SPH Fluids)
///
void computeViscosity();
void computeSurfaceTensionAcceleration();
///
/// \brief Compute surface tension and update velocities
/// Compute surface tension using Akinci et at. 2013 model
/// (Versatile Surface Tension and Adhesion for SPH Fluids)
/// \brief Update particle velocities from acceleration
///
void computeSurfaceTension();
void updateVelocity(const Real timestep);
///
/// \brief Compute viscosity (in the XSPH viscosity model, visocisty is computed after velocity update)
///
void computeViscosity();
///
/// \brief Move particles
///
void moveParticles(const Real timestep);
std::shared_ptr<PointSet> m_geometry;
SPHSimulationState m_simulationState;
/// A set containing all SPH models belonging to the same SPH solver
using SPHModelGroup = std::vector<std::shared_ptr<SPHModel>>;
///
/// \brief Initialize group of SPH models for a newly created SPH model
/// Each SPH model must associate with a model group, which consists of at least itself
/// Thus, this function MUST be called during SPH model initialization
///
static void initializeModelGroup(const std::shared_ptr<SPHModel>& sphModel);
///
/// \brief Return the set of SPH models which are in the same SPH solver
///
static const SPHModelGroup& getModelGroup(SPHModel* const sphModel);
///
/// \brief Setup model group for all SPH models in the same SPH solver
///
static void setupModelGroup(const std::vector<std::shared_ptr<SPHModel>>& models);
/// Map from a pointer of an SPH model to the set of its mode group
static std::unordered_map<SPHModel*, SPHModelGroup> s_mSPHModelGroups;
/// Turn this flag on/off to know whether an SPH model has been added to an SPH solver or not
/// This is used to check to avoid adding the same SPH model to multiple SPH solver
/// Thus, it should be access only in SPH solver during adding/removing model
bool m_connectedToSolver = false;
std::shared_ptr<PointSet> m_geometry; ///> Original pointset geometry
SPHSimulationState m_simulationState; ///> The simulation state of the SPH model
/// Time step size, which is computed from particle velocities at every time step, if using CFL condition
Real m_dt;
/// If using a fixed time step for simulation, the time step size will be fixed by this value
Real m_defaultDt = Real(1e-4);
/// The SPH kernels, which are initialized during model initialization
SPHSimulationKernels m_kernels;
Real m_dt; ///> time step size
Real m_defaultDt; ///> default time step size
/// Simulation parameters, which MUST be set upon instantiating an SPH model
std::shared_ptr<SPHModelConfig> m_modelParameters;
SPHSimulationKernels m_kernels; ///> SPH kernels (must be initialized during model initialization)
std::shared_ptr<SPHModelConfig> m_modelParameters; ///> SPH Model parameters (must be set before simulation)
std::shared_ptr<NeighborSearch> m_neighborSearcher; ///> Neighbor Search (must be initialized during model initialization)
/// Neighbor search data structure, which is initialized during model initialization
std::shared_ptr<NeighborSearch> m_neighborSearcher;
};
} // end namespace imstk