Commit 51caabd7 authored by Sreekanth Arikatla's avatar Sreekanth Arikatla
Browse files

ENH: Add support for point set based fluids

Modify PBD class to allow a point set to be used
for fluid simulations
parent 8da261b0
Pipeline #65954 passed with stage
......@@ -93,7 +93,7 @@ PbdInteractionPair::doNarrowPhaseCollision()
p2[0], p2[1], p2[2], prox1, prox2))
{
auto c = std::make_shared<PbdPointTriangleConstraint>();
c->initConstraint(dynaModel1, map1->getMapIdx(i),
c->initConstraint(dynaModel1, i,
dynaModel2, map2->getMapIdx(e[0]), map2->getMapIdx(e[1]), map2->getMapIdx(e[2]));
m_collisionConstraints.push_back(c);
}
......@@ -166,6 +166,38 @@ PbdInteractionPair::doNarrowPhaseCollision()
}
}
}
else if (g1->getType() == Geometry::Type::PointSet)
{
auto mesh1 = std::static_pointer_cast<Mesh>(g1);
// brute force, use BVH or spatial grid would be much better
// point
for (size_t i = 0; i < mesh1->getNumVertices(); ++i)
{
const auto p = mesh1->getVertexPosition(i);
auto elements = mesh2->getTrianglesVertices();
for (size_t j = 0; j < elements.size(); ++j)
{
auto& e = elements[j];
const Vec3d p0 = mesh2->getVertexPosition(e[0]);
const Vec3d p1 = mesh2->getVertexPosition(e[1]);
const Vec3d p2 = mesh2->getVertexPosition(e[2]);
if (testPointToTriAABB(p[0], p[1], p[2],
p0[0], p0[1], p0[2],
p1[0], p1[1], p1[2],
p2[0], p2[1], p2[2], prox1, prox2))
{
auto c = std::make_shared<PbdPointTriangleConstraint>();
auto mappedIndex1 = (map1) ? map1->getMapIdx(i) : i;
c->initConstraint(dynaModel1, mappedIndex1, dynaModel2, map2->getMapIdx(e[0]), map2->getMapIdx(e[1]), map2->getMapIdx(e[2]));
m_collisionConstraints.push_back(c);
}
}
}
}
else
{
auto mesh1 = std::static_pointer_cast<SurfaceMesh>(g1);
......@@ -175,7 +207,7 @@ PbdInteractionPair::doNarrowPhaseCollision()
for (size_t i = 0; i < mesh1->getNumVertices(); ++i)
{
const Vec3d p = mesh1->getVertexPosition(i);
std::vector<SurfaceMesh::TriangleArray> elements = mesh2->getTrianglesVertices();
auto elements = mesh2->getTrianglesVertices();
for (size_t j = 0; j < elements.size(); ++j)
{
......
......@@ -41,10 +41,6 @@ PbdConstantDensityConstraint::initConstraint(PbdModel& model, const double k)
m_numNeighbors.resize(np);
m_deltaPositions.resize(np);
m_neighbors.resize(np * m_maxNumNeighbors);
m_xPosIndexes.resize(np);
m_yPosIndexes.resize(np);
m_zPosIndexes.resize(np);
}
bool
......@@ -101,7 +97,7 @@ PbdConstantDensityConstraint::gradSpiky(const Vec3d &pi, const Vec3d &pj)
}
inline void
PbdConstantDensityConstraint::clearNeighbors(const size_t &np)
PbdConstantDensityConstraint::clearNeighbors(const size_t np)
{
m_numNeighbors.clear();
m_neighbors.clear();
......@@ -111,7 +107,7 @@ PbdConstantDensityConstraint::clearNeighbors(const size_t &np)
inline void
PbdConstantDensityConstraint::updateNeighborsBruteForce(const Vec3d &pi,
const size_t &index,
const size_t index,
const StdVectorOfVec3d &positions)
{
const double neighborRadius = m_maxDist;
......@@ -141,7 +137,7 @@ PbdConstantDensityConstraint::updateNeighborsBruteForce(const Vec3d &pi,
inline void
PbdConstantDensityConstraint::calculateDensityEstimate(const Vec3d &pi,
const size_t &index,
const size_t index,
const StdVectorOfVec3d &positions)
{
double densitySum = 0.0;
......@@ -155,7 +151,7 @@ PbdConstantDensityConstraint::calculateDensityEstimate(const Vec3d &pi,
inline void
PbdConstantDensityConstraint::calculateLambdaScalingFactor(const Vec3d &pi,
const size_t &index,
const size_t index,
const StdVectorOfVec3d &positions)
{
const double densityConstraint = (m_densities[index] / m_restDensity) - 1;
......@@ -170,7 +166,7 @@ PbdConstantDensityConstraint::calculateLambdaScalingFactor(const Vec3d &pi,
inline void
PbdConstantDensityConstraint::updatePositions(const Vec3d &pi,
const size_t &index,
const size_t index,
StdVectorOfVec3d &positions)
{
//Make sure the point is valid
......
......@@ -37,7 +37,7 @@ public:
///
/// \brief constructor
///
PbdConstantDensityConstraint() : PbdConstraint() { m_vertexIds.resize(1); }
PbdConstantDensityConstraint() : PbdConstraint() {}
///
/// \Constant Density Constraint Initialization
......@@ -73,27 +73,27 @@ private:
///
/// \brief Update the neighbors of each node using burte force search O(n*n)
///
void updateNeighborsBruteForce(const Vec3d &pi, const size_t &index, const StdVectorOfVec3d &positions);
void updateNeighborsBruteForce(const Vec3d &pi, const size_t index, const StdVectorOfVec3d &positions);
///
/// \brief Clear the list of neighbors
///
void clearNeighbors(const size_t &np);
void clearNeighbors(const size_t np);
///
/// \brief
///
void calculateDensityEstimate(const Vec3d &pi, const size_t &index, const StdVectorOfVec3d &positions);
void calculateDensityEstimate(const Vec3d &pi, const size_t index, const StdVectorOfVec3d &positions);
///
/// \brief
///
void calculateLambdaScalingFactor(const Vec3d &pi, const size_t &index, const StdVectorOfVec3d &positions);
void calculateLambdaScalingFactor(const Vec3d &pi, const size_t index, const StdVectorOfVec3d &positions);
///
/// \brief
///
void updatePositions(const Vec3d &pi, const size_t &index, StdVectorOfVec3d &positions);
void updatePositions(const Vec3d &pi, const size_t index, StdVectorOfVec3d &positions);
///
/// \brief Set/Get rest density
......@@ -121,10 +121,6 @@ private:
std::vector<Vec3d> m_deltaPositions; ///> delta positions
std::vector<int> m_neighbors; ///> index of neighbors
std::vector<int> m_numNeighbors; ///> number of neighbors
std::vector<int> m_xPosIndexes;
std::vector<int> m_yPosIndexes;
std::vector<int> m_zPosIndexes;
};
}
......
......@@ -96,7 +96,7 @@ PbdPointTriangleConstraint::solvePositionConstraint()
lambda = (l - dist)/lambda;
// LOG(INFO) << "Lambda:" << lambda <<" Normal:" << n[0] <<" " << n[1] <<" "<<n[2];
//LOG(INFO) << "Lambda:" << lambda <<" Normal:" << n[0] <<" " << n[1] <<" "<<n[2];
if (im0 > 0)
{
......
......@@ -401,7 +401,7 @@ PbdModel::initializeConstantDensityConstraint(const double& stiffness)
{
// check if constraint type matches the mesh type
if (m_mesh->getType() != Geometry::Type::SurfaceMesh && m_mesh->getType() != Geometry::Type::TetrahedralMesh && m_mesh->getType() != Geometry::Type::LineMesh &&
m_mesh->getType() != Geometry::Type::HexahedralMesh)
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.
return false;
......
......@@ -219,10 +219,10 @@ protected:
double m_lambda; ///> Lame constant
// Mass properties
std::vector<double> m_mass; ///> Mass of nodes
std::vector<double> m_invMass; ///> Inverse of mass of nodes
std::vector<double> m_mass; ///> Mass of nodes
std::vector<double> m_invMass; ///> Inverse of mass of nodes
double m_contactStiffness; ///> Contact stiffness for collisions
double m_contactStiffness = 1.; ///> Contact stiffness for collisions
Vec3d m_gravity; ///> Gravity
double m_viscousDampingCoeff = 0.01; ///> Viscous damping coefficient [0, 1]
......
......@@ -167,7 +167,7 @@ PbdObject::initialize(int nCons, ...)
m_pbdModel->setMaxNumIterations(va_arg(args,int));
}
if (m_physicsToCollidingGeomMap && m_collidingGeometry)
if (m_collidingGeometry)
{
m_pbdModel->setProximity(va_arg(args,double));
m_pbdModel->setContactStiffness(va_arg(args,double));
......
Supports Markdown
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