Commit 6e2eb4b6 authored by Nghia Truong's avatar Nghia Truong
Browse files

REFACT: Change logging and const correctness

parent 5fe9ea3a
......@@ -28,7 +28,7 @@ using namespace imstk;
///
/// \brief Generate a sphere-shape fluid object
///
StdVectorOfVec3d generateSphereShapeFluid(double particleRadius)
StdVectorOfVec3d generateSphereShapeFluid(const double particleRadius)
{
const double sphereRadius = 2.0;
const Vec3d sphereCenter(0, 1, 0);
......@@ -63,7 +63,7 @@ StdVectorOfVec3d generateSphereShapeFluid(double particleRadius)
///
/// \brief Generate a box-shape fluid object
///
StdVectorOfVec3d generateBoxShapeFluid(double particleRadius)
StdVectorOfVec3d generateBoxShapeFluid(const double particleRadius)
{
const double boxWidth = 4.0;
const Vec3d boxLowerCorner(-2, -3, -2);
......@@ -93,19 +93,15 @@ StdVectorOfVec3d getBunny(); // Defined in Bunny.cpp
///
/// \brief Generate a bunny-shape fluid object
///
StdVectorOfVec3d generateBunnyShapeFluid(double particleRadius)
StdVectorOfVec3d generateBunnyShapeFluid(const double particleRadius)
{
if(std::abs(particleRadius - 0.08) > 1e-6)
{
LOG(FATAL) << "Particle radius for this scene must be 0.08";
}
LOG_IF(FATAL, (std::abs(particleRadius - 0.08) > 1e-6)) << "Particle radius for this scene must be 0.08";
StdVectorOfVec3d particles = getBunny();
return particles;
}
std::shared_ptr<SPHObject> generateFluid(const std::shared_ptr<Scene>&scene, int sceneIdx, double particleRadius)
std::shared_ptr<SPHObject> generateFluid(const std::shared_ptr<Scene>&scene, int sceneIdx, const double particleRadius)
{
StdVectorOfVec3d particles;
switch(sceneIdx)
......
......@@ -26,11 +26,13 @@
using namespace imstk;
std::shared_ptr<SPHObject> generateFluid(const std::shared_ptr<Scene>&scene, int sceneIdx, double particleRadius);
std::shared_ptr<SPHObject> generateFluid(const std::shared_ptr<Scene>&scene, int sceneIdx, const 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
/// \brief Usage: ./SPHFluid [scene=<scene_id>] [threads=<num_threads>] [radius=<particle_radius>], where scene_id is in [1..3]
/// \brief Example: ./SPHFluid scene=1 threads=8 radius=0.01
///
int main(int argc, char* argv[])
{
......@@ -117,10 +119,6 @@ int main(int argc, char* argv[])
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);
......
......@@ -23,29 +23,32 @@
namespace imstk
{
void SPHCollisionHandling::setBoundaryFriction(Real frictionLength)
void SPHCollisionHandling::setBoundaryFriction(const Real friction)
{
if(frictionLength < 0.0 || frictionLength > 1.0)
m_BoundaryFriction = friction;
if(m_BoundaryFriction < 0.0 || m_BoundaryFriction > 1.0)
{
LOG(WARNING) << "Invalid frictionLength coefficient (it must be in [0, 1])";
if(frictionLength < 0)
if(m_BoundaryFriction < 0)
{
frictionLength = 0;
m_BoundaryFriction = 0;
}
else if(frictionLength > 1.0)
else if(m_BoundaryFriction > static_cast<Real>(1.0))
{
frictionLength = 1.0;
m_BoundaryFriction = static_cast<Real>(1.0);
}
}
m_BoundaryFriction = frictionLength;
}
void SPHCollisionHandling::computeContactForces()
{
const auto& SPHModel = m_SPHObject->getSPHModel();
assert(SPHModel);
auto& state = SPHModel->getState();
#if defined(DEBUG) || defined(_DEBUG) || !defined(NDEBUG)
LOG_IF(FATAL, !SPHModel) << "SPH model was not initialized";
#endif
auto& state = SPHModel->getState();
for(const auto& cd : m_colData.MAColData)
{
const auto pidx = cd.nodeId; // Fluid particle index
......
......@@ -42,7 +42,7 @@ public:
///
/// \brief Set the friction coefficient
///
void setBoundaryFriction(Real friction);
void setBoundaryFriction(const Real friction);
///
/// \brief Get the friction coefficient
......@@ -53,4 +53,4 @@ private:
std::shared_ptr<SPHObject> m_SPHObject;
Real m_BoundaryFriction = 0.1;
};
} // end namespace imstk
\ No newline at end of file
} // end namespace imstk
......@@ -68,8 +68,8 @@ public:
void getNeighbors(std::vector<std::vector<size_t>>& result, const StdVectorOfVec3r& setA, const StdVectorOfVec3r& setB);
private:
Real m_SearchRadius = Real(0);
Real m_SearchRadiusSqr = Real(0);
Real m_SearchRadius = 0.;
Real m_SearchRadiusSqr = 0.;
UniformSpatialGrid<std::vector<size_t>> m_Grid;
};
} // end namespace imstk
......@@ -49,7 +49,8 @@ public:
///
void initialize(const Vec3r& lowerCorner, const Vec3r& upperCorner, const Real cellSize)
{
assert(cellSize > 0);
LOG_IF(FATAL, (cellSize <= 0)) << "Invalid cell size";
m_LowerCorner = lowerCorner;
m_UpperCorner = upperCorner;
......@@ -63,13 +64,10 @@ public:
m_NTotalCells *= m_Resolution[i];
}
if(m_NTotalCells == 0)
{
LOG(FATAL) << "Invalid grid size: [" +
std::to_string(m_LowerCorner[0]) + ", " + std::to_string(m_LowerCorner[1]) + ", " + std::to_string(m_LowerCorner[2]) + "] => " +
std::to_string(m_UpperCorner[0]) + ", " + std::to_string(m_UpperCorner[1]) + ", " + std::to_string(m_UpperCorner[2]) + "], " +
"cellSize = " + std::to_string(m_CellSize);
}
LOG_IF(FATAL, (m_NTotalCells == 0)) << "Invalid grid size: [" +
std::to_string(m_LowerCorner[0]) + ", " + std::to_string(m_LowerCorner[1]) + ", " + std::to_string(m_LowerCorner[2]) + "] => " +
std::to_string(m_UpperCorner[0]) + ", " + std::to_string(m_UpperCorner[1]) + ", " + std::to_string(m_UpperCorner[2]) + "], " +
"cellSize = " + std::to_string(m_CellSize);
// cell data must be resized to equal to number of cells
m_CellData.resize(m_NTotalCells);
......@@ -175,13 +173,11 @@ private:
{
auto cellIdx = getCellIndexFromCoordinate<IndexType>(ppos);
#if defined(DEBUG) || defined(_DEBUG) || !defined(NDEBUG)
if(!isValidCellIndices(cellIdx[0], cellIdx[1], cellIdx[2]))
{
LOG(FATAL) << "Invalid cell indices: " +
std::to_string(cellIdx[0]) + "/" + std::to_string(m_Resolution[0]) + ", " +
std::to_string(cellIdx[1]) + "/" + std::to_string(m_Resolution[1]) + ", " +
std::to_string(cellIdx[2]) + "/" + std::to_string(m_Resolution[2]);
}
LOG_IF(FATAL, !isValidCellIndices(cellIdx[0], cellIdx[1], cellIdx[2])) <<
"Invalid cell indices: " +
std::to_string(cellIdx[0]) + "/" + std::to_string(m_Resolution[0]) + ", " +
std::to_string(cellIdx[1]) + "/" + std::to_string(m_Resolution[1]) + ", " +
std::to_string(cellIdx[2]) + "/" + std::to_string(m_Resolution[2]);
#endif
return getCellFlatIndexFrom3DIndices<IndexType>(cellIdx[0], cellIdx[1], cellIdx[2]);
}
......
......@@ -23,6 +23,7 @@
#pragma once
#include <cmath>
#include <g3log/g3log.hpp>
#include "imstkMath.h"
namespace imstk
......@@ -37,26 +38,26 @@ public:
void setRadius(const Real radius)
{
m_radius = radius;
const auto pi = static_cast<Real>(PI);
const auto h2 = m_radius * m_radius;
const auto h3 = h2 * m_radius;
// if constexpr (N == 2) {
if(N == 2)
{
m_k = Real(40.0 / 7.0) / (pi * h2);
m_l = Real(240.0 / 7.0) / (pi * h2);
m_k = Real(40.0 / 7.0) / (PI * h2);
m_l = Real(240.0 / 7.0) / (PI * h2);
}
else
{
m_k = Real(8.0) / (pi * h3);
m_l = Real(48.0) / (pi * h3);
m_k = Real(8.0) / (PI * h3);
m_l = Real(48.0) / (PI * h3);
}
m_W_zero = W(VecXr::Zero());
}
Real W(Real r) const
Real W(const Real r) const
{
auto res = Real(0);
Real res = 0.;
const auto q = r / m_radius;
if(q <= Real(1.0))
{
......@@ -118,27 +119,27 @@ public:
{
m_radius = radius;
m_radius2 = m_radius * m_radius;
const auto pi = static_cast<Real>(PI);
// if constexpr (N == 2) {
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));
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_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_W_zero = W(VecXr::Zero());
}
/**
* W(r,h) = (315/(64 pi h^9))(h^2-|r|^2)^3
* = (315/(64 pi h^9))(h^2-r*r)^3
* W(r,h) = (315/(64 PI h^9))(h^2-|r|^2)^3
* = (315/(64 PI h^9))(h^2-r*r)^3
*/
Real W(Real r) const
Real W(const Real r) const
{
const auto r2 = r * r;
return (r2 <= m_radius2) ? std::pow(m_radius2 - r2, 3) * m_k : Real(0);
......@@ -151,8 +152,8 @@ public:
}
/**
* grad(W(r,h)) = r(-945/(32 pi h^9))(h^2-|r|^2)^2
* = r(-945/(32 pi h^9))(h^2-r*r)^2
* grad(W(r,h)) = r(-945/(32 PI h^9))(h^2-|r|^2)^2
* = r(-945/(32 PI h^9))(h^2-r*r)^2
*/
VecXr gradW(const VecXr& r) const
{
......@@ -168,12 +169,12 @@ public:
}
/**
* laplacian(W(r,h)) = (-945/(32 pi h^9))(h^2-|r|^2)(-7|r|^2+3h^2)
* = (-945/(32 pi h^9))(h^2-r*r)(3 h^2-7 r*r)
* laplacian(W(r,h)) = (-945/(32 PI h^9))(h^2-|r|^2)(-7|r|^2+3h^2)
* = (-945/(32 PI h^9))(h^2-r*r)(3 h^2-7 r*r)
*/
Real laplacianW(const VecXr& r) const
{
auto res = Real(0);
Real res = 0.;
const auto r2 = r.squaredNorm();
if(r2 <= m_radius2)
{
......@@ -198,7 +199,7 @@ protected:
template<int N>
class SpikyKernel {
class SPIkyKernel {
using VecXr = Eigen::Matrix<Real, N, 1>;
public:
......@@ -206,27 +207,27 @@ public:
{
m_radius = radius;
m_radius2 = m_radius * m_radius;
const auto pi = Real(PI);
// if constexpr (N == 2) {
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);
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_k = Real(15.0) / (PI * radius6);
m_l = -Real(45.0) / (PI * radius6);
}
m_W_zero = W(VecXr::Zero());
}
/**
* W(r,h) = 15/(pi*h^6) * (h-r)^3
* W(r,h) = 15/(PI*h^6) * (h-r)^3
*/
Real W(Real r) const { return (r <= m_radius) ? std::pow(m_radius - r, 3) * m_k : Real(0); }
Real W(const Real r) const { return (r <= m_radius) ? std::pow(m_radius - r, 3) * m_k : Real(0); }
Real W(const VecXr& r) const
{
......@@ -235,7 +236,7 @@ public:
}
/**
* grad(W(r,h)) = -r(45/(pi*h^6) * (h-r)^2)
* grad(W(r,h)) = -r(45/(PI*h^6) * (h-r)^2)
*/
VecXr gradW(const VecXr& r) const
{
......@@ -272,27 +273,27 @@ public:
{
m_radius = radius;
m_radius2 = m_radius * m_radius;
const auto pi = static_cast<Real>(PI);
// if constexpr (N == 2) {
if(N == 2)
{
throw std::runtime_error("Error: unimplemented!");
LOG(FATAL) << "Unimplemented function";
}
else
{
m_k = Real(32.0) / (pi * std::pow(m_radius, 9));
m_k = Real(32.0) / (PI * std::pow(m_radius, 9));
m_c = std::pow(m_radius, 6) / Real(64.0);
}
m_W_zero = W(VecXr::Zero());
}
/**
* 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
* 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(Real r) const
Real W(const Real r) const
{
auto res = Real(0);
Real res = 0.;
const auto r2 = r * r;
if(r2 <= m_radius2)
{
......@@ -312,7 +313,7 @@ public:
Real W(const VecXr& r) const
{
auto res = Real(0);
Real res = 0.;
const auto r2 = r.squaredNorm();
if(r2 <= m_radius2)
{
......@@ -349,10 +350,11 @@ public:
{
m_radius = radius;
m_radius2 = m_radius * m_radius;
// if constexpr (N == 2) {
if(N == 2)
{
throw std::runtime_error("Error: unimplemented!");
LOG(FATAL) << "Unimplemented function";
}
else
{
......@@ -364,9 +366,9 @@ public:
/**
* W(r,h) = (0.007/h^3.25)(-4r^2/h + 6r -2h)^0.25 if h/2 < r <= h
*/
Real W(Real r) const
Real W(const Real r) const
{
auto res = Real(0);
Real res = 0.;
const auto r2 = r * r;
if(r2 <= m_radius2)
{
......@@ -381,7 +383,7 @@ public:
Real W(const VecXr& r) const
{
auto res = Real(0);
Real res = 0.;
const auto r2 = r.squaredNorm();
if(r2 <= m_radius2)
{
......@@ -417,7 +419,7 @@ public:
Real laplace(const VecXr& r) const
{
auto res = Real(0);
Real res = 0.;
const auto r2 = r.squaredNorm();
if(r2 <= m_radius2)
{
......@@ -444,20 +446,20 @@ public:
void initialize(const Real kernelRadius)
{
m_poly6.setRadius(kernelRadius);
m_spiky.setRadius(kernelRadius);
m_sPIky.setRadius(kernelRadius);
m_viscosity.setRadius(kernelRadius);
m_cohesion.setRadius(kernelRadius);
}
auto W_zero() const { return m_poly6.W_zero(); }
auto W(const Vec3r& r) const { return m_poly6.W(r); }
auto gradW(const Vec3r& r) const { return m_spiky.gradW(r); }
auto gradW(const Vec3r& r) const { return m_sPIky.gradW(r); }
auto laplace(const Vec3r& r) const { return m_viscosity.laplace(r); }
auto cohesionW(const Vec3r& r) const { return m_cohesion.W(r); }
protected:
SPH::Poly6Kernel<3> m_poly6;
SPH::SpikyKernel<3> m_spiky;
SPH::SPIkyKernel<3> m_sPIky;
SPH::ViscosityKernel<3> m_viscosity;
SPH::CohesionKernel<3> m_cohesion;
};
......
......@@ -35,10 +35,8 @@ SPHModelConfig::SPHModelConfig(const Real particleRadius)
void SPHModelConfig::initialize()
{
if(std::abs(m_ParticleRadius) < Real(1e-6) )
{
LOG(FATAL) << "Particle radius was not set properly";
}
LOG_IF(FATAL, (std::abs(m_ParticleRadius) < Real(1e-6))) << "Particle radius was not set properly";
m_ParticleRadiusSqr = m_ParticleRadius * m_ParticleRadius;
m_ParticleMass = Real(std::pow(Real(2.0) * m_ParticleRadius, 3)) * m_RestDensity * m_ParticleMassScale;
......@@ -53,10 +51,7 @@ void SPHModelConfig::initialize()
bool SPHModel::initialize()
{
if(!m_Geometry)
{
LOG(FATAL) << "Model geometry is not yet set! Cannot initialize without model geometry.";
}
LOG_IF(FATAL, !m_Geometry) << "Model geometry is not yet set! Cannot initialize without model geometry.";
// initialize positions and velocity of the particles
this->m_initialState = std::make_shared<SPHKinematicState>();
......@@ -163,7 +158,7 @@ void SPHModel::findParticleNeighbors()
searcher.insertPoints(getState().getPositions());
runLoop(getState().size(),
[&] (size_t p) {
[&] (const size_t p) {
const auto& ppos = getState().getPositions()[p];
searcher.getPointsInSphere(getState().getFluidNeighborLists()[p], ppos, m_Parameters->m_KernelRadius);
});
......@@ -174,7 +169,7 @@ void SPHModel::findParticleNeighbors()
searcher.insertPoints(getState().getBoundaryParticlePositions());
runLoop(getState().size(),
[&] (size_t p) {
[&] (const size_t p) {
const auto& ppos = getState().getPositions()[p];
searcher.getPointsInSphere(getState().getBoundaryNeighborLists()[p], ppos, m_Parameters->m_KernelRadius);
});
......@@ -195,7 +190,7 @@ void SPHModel::computeNeighborRelativePositions()
};
////////////////////////////////////////////////////////////////////////////////
runLoop(getState().size(),
[&](size_t p) {
[&] (const size_t p) {
const auto& ppos = getState().getPositions()[p];
auto& neighborInfo = getState().getNeighborInfo()[p];
neighborInfo.resize(0);
......@@ -217,7 +212,7 @@ void SPHModel::collectNeighborDensity()
// this is usefull because relative positions and densities are accessed together multiple times
// caching relative positions and densities therefore can reduce computation time significantly (tested)
runLoop(getState().size(),
[&](size_t p) {
[&] (const size_t p) {
auto& neighborInfo = getState().getNeighborInfo()[p];
assert(neighborInfo.size() >= 1); // neighbor of particle p should contain at least p index
if(neighborInfo.size() <= 1)
......@@ -238,7 +233,7 @@ void SPHModel::collectNeighborDensity()
void SPHModel::computeDensity()
{
runLoop(getState().size(),
[&](size_t p) {
[&] (const size_t p) {
const auto& neighborInfo = getState().getNeighborInfo()[p];
assert(neighborInfo.size() >= 1); // neighbor of particle p should contain at least p index
if(neighborInfo.size() <= 1)
......@@ -266,7 +261,7 @@ void SPHModel::normalizeDensity()
getState().getNormalizedDensities().resize(getState().size());
////////////////////////////////////////////////////////////////////////////////
runLoop(getState().size(),
[&](size_t p) {
[&] (const size_t p) {
auto& neighborInfo = getState().getNeighborInfo()[p];
assert(neighborInfo.size() >= 1); // neighbor of particle p should contain at least p index
if(neighborInfo.size() <= 1)
......@@ -316,7 +311,7 @@ void SPHModel::computePressureAcceleration()
};
runLoop(getState().size(),
[&](size_t p) {
[&] (const size_t p) {
Vec3r accel(0, 0, 0);
const auto& neighborInfo = getState().getNeighborInfo()[p];
assert(neighborInfo.size() >= 1); // neighbor of particle p should contain at least p index
......@@ -345,10 +340,10 @@ void SPHModel::computePressureAcceleration()
}
void SPHModel::updateVelocity(Real timestep)
void SPHModel::updateVelocity(const Real timestep)
{
runLoop(getState().size(),
[&](size_t p) {
[&] (const size_t p) {
getState().getVelocities()[p] += (m_Parameters->m_Gravity + getState().getAccelerations()[p]) * timestep;
});
}
......@@ -357,7 +352,7 @@ void SPHModel::updateVelocity(Real timestep)
void SPHModel::computeViscosity()
{
runLoop(getState().size(),
[&](size_t p) {
[&] (const size_t p) {
const auto& neighborInfo = getState().getNeighborInfo()[p];
assert(neighborInfo.size() >= 1); // neighbor of particle p should contain at least p index
if(neighborInfo.size() <= 1)
......@@ -397,7 +392,7 @@ void SPHModel::computeViscosity()
// add diffused velocity back to velocity, causing viscosity
runLoop(getState().size(),
[&](size_t p) {
[&] (const size_t p) {
getState().getVelocities()[p] += getState().getDiffuseVelocities()[p];
});
}
......@@ -406,7 +401,7 @@ void SPHModel::computeViscosity()
void SPHModel::computeNormal()
{
runLoop(getState().size(),
[&](size_t p) {
[&] (const size_t p) {
Vec3r n(0, 0, 0);
const auto& neighborInfo = getState().getNeighborInfo()[p];
assert(neighborInfo.size() >= 1); // neighbor of particle p should contain at least p index
......@@ -437,7 +432,7 @@ void SPHModel::computeSurfaceTension()
// Compute surface tension for each particle, using Akinci et at. 2013 model
runLoop(getState().size(),
[&](size_t p) {
[&] (const size_t p) {
const auto& fluidNeighborList = getState().getFluidNeighborLists()[p];
assert(fluidNeighborList.size() >= 1); // Neighbor of particle p should contain at least p index
if(fluidNeighborList.size() <= 1)
......@@ -482,10 +477,10 @@ void SPHModel::computeSurfaceTension()
}
void SPHModel::advect(Real timestep)
void SPHModel::advect(const Real timestep)
{
runLoop(getState().size(),
[&](size_t p) {
[&] (const size_t p) {
getState().getPositions()[p] += getState().getVelocities()[p] * timestep;
});
}
......