Commit cf0dbc3a authored by Sreekanth Arikatla's avatar Sreekanth Arikatla

Merge branch 'minorCleanUp' into 'master'

Minor clean up

See merge request iMSTK/iMSTK!350
parents 1ec0fe1d eecca310
......@@ -128,10 +128,6 @@ main()
scene->addLight(colorLight);
scene->addSceneObject(deformableObj);
// print UPS
/*auto ups = std::make_shared<UPSCounter>();
apiutils::printUPS(sdk->getSceneManager(scene), ups);*/
scene->getCamera()->setFocalPoint(0, -5, 5);
scene->getCamera()->setPosition(-15., -5.0, 15.0);
......
......@@ -136,7 +136,7 @@ generateFluid(const std::shared_ptr<Scene>& scene, const double particleRadius)
// Create a fluids object
auto fluidObj = std::make_shared<SPHObject>("Sphere");
// Create a visiual model
// Create a visual model
auto fluidVisualModel = std::make_shared<VisualModel>(fluidGeometry);
auto fluidMaterial = std::make_shared<RenderMaterial>();
fluidMaterial->setColor(Color::Blue);
......@@ -152,14 +152,14 @@ generateFluid(const std::shared_ptr<Scene>& scene, const double particleRadius)
sphParams->m_bNormalizeDensity = true;
if (SCENE_ID == 2) // highly viscous fluid
{
sphParams->m_RatioKernelOverParticleRadius = 6.0;
sphParams->m_ViscosityFluid = 0.5;
sphParams->m_SurfaceTensionStiffness = 5.0;
sphParams->m_kernelOverParticleRadiusRatio = 6.0;
sphParams->m_viscosityCoeff = 0.5;
sphParams->m_surfaceTensionStiffness = 5.0;
}
if (SCENE_ID == 3) // bunny-shaped fluid
{
sphParams->m_BoundaryFriction = 0.3;
sphParams->m_frictionBoundary= 0.3;
}
sphModel->configure(sphParams);
......
......@@ -95,7 +95,7 @@ BoneDrillingCH::erodeBone()
if (m_nodalDensity[cd.nodeIdx] <= 0.)
{
// TODO: Unused variable, maybe used in furture?
/// \todo Unused variable, maybe used in furture?
// lock.lock();
// m_erodedNodes.push_back(cd.nodeId);
// lock.unlock();
......
......@@ -89,7 +89,7 @@ private:
std::vector<double> m_nodalDensity; ///> Density of the bone
double m_initialBoneDensity = 1.0; ///> Density of the bone before the start of the drilling process
// std::vector<size_t> m_erodedNodes; // TODO: Unused variable
// std::vector<size_t> m_erodedNodes; /// \todo Unused variable
std::vector<bool> m_nodeRemovalStatus; ///> Keeps track of the removal status of the node
std::vector<std::vector<size_t>> m_nodalCardinalSet; ///> Keeps track of the removal status of the node
......
......@@ -144,6 +144,6 @@ PBDCollisionHandling::generatePBDConstraints()
m_PBDConstraints.push_back(c);
}
//TODO: generating PbdPointTriangleConstraint from the VTColData should be added
//\todo generating PbdPointTriangleConstraint from the VTColData should be added
}
}
......@@ -42,7 +42,7 @@ SPHCollisionHandling::processCollisionData()
LOG_IF(FATAL, (!SPHModel)) << "SPH model was not initialized";
#endif
const auto boundaryFriction = m_SPHObject->getSPHModel()->getParameters()->m_BoundaryFriction;
const auto boundaryFriction = m_SPHObject->getSPHModel()->getParameters()->m_frictionBoundary;
#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])";
......
......@@ -33,7 +33,7 @@ PbdConstantDensityConstraint::initConstraint(PbdModel& model, const double)
// constraint parameters
// Refer: Miller, et al 2003 "Particle-Based Fluid Simulation for Interactive Applications."
// TODO: Check if these numbers can be variable
/// \todo Check if these numbers can be variable
m_wPoly6Coeff = 315.0 / (64.0 * PI * pow(m_maxDist, 9));
m_wSpikyCoeff = 15.0 / (PI * pow(m_maxDist, 6));
......
......@@ -21,7 +21,7 @@
#pragma once
// TODO: Remove this in TBB 2019 Update 4: https://github.com/intel/tbb/blob/tbb_2019/CHANGES#L117
/// \todo Remove this in TBB 2019 Update 4: https://github.com/intel/tbb/blob/tbb_2019/CHANGES#L117
#define TBB_PREVIEW_GLOBAL_CONTROL 1
#include <tbb/tbb.h>
......
......@@ -23,6 +23,8 @@
#include <iostream>
/// \todo remove nameless union/struct in the future
#pragma warning(disable : 4201)
namespace imstk
{
///
......@@ -102,4 +104,5 @@ struct Color
static Color Orange;
static Color Yellow;
};
#pragma warning(default : 4201)
}
......@@ -30,7 +30,7 @@ namespace imstk
///
/// \struct stdSink
///
/// \brief
/// \brief A standard sink that prints the message to a standard output
///
struct stdSink
{
......
......@@ -101,12 +101,14 @@ public:
/// \brief Check if cell index in dimension d is valid (d = 0/1/2 => x/y/z dimension)
///
template<int d>
bool isValidCellIndex(const int idx) const { return idx >= 0 && static_cast<unsigned int>(idx) < m_Resolution[d]; }
bool isValidCellIndex(const int idx) const
{ return idx >= 0 && static_cast<unsigned int>(idx) < m_Resolution[d]; }
///
/// \brief Check if 3D cell indices are valid
///
bool isValidCellIndices(const int i, const int j, const int k) const { return isValidCellIndex<0>(i) && isValidCellIndex<1>(j) && isValidCellIndex<2>(k); }
bool isValidCellIndices(const int i, const int j, const int k) const
{ return isValidCellIndex<0>(i) && isValidCellIndex<1>(j) && isValidCellIndex<2>(k); }
///
/// \brief Get the 3D index (cell_x, cell_y, cell_z) of the cell containing the given positions
......@@ -136,53 +138,77 @@ public:
/// \brief Get data in a cell
/// \param A position in space
///
CellData& getCellData(const Vec3r& ppos) { return m_CellData[getCellLinearizedIndex < unsigned int > (ppos)]; }
CellData& getCellData(const Vec3r& ppos)
{
return m_CellData[getCellLinearizedIndex < unsigned int > (ppos)];
}
///
/// \brief Get data in a cell
/// \param A position in space
///
const CellData& getCellData(const Vec3r& ppos) const { return m_CellData[getCellLinearizedIndex < unsigned int > (ppos)]; }
const CellData& getCellData(const Vec3r& ppos) const
{ return m_CellData[getCellLinearizedIndex < unsigned int > (ppos)]; }
///
/// \brief Get data in a cell
/// \param A linearized index of cell
///
CellData& getCellData(size_t linearizedIdx) { assert(linearizedIdx < m_CellData.size()); return m_CellData[linearizedIdx]; }
CellData& getCellData(size_t linearizedIdx)
{
assert(linearizedIdx < m_CellData.size());
return m_CellData[linearizedIdx];
}
///
/// \brief Get data in a cell
/// \param A linearized index of cell
///
const CellData& getCellData(size_t linearizedIdx) const { assert(linearizedIdx < m_CellData.size()); return m_CellData[linearizedIdx]; }
const CellData& getCellData(size_t linearizedIdx) const
{
assert(linearizedIdx < m_CellData.size());
return m_CellData[linearizedIdx];
}
///
/// \brief Get data in a cell
/// \param 3D index of a cell
///
template<class IndexType>
CellData& getCellData(const std::array<IndexType, 3>& cellIdx) { return m_CellData[getCellLinearizedIndex(cellIdx[0], cellIdx[1], cellIdx[2])]; }
CellData& getCellData(const std::array<IndexType, 3>& cellIdx)
{
return m_CellData[getCellLinearizedIndex(cellIdx[0], cellIdx[1], cellIdx[2])];
}
///
/// \brief Get data in a cell
/// \param 3D index of a cell
///
template<class IndexType>
const CellData& getCellData(const std::array<IndexType, 3>& cellIdx) const { return m_CellData[getCellLinearizedIndex(cellIdx[0], cellIdx[1], cellIdx[2])]; }
const CellData& getCellData(const std::array<IndexType, 3>& cellIdx) const
{
return m_CellData[getCellLinearizedIndex(cellIdx[0], cellIdx[1], cellIdx[2])];
}
///
/// \brief Get data in a cell
/// \param 3D index of a cell
///
template<class IndexType>
CellData& getCellData(const IndexType i, const IndexType j, const IndexType k) { return m_CellData[getCellLinearizedIndex(i, j, k)]; }
CellData& getCellData(const IndexType i, const IndexType j, const IndexType k)
{
return m_CellData[getCellLinearizedIndex(i, j, k)];
}
///
/// \brief Get data in a cell
/// \param 3D index of a cell
///
template<class IndexType>
const CellData& getCellData(const IndexType i, const IndexType j, const IndexType k) const { return m_CellData[getCellLinearizedIndex(i, j, k)]; }
const CellData& getCellData(const IndexType i, const IndexType j, const IndexType k) const
{
return m_CellData[getCellLinearizedIndex(i, j, k)];
}
///
/// \brief Apply a function to all cell data
......
......@@ -23,4 +23,4 @@ imstk_add_library( DynamicalModels
#-----------------------------------------------------------------------------
if( iMSTK_BUILD_TESTING )
add_subdirectory( Testing )
endif()
endif()
\ No newline at end of file
#-----------------------------------------------------------------------------
# Create target
#-----------------------------------------------------------------------------
include(imstkAddLibrary)
imstk_add_library( ForceModel
DEPENDS
Core
VegaFEM::massSpringSystem
VegaFEM::corotationalLinearFEM
VegaFEM::isotropicHyperelasticFEM
VegaFEM::forceModel
VegaFEM::stvk
VegaFEM::graph
VegaFEM::volumetricMesh
)
#-----------------------------------------------------------------------------
# Testing
#-----------------------------------------------------------------------------
if( iMSTK_BUILD_TESTING )
add_subdirectory( Testing )
endif()
......@@ -31,7 +31,7 @@
namespace imstk
{
// TODO: Move to appropriate place
/// \todo Move to appropriate place
enum class ForceModelType
{
StVK,
......@@ -45,7 +45,7 @@ enum class ForceModelType
none
};
// TODO: Move to appropriate place
/// \todo Move to appropriate place
enum class HyperElasticMaterialType
{
StVK,
......@@ -91,10 +91,10 @@ enum class HyperElasticMaterialType
/// numberOfThreads Number of threads spawned by the force model
/// [default = 0]
///
/// TODO: Convert this to input through JSON format
/// \todo Convert this to input through JSON format
class ForceModelConfig
{
// TODO: Do this in a better way
/// \todo Do this in a better way
struct customOptionsList
{
char femMethod[256];
......
......@@ -36,8 +36,6 @@ enum class DynamicalModelType
elastoDynamics,
positionBasedDynamics,
SPH,
NavierStokes,
HeatEquation,
none
};
......
......@@ -250,7 +250,7 @@ FEMDeformableBodyModel::initializeMassMatrix(const bool saveToDisk /*= false*/)
this->initializeEigenMatrixFromVegaMatrix(*vegaMatrix, m_M);
// TODO: Add option to save mass matrix to file
/// \todo Add option to save mass matrix to file
return true;
}
......
......@@ -393,7 +393,8 @@ PbdModel::initializeConstantDensityConstraint(const double stiffness)
&& m_mesh->getType() != Geometry::Type::HexahedralMesh
&& m_mesh->getType() != Geometry::Type::PointSet)
{
LOG(WARNING) << "Constant constraint should come with a mesh"; //TODO: Really only need a point cloud, so may need to change this.
//\todo Really only need a point cloud, so may need to change this.
LOG(WARNING) << "Constant constraint should come with a mesh!";
return false;
}
......
......@@ -432,7 +432,8 @@ protected:
} // end namespace SPH
///
/// \brief Struct contains SPH kernels for time integration, using different kernel for different purposes
/// \brief Class contains SPH kernels for time integration,
/// using different kernel for different purposes
///
class SPHSimulationKernels
{
......
......@@ -27,51 +27,59 @@ namespace imstk
{
SPHModelConfig::SPHModelConfig(const Real particleRadius)
{
m_ParticleRadius = particleRadius;
if (std::abs(particleRadius) > Real(1.e-6))
{
LOG_IF(WARNING, (particleRadius < 0)) << "Particle radius supplied is negative! Using absolute value of the supplied radius.";
m_particleRadius = std::abs(particleRadius);
}
else
{
LOG(WARNING) << "Particle radius too small! Setting to 1.e-6";
m_particleRadius = 1.e-6;
}
initialize();
}
void
SPHModelConfig::initialize()
{
LOG_IF(FATAL, (std::abs(m_ParticleRadius) < Real(1e-6))) << "Particle radius was not set properly";
m_ParticleRadiusSqr = m_ParticleRadius * m_ParticleRadius;
// Compute the derived quantities
m_particleRadiusSqr = m_particleRadius * m_particleRadius;
m_ParticleMass = Real(std::pow(Real(2.0) * m_ParticleRadius, 3)) * m_RestDensity * m_ParticleMassScale;
m_RestDensitySqr = m_RestDensity * m_RestDensity;
m_RestDensityInv = Real(1) / m_RestDensity;
m_particleMass = Real(std::pow(Real(2.0) * m_particleRadius, 3)) * m_restDensity * m_particleMassScale;
m_restDensitySqr = m_restDensity * m_restDensity;
m_restDensityInv = Real(1) / m_restDensity;
m_KernelRadius = m_ParticleRadius * m_RatioKernelOverParticleRadius;
m_KernelRadiusSqr = m_KernelRadius * m_KernelRadius;
m_kernelRadius = m_particleRadius * m_kernelOverParticleRadiusRatio;
m_kernelRadiusSqr = m_kernelRadius * m_kernelRadius;
}
bool
SPHModel::initialize()
{
LOG_IF(FATAL, (!m_Geometry)) << "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>();
this->m_currentState = std::make_shared<SPHKinematicState>();
// Set particle positions and zero default velocities
// TODO: set particle data with given (non-zero) velocities
this->m_initialState->setParticleData(m_Geometry->getVertexPositions());
/// \todo set particle data with given (non-zero) velocities
this->m_initialState->setParticleData(m_geometry->getVertexPositions());
this->m_currentState->setState(this->m_initialState);
// Attach current state to simulation state
m_SimulationState.setKinematicState(this->m_currentState);
m_simulationState.setKinematicState(this->m_currentState);
// Initialize (allocate memory for) simulation data such as density, acceleration etc.
m_SimulationState.initializeData();
m_simulationState.initializeData();
// Initialize simulation dependent parameters and kernel data
m_Kernels.initialize(m_Parameters->m_KernelRadius);
m_kernels.initialize(m_modelParameters->m_kernelRadius);
// Initialize neighbor searcher
m_NeighborSearcher = std::make_shared<NeighborSearch>(m_Parameters->m_NeighborSearchMethod,
m_Parameters->m_KernelRadius);
m_neighborSearcher = std::make_shared<NeighborSearch>(m_modelParameters->m_NeighborSearchMethod,
m_modelParameters->m_kernelRadius);
return true;
}
......@@ -95,7 +103,7 @@ SPHModel::advanceTimeStep()
void
SPHModel::computeTimeStepSize()
{
m_dt = (this->m_timeStepSizeType == TimeSteppingType::fixed) ? m_DefaultDt : computeCFLTimeStepSize();
m_dt = (this->m_timeStepSizeType == TimeSteppingType::fixed) ? m_defaultDt : computeCFLTimeStepSize();
}
Real
......@@ -104,16 +112,18 @@ SPHModel::computeCFLTimeStepSize()
auto maxVel = ParallelUtils::findMaxL2Norm(getState().getVelocities());
// dt = CFL * 2r / max{|| v ||}
Real timestep = maxVel > Real(1e-6) ? m_Parameters->m_CFLFactor * (Real(2.0) * m_Parameters->m_ParticleRadius / maxVel) : m_Parameters->m_MaxTimestep;
Real timestep = maxVel > Real(1e-6) ?
m_modelParameters->m_CFLFactor * (Real(2.0) * m_modelParameters->m_particleRadius / maxVel) :
m_modelParameters->m_maxTimestep;
// clamp the time step size to be within a given range
if (timestep > m_Parameters->m_MaxTimestep)
if (timestep > m_modelParameters->m_maxTimestep)
{
timestep = m_Parameters->m_MaxTimestep;
timestep = m_modelParameters->m_maxTimestep;
}
else if (timestep < m_Parameters->m_MinTimestep)
else if (timestep < m_modelParameters->m_minTimestep)
{
timestep = m_Parameters->m_MinTimestep;
timestep = m_modelParameters->m_minTimestep;
}
return timestep;
......@@ -122,11 +132,12 @@ SPHModel::computeCFLTimeStepSize()
void
SPHModel::findParticleNeighbors()
{
m_NeighborSearcher->getNeighbors(getState().getFluidNeighborLists(), getState().getPositions());
if (m_Parameters->m_bDensityWithBoundary) // if considering boundary particles for computing fluid density
m_neighborSearcher->getNeighbors(getState().getFluidNeighborLists(), getState().getPositions());
if (m_modelParameters->m_bDensityWithBoundary) // if considering boundary particles for computing fluid density
{
m_NeighborSearcher->getNeighbors(getState().getBoundaryNeighborLists(), getState().getPositions(),
getState().getBoundaryParticlePositions());
m_neighborSearcher->getNeighbors(getState().getBoundaryNeighborLists(),
getState().getPositions(),
getState().getBoundaryParticlePositions());
}
}
......@@ -139,7 +150,7 @@ SPHModel::computeNeighborRelativePositions()
{
const Vec3r& qpos = allPositions[q];
const Vec3r r = ppos - qpos;
neighborInfo.push_back({ r, m_Parameters->m_RestDensity });
neighborInfo.push_back({ r, m_modelParameters->m_restDensity });
}
};
......@@ -152,7 +163,7 @@ SPHModel::computeNeighborRelativePositions()
computeRelativePositions(ppos, getState().getFluidNeighborLists()[p], getState().getPositions(), neighborInfo);
// if considering boundary particles then also cache relative positions with them
if (m_Parameters->m_bDensityWithBoundary)
if (m_modelParameters->m_bDensityWithBoundary)
{
computeRelativePositions(ppos, getState().getBoundaryNeighborLists()[p], getState().getBoundaryParticlePositions(), neighborInfo);
}
......@@ -163,7 +174,7 @@ void
SPHModel::collectNeighborDensity()
{
// after computing particle densities, cache them into neighborInfo variable, next to relative positions
// this is usefull because relative positions and densities are accessed together multiple times
// this is useful because relative positions and densities are accessed together multiple times
// caching relative positions and densities therefore can reduce computation time significantly (tested)
ParallelUtils::parallelFor(getState().getNumParticles(),
[&](const size_t p) {
......@@ -196,9 +207,9 @@ SPHModel::computeDensity()
Real pdensity = 0;
for (const auto& qInfo : neighborInfo)
{
pdensity += m_Kernels.W(qInfo.xpq);
pdensity += m_kernels.W(qInfo.xpq);
}
pdensity *= m_Parameters->m_ParticleMass;
pdensity *= m_modelParameters->m_particleMass;
getState().getDensities()[p] = pdensity;
});
}
......@@ -206,7 +217,7 @@ SPHModel::computeDensity()
void
SPHModel::normalizeDensity()
{
if (!m_Parameters->m_bNormalizeDensity)
if (!m_modelParameters->m_bNormalizeDensity)
{
return;
}
......@@ -230,10 +241,10 @@ SPHModel::normalizeDensity()
// because we're not done with density computation, qInfo does not contain desity of particle q yet
const auto q = fluidNeighborList[i];
const auto qdensity = getState().getDensities()[q];
tmp += m_Kernels.W(qInfo.xpq) / qdensity;
tmp += m_kernels.W(qInfo.xpq) / qdensity;
}
if (m_Parameters->m_bDensityWithBoundary)
if (m_modelParameters->m_bDensityWithBoundary)
{
const auto& BDNeighborList = getState().getBoundaryNeighborLists()[p];
#if defined(DEBUG) || defined(_DEBUG) || !defined(NDEBUG)
......@@ -243,11 +254,11 @@ SPHModel::normalizeDensity()
for (size_t i = fluidNeighborList.size(); i < neighborInfo.size(); ++i)
{
const auto& qInfo = neighborInfo[i];
tmp += m_Kernels.W(qInfo.xpq) / m_Parameters->m_RestDensity; // density of boundary particle is set to rest density
tmp += m_kernels.W(qInfo.xpq) / m_modelParameters->m_restDensity; // density of boundary particle is set to rest density
}
}
getState().getNormalizedDensities()[p] = getState().getDensities()[p] / (tmp * m_Parameters->m_ParticleMass);
getState().getNormalizedDensities()[p] = getState().getDensities()[p] / (tmp * m_modelParameters->m_particleMass);
});
// put normalized densities to densities
......@@ -258,7 +269,7 @@ void
SPHModel::computePressureAcceleration()
{
auto particlePressure = [&](const Real density) {
const Real error = std::pow(density / m_Parameters->m_RestDensity, 7) - Real(1);
const Real error = std::pow(density / m_modelParameters->m_restDensity, 7) - Real(1);
// clamp pressure error to zero to maintain stability
return error > Real(0) ? error : Real(0);
};
......@@ -284,10 +295,10 @@ SPHModel::computePressureAcceleration()
const auto qpressure = particlePressure(qdensity);
// pressure forces
accel += -(ppressure / (pdensity * pdensity) + qpressure / (qdensity * qdensity)) * m_Kernels.gradW(r);
accel += -(ppressure / (pdensity * pdensity) + qpressure / (qdensity * qdensity)) * m_kernels.gradW(r);
}
accel *= m_Parameters->m_PressureStiffness * m_Parameters->m_ParticleMass;
accel *= m_modelParameters->m_pressureStiffness * m_modelParameters->m_particleMass;
getState().getAccelerations()[p] = accel;
});
}
......@@ -297,7 +308,7 @@ SPHModel::updateVelocity(Real timestep)
{
ParallelUtils::parallelFor(getState().getNumParticles(),
[&](const size_t p) {
getState().getVelocities()[p] += (m_Parameters->m_Gravity + getState().getAccelerations()[p]) * timestep;
getState().getVelocities()[p] += (m_modelParameters->m_gravity + getState().getAccelerations()[p]) * timestep;
});
}
......@@ -324,23 +335,23 @@ SPHModel::computeViscosity()
const auto& qInfo = neighborInfo[i];
const auto r = qInfo.xpq;
const auto qdensity = qInfo.density;
diffuseFluid += (Real(1.0) / qdensity) * m_Kernels.W(r) * (qvel - pvel);
diffuseFluid += (Real(1.0) / qdensity) * m_kernels.W(r) * (qvel - pvel);
}
diffuseFluid *= m_Parameters->m_ViscosityFluid;
diffuseFluid *= m_modelParameters->m_viscosityCoeff;
Vec3r diffuseBoundary(0, 0, 0);
if (m_Parameters->m_bDensityWithBoundary)
if (m_modelParameters->m_bDensityWithBoundary)
{
for (size_t i = fluidNeighborList.size(); i < neighborInfo.size(); ++i)
{
const auto& qInfo = neighborInfo[i];
const auto r = qInfo.xpq;
diffuseBoundary -= m_Parameters->m_RestDensityInv * m_Kernels.W(r) * pvel;
diffuseBoundary -= m_modelParameters->m_restDensityInv * m_kernels.W(r) * pvel;
}
diffuseBoundary *= m_Parameters->m_ViscosityBoundary;
diffuseBoundary *= m_modelParameters->m_viscosityBoundary;
}
getState().getDiffuseVelocities()[p] = (diffuseFluid + diffuseBoundary) * m_Parameters->m_ParticleMass;
getState().getDiffuseVelocities()[p] = (diffuseFluid + diffuseBoundary) * m_modelParameters->m_particleMass;
});
// add diffused velocity back to velocity, causing viscosity
......@@ -350,11 +361,10 @@ SPHModel::computeViscosity()
});
}
// Compute surface tension Akinci et at. 2013 model (Versatile Surface Tension and Adhesion for SPH Fluids)
void
SPHModel::computeSurfaceTension()
{
// Firstly compute surface normal for all particles
// First, compute surface normal for all particles
ParallelUtils::parallelFor(getState().getNumParticles(),
[&](const size_t p) {
Vec3r n(0, 0, 0);
......@@ -370,14 +380,14 @@ SPHModel::computeSurfaceTension()
const auto& qInfo = neighborInfo[i];
const auto r = qInfo.xpq;
const auto qdensity = qInfo.density;
n += (Real