Commit 66e54157 authored by Nghia Truong's avatar Nghia Truong
Browse files

REFAC: Rewrite all template classes to use Real as type aliasing

parent 1b108566
......@@ -105,7 +105,7 @@ StdVectorOfVec3d generateBunnyShapeFluid(double particleRadius)
}
std::shared_ptr<SPHObjectD> generateFluid(const std::shared_ptr<Scene>&scene, int sceneIdx, double particleRadius)
std::shared_ptr<SPHObject> generateFluid(const std::shared_ptr<Scene>&scene, int sceneIdx, double particleRadius)
{
StdVectorOfVec3d particles;
switch(sceneIdx)
......@@ -130,7 +130,7 @@ std::shared_ptr<SPHObjectD> generateFluid(const std::shared_ptr<Scene>&scene, in
fluidGeometry->initialize(particles);
// Create a fluids object
auto fluidObj = std::make_shared<SPHObjectD>("Sphere");
auto fluidObj = std::make_shared<SPHObject>("Sphere");
// Create a visiual model
auto fluidVisualModel = std::make_shared<VisualModel>(fluidGeometry);
......@@ -140,11 +140,11 @@ std::shared_ptr<SPHObjectD> generateFluid(const std::shared_ptr<Scene>&scene, in
fluidVisualModel->setRenderMaterial(fluidMaterial);
// Create a physics model
auto sphModel = std::make_shared<SPHModelD>();
auto sphModel = std::make_shared<SPHModel>();
sphModel->setModelGeometry(fluidGeometry);
// configure model
auto sphParams = std::make_shared<SPHParametersD>(particleRadius);
auto sphParams = std::make_shared<SPHModelConfig>(particleRadius);
sphParams->m_bNormalizeDensity = true;
if(sceneIdx == 2) // highly viscous fluid
{
......@@ -164,7 +164,7 @@ std::shared_ptr<SPHObjectD> generateFluid(const std::shared_ptr<Scene>&scene, in
scene->addSceneObject(fluidObj);
// Configure the solver
auto sphSolver = std::make_shared<SPHSolverD>();
auto sphSolver = std::make_shared<SPHSolver>();
sphSolver->setSPHObject(fluidObj);
scene->addNonlinearSolver(sphSolver);
......
......@@ -26,7 +26,7 @@
using namespace imstk;
std::shared_ptr<SPHObjectD> generateFluid(const std::shared_ptr<Scene>&scene, int sceneIdx, double particleRadius);
std::shared_ptr<SPHObject> 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);
///
......
......@@ -339,7 +339,7 @@ void neighborSearch_GridBased(Data& verts, std::vector<std::vector<size_t>>& nei
}
const double radius = 4.000000000000001 * PARTICLE_RADIUS;
static GridBasedNeighborSearch<double> gridSearch;
static GridBasedNeighborSearch gridSearch;
gridSearch.setSearchRadius(radius);
gridSearch.getNeighbors(neighbors, verts);
......
......@@ -83,8 +83,7 @@ CollisionHandling::make_collision_handling(const Type& type,
return std::make_shared<BoneDrillingCH>(side, colData, objA, objB);
case Type::SPH:
return std::make_shared<SPHCollisionHandling<double>>(side, colData, objA);
// TODO: instead of using <double>, use SPHCollisionHandling<Real>
return std::make_shared<SPHCollisionHandling>(side, colData, objA);
default:
LOG(WARNING) << "CollisionHandling::make_collision_handling error: type not implemented.";
......
......@@ -23,8 +23,7 @@
namespace imstk
{
template<class Real>
void SPHCollisionHandling<Real>::setBoundaryFriction(Real frictionLength)
void SPHCollisionHandling::setBoundaryFriction(Real frictionLength)
{
if(frictionLength < 0.0 || frictionLength > 1.0)
{
......@@ -41,14 +40,11 @@ void SPHCollisionHandling<Real>::setBoundaryFriction(Real frictionLength)
m_BoundaryFriction = frictionLength;
}
////////////////////////////////////////////////////////////////////////////////
template<class Real>
void SPHCollisionHandling<Real>::computeContactForces()
void SPHCollisionHandling::computeContactForces()
{
const auto& SPHModel = m_SPHObject->getSPHModel();
assert(SPHModel);
auto& state = SPHModel->getState();
using Vec3r = Eigen::Matrix<Real, 3, 1>;
for(const auto& cd : m_colData.MAColData)
{
......@@ -84,5 +80,3 @@ void SPHCollisionHandling<Real>::computeContactForces()
}
}
} // end namespace imstk
template class imstk::SPHCollisionHandling<double>;
......@@ -25,12 +25,11 @@
namespace imstk
{
template<class Real>
class SPHCollisionHandling : public CollisionHandling
{
public:
SPHCollisionHandling(const Side& side, const CollisionData& colData, const std::shared_ptr<CollidingObject>& obj) :
CollisionHandling(Type::SPH, side, colData), m_SPHObject(std::dynamic_pointer_cast<SPHObject<Real>>(obj)) {}
CollisionHandling(Type::SPH, side, colData), m_SPHObject(std::dynamic_pointer_cast<SPHObject>(obj)) {}
SPHCollisionHandling() = delete;
virtual ~SPHCollisionHandling() = default;
......@@ -51,7 +50,7 @@ public:
Real getBoundaryFriction() const { return m_BoundaryFriction; }
private:
std::shared_ptr<SPHObject<Real>> m_SPHObject;
std::shared_ptr<SPHObject> m_SPHObject;
Real m_BoundaryFriction = 0.1;
};
} // end namespace imstk
\ No newline at end of file
......@@ -30,6 +30,16 @@
namespace imstk
{
// Define Real type and dependent types
using Real = double;
using Vec2r = Eigen::Matrix<Real, 2, 1>;
using Vec3r = Eigen::Matrix<Real, 3, 1>;
using Vec4r = Eigen::Matrix<Real, 4, 1>;
using StdVectorOfReal = std::vector<Real>;
using StdVectorOfVec2r = std::vector<Vec2r, Eigen::aligned_allocator<Vec2r>>;
using StdVectorOfVec3r = std::vector<Vec3r, Eigen::aligned_allocator<Vec3r>>;
using StdVectorOfVec4r = std::vector<Vec4r, Eigen::aligned_allocator<Vec4r>>;
// 2D vector
using Vec2f = Eigen::Vector2f;
using Vec2d = Eigen::Vector2d;
......@@ -100,25 +110,29 @@ using AffineTransform3d = Eigen::Affine3d;
#define WORLD_ORIGIN Vec3d::Zero()
/// Some commonly used math constants
const double PI = 3.14159265358979323846;
const double PI_2 = 1.57079632679489661923;
const double PI_4 = 0.785398163397448309616;
const double INV_1_PI = 0.318309886183790671538;
const double INV_2_PI = 0.636619772367581343076;
const double TWO_OVER_SQRTPI = 1.12837916709551257390;
const double SQRT2 = 1.41421356237309504880;
const double SQRT1_2 = 0.707106781186547524401;
const double NLOG_E = 2.71828182845904523536;
const double LOG2E = 1.44269504088896340736;
const double LOG10E = 0.434294481903251827651;
const double LN2 = 0.693147180559945309417;
const double LN10 = 2.30258509299404568402;
const double MAX_D = std::numeric_limits<double>::max();
const double MIN_D = std::numeric_limits<double>::min();
const double MAX_F = std::numeric_limits<float>::max();
const double MIN_F = std::numeric_limits<float>::min();
const double MACHINE_PRECISION = std::numeric_limits<double>::denorm_min();
static constexpr Real PI = 3.14159265358979323846;
static constexpr Real PI_2 = 1.57079632679489661923;
static constexpr Real PI_4 = 0.785398163397448309616;
static constexpr Real INV_1_PI = 0.318309886183790671538;
static constexpr Real INV_2_PI = 0.636619772367581343076;
static constexpr Real TWO_OVER_SQRTPI = 1.12837916709551257390;
static constexpr Real SQRT2 = 1.41421356237309504880;
static constexpr Real SQRT1_2 = 0.707106781186547524401;
static constexpr Real NLOG_E = 2.71828182845904523536;
static constexpr Real LOG2E = 1.44269504088896340736;
static constexpr Real LOG10E = 0.434294481903251827651;
static constexpr Real LN2 = 0.693147180559945309417;
static constexpr Real LN10 = 2.30258509299404568402;
static constexpr double MAX_REAL = std::numeric_limits<Real>::max();
static constexpr double MIN_REAL = std::numeric_limits<Real>::min();
static constexpr double MAX_D = std::numeric_limits<double>::max();
static constexpr double MIN_D = std::numeric_limits<double>::min();
static constexpr float MAX_F = std::numeric_limits<float>::max();
static constexpr float MIN_F = std::numeric_limits<float>::min();
static constexpr Real MACHINE_PRECISION = std::numeric_limits<Real>::denorm_min();
static constexpr double MACHINE_PRECISION_D = std::numeric_limits<double>::denorm_min();
static constexpr double MACHINE_PRECISION_F = std::numeric_limits<double>::denorm_min();
}
#endif // ifndef imstkMath_h
......@@ -28,9 +28,6 @@
using namespace imstk;
using Real = double;
using Vec3r = Vec3d;
#define SPHERE_RADIUS Real(1)
#define SPHERE_CENTER Vec3r(0, 0, 0)
#define PARTICLE_RADIUS Real(0.05)
......@@ -91,7 +88,7 @@ void neighborSearch_GridBased(Data& verts, std::vector<std::vector<size_t>>& nei
{
neighbors.resize(verts.size());
const Real radius = Real(4.000000000000001) * PARTICLE_RADIUS;
static GridBasedNeighborSearch<Real> gridSearch;
static GridBasedNeighborSearch gridSearch;
gridSearch.setSearchRadius(radius);
gridSearch.getNeighbors(neighbors, verts);
}
......@@ -157,7 +154,7 @@ void neighborSearch_GridBased(Data& setA, Data& setB, std::vector<std::vector<si
{
neighbors.resize(setA.size());
const Real radius = Real(4.000000000000001) * PARTICLE_RADIUS;
static GridBasedNeighborSearch<Real> gridSearch;
static GridBasedNeighborSearch gridSearch;
gridSearch.setSearchRadius(radius);
gridSearch.getNeighbors(neighbors, setA, setB);
}
......
......@@ -24,30 +24,26 @@
namespace imstk
{
template<class Real>
void GridBasedNeighborSearch<Real>::setSearchRadius(const Real radius)
void GridBasedNeighborSearch::setSearchRadius(const Real radius)
{
m_SearchRadius = radius;
m_SearchRadiusSqr = radius * radius;
}
template<class Real>
std::vector<std::vector<size_t>>
GridBasedNeighborSearch<Real>::getNeighbors(const StdVT_Vec3r& points)
GridBasedNeighborSearch::getNeighbors(const StdVectorOfVec3r& points)
{
std::vector<std::vector<size_t>> result;
getNeighbors(result, points, points);
return result;
}
template<class Real>
void GridBasedNeighborSearch<Real>::getNeighbors(std::vector<std::vector<size_t>>&result, const StdVT_Vec3r& points)
void GridBasedNeighborSearch::getNeighbors(std::vector<std::vector<size_t>>&result, const StdVectorOfVec3r& points)
{
getNeighbors(result, points, points);
}
template<class Real>
void GridBasedNeighborSearch<Real>::getNeighbors(std::vector<std::vector<size_t>>&result, const StdVT_Vec3r& setA, const StdVT_Vec3r& setB)
void GridBasedNeighborSearch::getNeighbors(std::vector<std::vector<size_t>>&result, const StdVectorOfVec3r& setA, const StdVectorOfVec3r& setB)
{
if(std::abs(m_SearchRadius) < Real(1e-8))
{
......@@ -55,9 +51,8 @@ void GridBasedNeighborSearch<Real>::getNeighbors(std::vector<std::vector<size_t>
}
// firstly compute the bounding box of points in setB
const auto hugeNumber = std::numeric_limits<Real>::max();
Vec3r lowerCorner(hugeNumber, hugeNumber, hugeNumber);
Vec3r upperCorner(-hugeNumber, -hugeNumber, -hugeNumber);
Vec3r lowerCorner(MAX_REAL, MAX_REAL, MAX_REAL);
Vec3r upperCorner(-MAX_REAL, -MAX_REAL, -MAX_REAL);
for(auto& ppos : setB)
{
......@@ -133,6 +128,3 @@ void GridBasedNeighborSearch<Real>::getNeighbors(std::vector<std::vector<size_t>
}
}
} // end namespace imstk
// Explicit instantiate class GridBasedNeighborSearch
template class imstk::GridBasedNeighborSearch<double>;
\ No newline at end of file
......@@ -28,12 +28,8 @@ namespace imstk
///
/// \brief Class for searching neighbors using regular grid
///
template<class Real>
class GridBasedNeighborSearch
{
using Vec3r = Eigen::Matrix<Real, 3, 1>;
using StdVT_Vec3r = std::vector<Vec3r, Eigen::aligned_allocator<Vec3r>>;
public:
GridBasedNeighborSearch() = default;
......@@ -54,14 +50,14 @@ public:
/// \param points The given points to search for neighbors
/// \return List of list of neighbor indices for each point
///
std::vector<std::vector<size_t>> getNeighbors(const StdVT_Vec3r& points);
std::vector<std::vector<size_t>> getNeighbors(const StdVectorOfVec3r& points);
///
/// \brief Search neighbors for each point within the search radius
/// \param result The list of lists of neighbor indices for each point
/// \param points The given points to search for neighbors
///
void getNeighbors(std::vector<std::vector<size_t>>& result, const StdVT_Vec3r& points);
void getNeighbors(std::vector<std::vector<size_t>>& result, const StdVectorOfVec3r& points);
///
/// \brief Search neighbors from setB for each point in setA within the search radius. SetA and setB can be different.
......@@ -69,11 +65,11 @@ public:
/// \param setA The point set for which performing neighbor search
/// \param setB The point set where neighbor indices will be collected
///
void getNeighbors(std::vector<std::vector<size_t>>& result, const StdVT_Vec3r& setA, const StdVT_Vec3r& setB);
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);
UniformSpatialGrid<Real, std::vector<size_t>> m_Grid;
UniformSpatialGrid<std::vector<size_t>> m_Grid;
};
} // end namespace imstk
\ No newline at end of file
} // end namespace imstk
/*=========================================================================
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 "imstkUniformSpatialGrid.h"
namespace imstk
{
template<class Real, class CellData>
void UniformSpatialGrid<Real, CellData>::initialize(const Vec3r& lowerCorner, const Vec3r& upperCorner, const Real cellSize)
{
assert(cellSize > 0);
m_LowerCorner = lowerCorner;
m_UpperCorner = upperCorner;
m_CellSize = cellSize;
m_InvCellSize = Real(1.0) / m_CellSize;
m_NTotalCells = 1u;
for(int i = 0; i < 3; ++i)
{
m_Resolution[i] = static_cast<unsigned int>(std::ceil((m_UpperCorner[i] - m_LowerCorner[i]) / m_CellSize));
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);
}
// cell data must be resized to equal to number of cells
m_CellData.resize(m_NTotalCells);
}
} // end namespace imstk
// Explicit instantiate class UniformSpatialGrid
template class imstk::UniformSpatialGrid<double, std::vector<size_t>>;
\ No newline at end of file
......@@ -34,11 +34,9 @@ namespace imstk
/// \param Real The floating point type for computation (must be float or double
/// \param CellData The data type that each cell will hold an instance of it
///
template<class Real, class CellData>
template<class CellData>
class UniformSpatialGrid
{
using Vec3r = Eigen::Matrix<Real, 3, 1>;
public:
UniformSpatialGrid() = default;
UniformSpatialGrid(const Vec3r& lowerCorner, const Vec3r& upperCorner, Real cellSize) { initialize(lowerCorner, upperCorner, cellSize); }
......@@ -49,11 +47,37 @@ public:
/// \param upperCorner The upper corner of the grid
/// \param cellSize the edge length of grid cell
///
void initialize(const Vec3r& lowerCorner, const Vec3r& upperCorner, const Real cellSize);
void initialize(const Vec3r& lowerCorner, const Vec3r& upperCorner, const Real cellSize)
{
assert(cellSize > 0);
m_LowerCorner = lowerCorner;
m_UpperCorner = upperCorner;
m_CellSize = cellSize;
m_InvCellSize = Real(1.0) / m_CellSize;
m_NTotalCells = 1u;
for(int i = 0; i < 3; ++i)
{
m_Resolution[i] = static_cast<unsigned int>(std::ceil((m_UpperCorner[i] - m_LowerCorner[i]) / m_CellSize));
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);
}
// cell data must be resized to equal to number of cells
m_CellData.resize(m_NTotalCells);
}
///
/// \brief Get number of grid cell in 3 dimensions: (num_cell_x, num_cell_y, num_cell_z)
//
///
auto getResolution() const { return m_Resolution; }
///
......@@ -172,4 +196,4 @@ private:
std::vector<CellData> m_CellData;
};
} // end namespace imstk
\ No newline at end of file
} // end namespace imstk
......@@ -29,7 +29,7 @@ namespace imstk
{
namespace SPH
{
template<int N, class Real>
template<int N>
class CubicKernel {
using VecXr = Eigen::Matrix<Real, N, 1>;
......@@ -109,7 +109,7 @@ protected:
};
template<int N, class Real>
template<int N>
class Poly6Kernel {
using VecXr = Eigen::Matrix<Real, N, 1>;
......@@ -197,7 +197,7 @@ protected:
};
template<int N, class Real>
template<int N>
class SpikyKernel {
using VecXr = Eigen::Matrix<Real, N, 1>;
......@@ -263,7 +263,7 @@ protected:
};
template<int N, class Real>
template<int N>
class CohesionKernel {
using VecXr = Eigen::Matrix<Real, N, 1>;
......@@ -340,8 +340,7 @@ protected:
Real m_W_zero;
};
template<int N, class Real>
template<int N>
class AdhesionKernel {
using VecXr = Eigen::Matrix<Real, N, 1>;
......@@ -404,7 +403,7 @@ protected:
Real m_W_zero;
};
template<int N, class Real>
template<int N>
class ViscosityKernel {
using VecXr = Eigen::Matrix<Real, N, 1>;
......@@ -438,11 +437,9 @@ protected:
///
/// \brief Struct contains SPH kernels for time integration, using different kernel for different purposes
///
template<class Real>
class SPHSimulationKernels
{
using Vec3r = Eigen::Matrix<Real, 3, 1>;
public:
void initialize(const Real kernelRadius)
{
......@@ -459,9 +456,9 @@ public:
auto cohesionW(const Vec3r& r) const { return m_cohesion.W(r); }
protected:
SPH::Poly6Kernel<3, Real> m_poly6;
SPH::SpikyKernel<3, Real> m_spiky;
SPH::ViscosityKernel<3, Real> m_viscosity;
SPH::CohesionKernel<3, Real> m_cohesion;
SPH::Poly6Kernel<3> m_poly6;
SPH::SpikyKernel<3> m_spiky;
SPH::ViscosityKernel<3> m_viscosity;
SPH::CohesionKernel<3> m_cohesion;
};
} // end namespace imstk
......@@ -25,15 +25,15 @@
namespace imstk
{
// SPHModelConfig implementation ===>
template<class Real>
SPHModelConfig<Real>::SPHModelConfig(const Real particleRadius)
SPHModelConfig::SPHModelConfig(const Real particleRadius)
{
m_ParticleRadius = particleRadius;
initialize();
}
template<class Real>
void SPHModelConfig<Real>::initialize()
void SPHModelConfig::initialize()
{
if(std::abs(m_ParticleRadius) < Real(1e-6) )
{
......@@ -50,8 +50,8 @@ void SPHModelConfig<Real>::initialize()
}
// SPHModel implementation ===>
template<class Real>
bool SPHModel<Real>::initialize()
bool SPHModel::initialize()
{
if(!m_Geometry)
{
......@@ -59,8 +59,8 @@ bool SPHModel<Real>::initialize()
}
// initialize positions and velocity of the particles
this->m_initialState = std::make_shared<SPHKinematicState<Real>>();
this->m_currentState = std::make_shared<SPHKinematicState<Real>>();
this->m_initialState = std::make_shared<SPHKinematicState>();
this->m_currentState = std::make_shared<SPHKinematicState>();
// Set particle positions and zero default velocities
// TODO: set particle data with given (non-zero) velocities
......@@ -91,8 +91,8 @@ bool SPHModel<Real>::initialize()
return true;
}
template<class Real>
void SPHModel<Real>::simulationTimeStep()
void SPHModel::simulationTimeStep()
{
findParticleNeighbors();
computeNeighborRelativePositions();
......@@ -107,14 +107,14 @@ void SPHModel<Real>::simulationTimeStep()
advect(getTimeStep());
}
template<class Real>
void SPHModel<Real>::computeTimeStepSize()
void SPHModel::computeTimeStepSize()
{
m_dt = (this->m_timeStepSizeType == TimeSteppingType::fixed) ? m_DefaultDt : computeCFLTimeStepSize();
}
template<class Real>
Real SPHModel<Real>::computeCFLTimeStepSize()
Real SPHModel::computeCFLTimeStepSize()
{
Real maxVel = 0;
for(const auto& vel : getState().getVelocities())
......@@ -143,8 +143,8 @@ Real SPHModel<Real>::computeCFLTimeStepSize()
return timestep;
}
template<class Real>
void SPHModel<Real>::findParticleNeighbors()
void SPHModel::findParticleNeighbors()
{
if(m_Parameters->m_NeighborSearchMethod == NeighborSearchMethod::GridBased)
{
......@@ -182,8 +182,8 @@ void SPHModel<Real>::findParticleNeighbors()
}