Commit d3f7e50b authored by Sreekanth Arikatla's avatar Sreekanth Arikatla

Merge branch 'FixSPHBoundaryCondition' into 'master'

Fix SPH Boundary Condition

See merge request !343
parents e4db0984 1590733a
Pipeline #141336 passed with stage
in 0 seconds
......@@ -28,15 +28,16 @@ using namespace imstk;
///
/// \brief Generate a sphere-shape fluid object
///
StdVectorOfVec3d generateSphereShapeFluid(const double particleRadius)
StdVectorOfVec3d
generateSphereShapeFluid(const double particleRadius)
{
const double sphereRadius = 2.0;
const Vec3d sphereCenter(0, 1, 0);
const Vec3d sphereCenter(0, 1, 0);
const auto sphereRadiusSqr = sphereRadius * sphereRadius;
const auto spacing = 2.0 * particleRadius;
const auto N = static_cast<size_t>(2.0 * sphereRadius / spacing); // Maximum number of particles in each dimension
const Vec3d lcorner = sphereCenter - Vec3d(sphereRadius, sphereRadius, sphereRadius); // Cannot use auto here, due to Eigen bug
const auto sphereRadiusSqr = sphereRadius * sphereRadius;
const auto spacing = 2.0 * particleRadius;
const auto N = static_cast<size_t>(2.0 * sphereRadius / spacing); // Maximum number of particles in each dimension
const Vec3d lcorner = sphereCenter - Vec3d(sphereRadius, sphereRadius, sphereRadius); // Cannot use auto here, due to Eigen bug
StdVectorOfVec3d particles;
particles.reserve(N * N * N);
......@@ -63,10 +64,11 @@ StdVectorOfVec3d generateSphereShapeFluid(const double particleRadius)
///
/// \brief Generate a box-shape fluid object
///
StdVectorOfVec3d generateBoxShapeFluid(const double particleRadius)
StdVectorOfVec3d
generateBoxShapeFluid(const double particleRadius)
{
const double boxWidth = 4.0;
const Vec3d boxLowerCorner(-2, -3, -2);
const Vec3d boxLowerCorner(-2, -3, -2);
const auto spacing = 2.0 * particleRadius;
const auto N = static_cast<size_t>(boxWidth / spacing);
......@@ -95,7 +97,8 @@ StdVectorOfVec3d getBunny(); // Defined in Bunny.cpp
///
/// \brief Generate a bunny-shape fluid object
///
StdVectorOfVec3d generateBunnyShapeFluid(const double particleRadius)
StdVectorOfVec3d
generateBunnyShapeFluid(const double particleRadius)
{
LOG_IF(FATAL, (std::abs(particleRadius - 0.08) > 1e-6)) << "Particle radius for this scene must be 0.08";
StdVectorOfVec3d particles;
......@@ -105,7 +108,8 @@ StdVectorOfVec3d generateBunnyShapeFluid(const double particleRadius)
return particles;
}
std::shared_ptr<SPHObject> generateFluid(const std::shared_ptr<Scene>&scene, const double particleRadius)
std::shared_ptr<SPHObject>
generateFluid(const std::shared_ptr<Scene>& scene, const double particleRadius)
{
StdVectorOfVec3d particles;
switch (SCENE_ID)
......@@ -134,7 +138,7 @@ std::shared_ptr<SPHObject> generateFluid(const std::shared_ptr<Scene>&scene, con
// Create a visiual model
auto fluidVisualModel = std::make_shared<VisualModel>(fluidGeometry);
auto fluidMaterial = std::make_shared<RenderMaterial>();
auto fluidMaterial = std::make_shared<RenderMaterial>();
fluidMaterial->setColor(Color::Blue);
fluidMaterial->setSphereGlyphSize(particleRadius);
fluidVisualModel->setRenderMaterial(fluidMaterial);
......@@ -149,10 +153,15 @@ std::shared_ptr<SPHObject> generateFluid(const std::shared_ptr<Scene>&scene, con
if (SCENE_ID == 2) // highly viscous fluid
{
sphParams->m_RatioKernelOverParticleRadius = 6.0;
sphParams->m_ViscosityFluid = 0.5;
sphParams->m_ViscosityFluid = 0.5;
sphParams->m_SurfaceTensionStiffness = 5.0;
}
if (SCENE_ID == 3) // bunny-shaped fluid
{
sphParams->m_BoundaryFriction = 0.3;
}
sphModel->configure(sphParams);
sphModel->setTimeStepSizeType(TimeSteppingType::realTime);
......
......@@ -34,35 +34,21 @@ SPHCollisionHandling::SPHCollisionHandling(const CollisionHandling::Side&
LOG_IF(FATAL, (!m_SPHObject)) << "Invalid SPH object";
}
void
SPHCollisionHandling::setBoundaryFriction(const Real friction)
{
m_BoundaryFriction = friction;
if (m_BoundaryFriction < 0.0 || m_BoundaryFriction > 1.0)
{
LOG(WARNING) << "Invalid frictionLength coefficient (it must be in [0, 1])";
if (m_BoundaryFriction < 0)
{
m_BoundaryFriction = 0;
}
else if (m_BoundaryFriction > static_cast<Real>(1.0))
{
m_BoundaryFriction = static_cast<Real>(1.0);
}
}
}
void
SPHCollisionHandling::processCollisionData()
{
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
auto& state = SPHModel->getState();
const auto boundaryFriction = m_SPHObject->getSPHModel()->getParameters()->m_BoundaryFriction;
#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();
ParallelUtils::parallelFor(m_colData->MAColData.size(),
[&](const size_t idx)
{
......@@ -70,33 +56,40 @@ SPHCollisionHandling::processCollisionData()
const auto pidx = cd.nodeId; // Fluid particle index
auto n = cd.penetrationVector; // This vector should point into solid object
// Correct particle position:
state.getPositions()[pidx] -= n; // This should not cause race condition, since pidx is distinct
// Correct particle position
state.getPositions()[pidx] -= n;
// Correct particle velocity: slip boundary condition with friction
auto& vel = state.getVelocities()[pidx];
const auto nLengthSqr = n.squaredNorm();
if (nLengthSqr < Real(1e-20)) // Normalize n
{
return; // Too little penetration: ignore
}
n /= std::sqrt(nLengthSqr);
const auto vn = vel.dot(n);
vel -= n * vn; // From now, vel is parallel with the solid surface
if (m_BoundaryFriction > Real(1e-20))
// Correct particle velocity: slip boundary condition with friction
const auto oldVel = state.getVelocities()[pidx];
const auto vn = oldVel.dot(n);
// If particle is escaping the boundary, ignore it
if (vn > 0)
{
const auto velLength = vel.norm();
const auto frictionLength = -vn * m_BoundaryFriction;
if (frictionLength < velLength && velLength > Real(1e-10))
{
vel -= (vel / velLength) * frictionLength; // Subtract a friction from velocity, which is proportional to the amount of penetration
}
else
Vec3r correctedVel = oldVel - vn * n; // From now, vel is parallel with the solid surface
if (boundaryFriction > Real(1e-20))
{
vel = Vec3r::Zero();
const auto velLength = correctedVel.norm();
const auto frictionLength = vn * boundaryFriction; // This is always positive
if (frictionLength < velLength && velLength > Real(1e-10))
{
correctedVel -= (correctedVel / velLength) * frictionLength; // Subtract a friction from velocity, which is proportional to the amount of penetration
}
else
{
correctedVel = Vec3r::Zero();
}
}
state.getVelocities()[pidx] = correctedVel;
}
});
}
......
......@@ -33,25 +33,12 @@ public:
SPHCollisionHandling() = delete;
virtual ~SPHCollisionHandling() override = default;
///
/// \brief Compute forces based on collision data
///
virtual void processCollisionData() override;
///
/// \brief Set the friction coefficient
///
void setBoundaryFriction(const Real friction);
///
/// \brief Get the friction coefficient
///
Real getBoundaryFriction() const { return m_BoundaryFriction; }
private:
std::shared_ptr<SPHObject> m_SPHObject;
Real m_BoundaryFriction = Real(0.1);
};
} // end namespace imstk
......@@ -76,6 +76,7 @@ public:
Real m_ViscosityFluid = Real(1e-2);
Real m_ViscosityBoundary = Real(1e-5);
Real m_SurfaceTensionStiffness = Real(1);
Real m_BoundaryFriction = Real(0.1);
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment