Commit 1b108566 authored by Nghia Truong's avatar Nghia Truong
Browse files

ENH: Implement SPH fluid simulation

parent b13d5123
This diff is collapsed.
###########################################################################
#
# Copyright (c) Kitware, Inc.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0.txt
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
###########################################################################
project(Example-SPHFluid)
#-----------------------------------------------------------------------------
# Create executable
#-----------------------------------------------------------------------------
add_executable(${PROJECT_NAME} Main.cpp Fluid.cpp Solid.cpp Bunny.cpp)
#-----------------------------------------------------------------------------
# Link libraries to executable
#-----------------------------------------------------------------------------
target_link_libraries(${PROJECT_NAME} SimulationManager)
#-----------------------------------------------------------------------------
# Add shaders
#-----------------------------------------------------------------------------
include(imstkCopyAndCompileShaders)
CopyAndCompileShaders()
#-----------------------------------------------------------------------------
# Associate external data
#-----------------------------------------------------------------------------
/*=========================================================================
Library: iMSTK
Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
& Imaging in Medicine, Rensselaer Polytechnic Institute.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0.txt
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=========================================================================*/
#include "imstkSimulationManager.h"
#include "imstkSPHObject.h"
#include "imstkSPHSolver.h"
using namespace imstk;
///
/// \brief Generate a sphere-shape fluid object
///
StdVectorOfVec3d generateSphereShapeFluid(double particleRadius)
{
const double sphereRadius = 2.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
StdVectorOfVec3d particles;
particles.reserve(N * N * N);
for(size_t i = 0; i < N; ++i)
{
for(size_t j = 0; j < N; ++j)
{
for(size_t k = 0; k < N; ++k)
{
Vec3d ppos = lcorner + Vec3d(spacing * double(i), spacing * double(j), spacing * double(k));
Vec3d cx = ppos - sphereCenter;
if(cx.squaredNorm() < sphereRadiusSqr)
{
particles.push_back(ppos);
}
}
}
}
return particles;
}
///
/// \brief Generate a box-shape fluid object
///
StdVectorOfVec3d generateBoxShapeFluid(double particleRadius)
{
const double boxWidth = 4.0;
const Vec3d boxLowerCorner(-2, -3, -2);
const auto spacing = 2.0 * particleRadius;
const auto N = static_cast<size_t>(boxWidth / spacing);
StdVectorOfVec3d particles;
particles.reserve(N * N * N);
for(size_t i = 0; i < N; ++i)
{
for(size_t j = 0; j < N; ++j)
{
for(size_t k = 0; k < N; ++k)
{
Vec3d ppos = boxLowerCorner + Vec3d(spacing * double(i), spacing * double(j), spacing * double(k));
particles.push_back(ppos);
}
}
}
return particles;
}
StdVectorOfVec3d getBunny(); // Defined in Bunny.cpp
///
/// \brief Generate a bunny-shape fluid object
///
StdVectorOfVec3d generateBunnyShapeFluid(double particleRadius)
{
if(std::abs(particleRadius - 0.08) > 1e-6)
{
LOG(FATAL) << "Particle radius for this scene must be 0.08";
}
StdVectorOfVec3d particles = getBunny();
return particles;
}
std::shared_ptr<SPHObjectD> generateFluid(const std::shared_ptr<Scene>&scene, int sceneIdx, double particleRadius)
{
StdVectorOfVec3d particles;
switch(sceneIdx)
{
case 1:
particles = generateSphereShapeFluid(particleRadius);
break;
case 2:
particles = generateBoxShapeFluid(particleRadius);
break;
case 3:
particles = generateBunnyShapeFluid(particleRadius);
break;
default:
LOG(FATAL) << "Invalid scene index";
}
LOG(INFO) << "Number of particles: " << particles.size();
// Create a geometry object
auto fluidGeometry = std::make_shared<PointSet>();
fluidGeometry->initialize(particles);
// Create a fluids object
auto fluidObj = std::make_shared<SPHObjectD>("Sphere");
// Create a visiual model
auto fluidVisualModel = std::make_shared<VisualModel>(fluidGeometry);
auto fluidMaterial = std::make_shared<RenderMaterial>();
fluidMaterial->setColor(Color::Blue);
fluidMaterial->setSphereGlyphSize(particleRadius);
fluidVisualModel->setRenderMaterial(fluidMaterial);
// Create a physics model
auto sphModel = std::make_shared<SPHModelD>();
sphModel->setModelGeometry(fluidGeometry);
// configure model
auto sphParams = std::make_shared<SPHParametersD>(particleRadius);
sphParams->m_bNormalizeDensity = true;
if(sceneIdx == 2) // highly viscous fluid
{
sphParams->m_RatioKernelOverParticleRadius = 6.0;
sphParams->m_ViscosityFluid = 0.5;
sphParams->m_SurfaceTensionStiffness = 5.0;
}
sphModel->configure(sphParams);
sphModel->setTimeStepSizeType(TimeSteppingType::realTime);
// Add the component models
fluidObj->addVisualModel(fluidVisualModel);
fluidObj->setCollidingGeometry(fluidGeometry);
fluidObj->setDynamicalModel(sphModel);
fluidObj->setPhysicsGeometry(fluidGeometry); // TODO: Look into API duplication and resulting conflicts
scene->addSceneObject(fluidObj);
// Configure the solver
auto sphSolver = std::make_shared<SPHSolverD>();
sphSolver->setSPHObject(fluidObj);
scene->addNonlinearSolver(sphSolver);
return fluidObj;
}
/*=========================================================================
Library: iMSTK
Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
& Imaging in Medicine, Rensselaer Polytechnic Institute.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0.txt
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=========================================================================*/
#include "imstkSimulationManager.h"
#include "imstkSPHObject.h"
#include "imstkAPIUtilities.h"
using namespace imstk;
std::shared_ptr<SPHObjectD> generateFluid(const std::shared_ptr<Scene>&scene, int sceneIdx, double particleRadius);
std::vector<std::shared_ptr<CollidingObject>> generateSolids(const std::shared_ptr<Scene>& scene, int sceneIdx);
///
/// \brief This example demonstrates the fluid simulation using SPH
///
int main(int argc, char* argv[])
{
// SimulationManager must be created first
auto sdk = std::make_shared<SimulationManager>(0);
int threads = -1;
int sceneIdx = 1;
double particleRadius = 0.1;
// Parse command line arguments
for(int i = 1; i < argc; ++i)
{
auto param = std::string(argv[i]);
if(param.find("threads") != std::string::npos &&
param.find_first_of("=") != std::string::npos)
{
threads = std::stoi(param.substr(param.find_first_of("=") + 1));
}
else if(param.find("scene") != std::string::npos &&
param.find_first_of("=") != std::string::npos)
{
sceneIdx = std::stoi(param.substr(param.find_first_of("=") + 1));
if(sceneIdx < 1 )
{
sceneIdx = 1;
}
else if(sceneIdx > 3)
{
sceneIdx = 3;
}
LOG(INFO) << "Scene ID: " << sceneIdx;
}
else if(param.find("radius") != std::string::npos &&
param.find_first_of("=") != std::string::npos)
{
particleRadius = std::stod(param.substr(param.find_first_of("=") + 1));
LOG(INFO) << "Particle radius: " << particleRadius;
}
}
// Particle in this scene is pre-generated using particle radius 0.08
if(sceneIdx == 3)
{
particleRadius = 0.08;
}
auto scene = sdk->createNewScene("SPH Fluid");
// Generate fluid and solid objects
auto fluidObj = generateFluid(scene, sceneIdx, particleRadius);
auto solids = generateSolids(scene, sceneIdx);
// Collision between fluid and solid objects
auto colGraph = scene->getCollisionGraph();
for(auto& solid: solids)
{
if(std::dynamic_pointer_cast<Plane>(solid->getCollidingGeometry()))
{
colGraph->addInteractionPair(fluidObj, solid,
CollisionDetection::Type::PointSetToPlane,
CollisionHandling::Type::SPH,
CollisionHandling::Type::None);
}
else if(std::dynamic_pointer_cast<Sphere>(solid->getCollidingGeometry()))
{
colGraph->addInteractionPair(fluidObj, solid,
CollisionDetection::Type::PointSetToSphere,
CollisionHandling::Type::SPH,
CollisionHandling::Type::None);
}
else
{
LOG(FATAL) << "Invalid collision object";
}
}
// configure camera
scene->getCamera()->setPosition(0, 10.0, 15.0);
// configure light (white)
auto whiteLight = std::make_shared<DirectionalLight>("whiteLight");
whiteLight->setFocalPoint(Vec3d(5, -8, -5));
whiteLight->setIntensity(7);
scene->addLight(whiteLight);
// print UPS
auto ups = std::make_shared<UPSCounter>();
apiutils::printUPS(sdk->getSceneManager(scene), ups);
sdk->setActiveScene(scene);
sdk->startSimulation(SimulationStatus::PAUSED);
return 0;
}
/*=========================================================================
Library: iMSTK
Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
& Imaging in Medicine, Rensselaer Polytechnic Institute.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0.txt
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=========================================================================*/
#include "imstkSimulationManager.h"
#include "imstkPlane.h"
#include "imstkSphere.h"
using namespace imstk;
///
/// \brief Generate two planes and a solid sphere
///
std::vector<std::shared_ptr<CollidingObject>> generateSolidsScene1(const std::shared_ptr<Scene>& scene)
{
std::vector<std::shared_ptr<CollidingObject>> solids;
{
auto geometry = std::make_shared<Plane>();
geometry->setWidth(40);
geometry->setPosition(0, -6, 0);
geometry->setNormal(Vec3d(0, 1, -0.5));
auto visualModel = std::make_shared<VisualModel>(geometry);
auto material = std::make_shared<RenderMaterial>();
material->setColor(Color::DarkGray);
visualModel->setRenderMaterial(material);
auto obj = std::make_shared<CollidingObject>("Floor");
obj->addVisualModel(visualModel);
obj->setCollidingGeometry(geometry);
scene->addSceneObject(obj);
solids.push_back(obj);
}
{
auto geometry = std::make_shared<Plane>();
geometry->setWidth(40);
geometry->setPosition(0, -6, 0);
geometry->setNormal(Vec3d(0, 1, 1));
auto visualModel = std::make_shared<VisualModel>(geometry);
auto material = std::make_shared<RenderMaterial>();
material->setColor(Color::LightGray);
visualModel->setRenderMaterial(material);
auto obj = std::make_shared<CollidingObject>("Back Plane");
obj->addVisualModel(visualModel);
obj->setCollidingGeometry(geometry);
scene->addSceneObject(obj);
solids.push_back(obj);
}
{
auto geometry = std::make_shared<Sphere>();
geometry->setRadius(2);
geometry->setPosition(0, -6, 0);
auto visualModel = std::make_shared<VisualModel>(geometry);
auto material = std::make_shared<RenderMaterial>();
material->setColor(Color::Red);
visualModel->setRenderMaterial(material);
auto obj = std::make_shared<CollidingObject>("Sphere on Floor");
obj->addVisualModel(visualModel);
obj->setCollidingGeometry(geometry);
scene->addSceneObject(obj);
solids.push_back(obj);
}
return solids;
}
///
/// \brief Generate two planes
///
std::vector<std::shared_ptr<CollidingObject>> generateSolidsScene2(const std::shared_ptr<Scene>& scene)
{
std::vector<std::shared_ptr<CollidingObject>> solids;
{
auto geometry = std::make_shared<Plane>();
geometry->setWidth(40);
geometry->setPosition(0, -6, 0);
geometry->setNormal(Vec3d(0, 1, -0.5));
auto visualModel = std::make_shared<VisualModel>(geometry);
auto material = std::make_shared<RenderMaterial>();
material->setColor(Color::DarkGray);
visualModel->setRenderMaterial(material);
auto obj = std::make_shared<CollidingObject>("Floor");
obj->addVisualModel(visualModel);
obj->setCollidingGeometry(geometry);
scene->addSceneObject(obj);
solids.push_back(obj);
}
{
auto geometry = std::make_shared<Plane>();
geometry->setWidth(40);
geometry->setPosition(0, -6, 0);
geometry->setNormal(Vec3d(0, 1, 1));
auto visualModel = std::make_shared<VisualModel>(geometry);
auto material = std::make_shared<RenderMaterial>();
material->setColor(Color::LightGray);
visualModel->setRenderMaterial(material);
auto obj = std::make_shared<CollidingObject>("Back Plane");
obj->addVisualModel(visualModel);
obj->setCollidingGeometry(geometry);
scene->addSceneObject(obj);
solids.push_back(obj);
}
return solids;
}
///
/// \brief Generate an open box by 5 planes: 1 floor and 4 walls
///
std::vector<std::shared_ptr<CollidingObject>> generateSolidsScene3(const std::shared_ptr<Scene>& scene)
{
std::vector<std::shared_ptr<CollidingObject>> solids;
{
auto geometry = std::make_shared<Plane>();
geometry->setWidth(14);
geometry->setPosition(0, -6, 0);
geometry->setNormal(Vec3d(0, 1, 0));
auto visualModel = std::make_shared<VisualModel>(geometry);
auto material = std::make_shared<RenderMaterial>();
material->setColor(Color(0.2, 0.2, 0.2, 1.0));
visualModel->setRenderMaterial(material);
auto obj = std::make_shared<CollidingObject>("Floor");
obj->addVisualModel(visualModel);
obj->setCollidingGeometry(geometry);
scene->addSceneObject(obj);
solids.push_back(obj);
}
{
auto geometry = std::make_shared<Plane>();
geometry->setWidth(14);
geometry->setPosition(0, 0, -7);
geometry->setNormal(Vec3d(0, 0, 1));
auto visualModel = std::make_shared<VisualModel>(geometry);
auto material = std::make_shared<RenderMaterial>();
material->setColor(Color::LightGray);
visualModel->setRenderMaterial(material);
auto obj = std::make_shared<CollidingObject>("Back Wall");
obj->addVisualModel(visualModel);
obj->setCollidingGeometry(geometry);
scene->addSceneObject(obj);
solids.push_back(obj);
}
{
auto geometry = std::make_shared<Plane>();
geometry->setWidth(14);
geometry->setPosition(0, 0, 7);
geometry->setNormal(Vec3d(0, 0, -1));
auto visualModel = std::make_shared<VisualModel>(geometry);
auto material = std::make_shared<RenderMaterial>();
material->setColor(Color::LightGray);
visualModel->setRenderMaterial(material);
auto obj = std::make_shared<CollidingObject>("Front Wall");
obj->addVisualModel(visualModel);
obj->setCollidingGeometry(geometry);
scene->addSceneObject(obj);
solids.push_back(obj);
}
{
auto geometry = std::make_shared<Plane>();
geometry->setWidth(14);
geometry->setPosition(7, 0, 0);
geometry->setNormal(Vec3d(-1, 0, 0));
auto visualModel = std::make_shared<VisualModel>(geometry);
auto material = std::make_shared<RenderMaterial>();
material->setColor(Color::LightGray);
visualModel->setRenderMaterial(material);
auto obj = std::make_shared<CollidingObject>("Left Wall");
obj->addVisualModel(visualModel);
obj->setCollidingGeometry(geometry);
scene->addSceneObject(obj);
solids.push_back(obj);
}
{
auto geometry = std::make_shared<Plane>();
geometry->setWidth(14);
geometry->setPosition(-7, 0, 0);
geometry->setNormal(Vec3d(1, 0, 0));
auto visualModel = std::make_shared<VisualModel>(geometry);
auto material = std::make_shared<RenderMaterial>();
material->setColor(Color::LightGray);
visualModel->setRenderMaterial(material);
auto obj = std::make_shared<CollidingObject>("Right Wall");
obj->addVisualModel(visualModel);
obj->setCollidingGeometry(geometry);
scene->addSceneObject(obj);
solids.push_back(obj);
}
return solids;
}
std::vector<std::shared_ptr<CollidingObject>> generateSolids(const std::shared_ptr<Scene>& scene, int sceneIdx)
{
switch(sceneIdx)
{
case 1:
return generateSolidsScene1(scene);
case 2:
return generateSolidsScene2(scene);
case 3:
return generateSolidsScene3(scene);
default:
LOG(FATAL) << "Invalid scene index";
return {}; // To avoid warning
}
}
......@@ -45,10 +45,10 @@ PointSetToPlaneCD::computeCollisionData()
size_t nodeId = 0;
for (const auto& p : m_pointSet->getVertexPositions())
{
auto peneDistance = (planePos - p).dot(planeNormal);
auto peneDistance = (p - planePos).dot(planeNormal);
if (peneDistance <= 0.0)
{
m_colData.MAColData.push_back({ nodeId, planeNormal * -peneDistance });
m_colData.MAColData.push_back({ nodeId, planeNormal * peneDistance });
}
nodeId++;
}
......
......@@ -26,6 +26,7 @@
#include "imstkPickingCH.h"
#include "imstkDeformableObject.h"
#include "imstkBoneDrillingCH.h"
#include "imstkSPHCollisionHandling.h"