Skip to content

GitLab

Commit e318c359 authored by Sreekanth Arikatla's avatar Sreekanth Arikatla
Browse files

Merge branch 'BendConstraintStride' into 'master' 

ENH: BendConstraintStride and PbdConstraintFunctor warnings fixed

See merge request !672
parents 5e7c7499 a55301f1
Pipeline #247044 failed with stage
in 0 seconds
Hide whitespace changes
Inline Side-by-side
Source/Constraint/Testing/imstkPbdBendConstraintTest.cpp 0 → 100644
View file @ e318c359
/*=========================================================================
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 "gtest/gtest.h"
#include "imstkPbdBendConstraint.h"
using namespace imstk;
///
/// \brief Test that two connecting line segments unfold
///
TEST(imstkPbdBendConstraintTest, TestConvergence1)
{
PbdBendConstraint constraint;
// Straight line upon initialization
VecDataArray<double, 3> vertices(3);
vertices[0] = Vec3d(0.0, 0.0, 0.0);
vertices[1] = Vec3d(0.5, 0.0, 0.0);
vertices[2] = Vec3d(1.0, 0.0, 0.0);
DataArray<double> invMasses(3);
invMasses[0] = 1.0;
invMasses[1] = 0.0; // Center doesn't move
invMasses[2] = 1.0;
constraint.initConstraint(vertices, 0, 1, 2, 1e20);
// Modify it so the line segments look like \/
vertices[0][1] = 0.1;
vertices[2][1] = 0.1;
for (int i = 0; i < 500; i++)
{
constraint.projectConstraint(invMasses, 0.01, PbdConstraint::SolverType::xPBD, vertices);
}
// Should resolve back to a flat line
EXPECT_NEAR(vertices[0][1], 0.0, IMSTK_DOUBLE_EPS);
EXPECT_NEAR(vertices[2][1], 0.0, IMSTK_DOUBLE_EPS);
}
\ No newline at end of file
Source/Constraint/Testing/imstkPbdPointEdgeConstraintTest.cpp
View file @ e318c359
......@@ -25,34 +25,27 @@
using namespace imstk;
///
/// \brief TODO
///
class imstkPbdPointEdgeConstraintTest : public ::testing::Test
{
protected:
PbdPointEdgeConstraint m_constraint;
};
///
/// \brief Test that a point and edge meet on touching
///
TEST_F(imstkPbdPointEdgeConstraintTest, TestConvergence1)
TEST(imstkPbdPointEdgeConstraintTest, TestConvergence1)
{
PbdPointEdgeConstraint constraint;
Vec3d a = Vec3d(-0.5, 0.0, 0.0);
Vec3d b = Vec3d(0.5, 0.0, 0.0);
Vec3d x = Vec3d(0.0, -1.0, 0.0);
Vec3d zeroVelocity = Vec3d::Zero();
m_constraint.initConstraint(
constraint.initConstraint(
{ &x, 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
m_constraint.solvePosition();
constraint.solvePosition();
}
EXPECT_NEAR(x[1], a[1], IMSTK_DOUBLE_EPS);
......@@ -63,22 +56,24 @@ TEST_F(imstkPbdPointEdgeConstraintTest, TestConvergence1)
///
/// \brief Test that a point and edge meet on touching
///
TEST_F(imstkPbdPointEdgeConstraintTest, TestConvergence2)
TEST(imstkPbdPointEdgeConstraintTest, TestConvergence2)
{
PbdPointEdgeConstraint constraint;
Vec3d a = Vec3d(-0.5, 0.0, 0.0);
Vec3d b = Vec3d(0.5, 0.0, 0.0);
Vec3d x = Vec3d(0.0, 1.0, 0.0);
Vec3d zeroVelocity = Vec3d::Zero();
m_constraint.initConstraint(
constraint.initConstraint(
{ &x, 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
m_constraint.solvePosition();
constraint.solvePosition();
}
EXPECT_NEAR(x[1], a[1], IMSTK_DOUBLE_EPS);
......@@ -89,8 +84,10 @@ TEST_F(imstkPbdPointEdgeConstraintTest, TestConvergence2)
///
/// \brief Test that a point not within bounds of edge does not move (left)
///
TEST_F(imstkPbdPointEdgeConstraintTest, TestNoConvergence1)
TEST(imstkPbdPointEdgeConstraintTest, TestNoConvergence1)
{
PbdPointEdgeConstraint constraint;
Vec3d a = Vec3d(-0.5, 0.0, 0.0);
const Vec3d aInit = a;
Vec3d b = Vec3d(0.5, 0.0, 0.0);
......@@ -100,14 +97,14 @@ TEST_F(imstkPbdPointEdgeConstraintTest, TestNoConvergence1)
const Vec3d xInit = x;
Vec3d zeroVelocity = Vec3d::Zero();
m_constraint.initConstraint(
constraint.initConstraint(
{ &x, 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
m_constraint.solvePosition();
constraint.solvePosition();
}
EXPECT_EQ(xInit[0], x[0]);
......@@ -126,8 +123,10 @@ TEST_F(imstkPbdPointEdgeConstraintTest, TestNoConvergence1)
///
/// \brief Test that a point not within bounds of edge does not move (right)
///
TEST_F(imstkPbdPointEdgeConstraintTest, TestNoConvergence2)
TEST(imstkPbdPointEdgeConstraintTest, TestNoConvergence2)
{
PbdPointEdgeConstraint constraint;
Vec3d a = Vec3d(-0.5, 0.0, 0.0);
const Vec3d aInit = a;
Vec3d b = Vec3d(0.5, 0.0, 0.0);
......@@ -137,14 +136,14 @@ TEST_F(imstkPbdPointEdgeConstraintTest, TestNoConvergence2)
const Vec3d xInit = x;
Vec3d zeroVelocity = Vec3d::Zero();
m_constraint.initConstraint(
constraint.initConstraint(
{ &x, 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
m_constraint.solvePosition();
constraint.solvePosition();
}
EXPECT_EQ(xInit[0], x[0]);
......@@ -158,17 +157,4 @@ TEST_F(imstkPbdPointEdgeConstraintTest, TestNoConvergence2)
EXPECT_EQ(bInit[0], b[0]);
EXPECT_EQ(bInit[1], b[1]);
EXPECT_EQ(bInit[2], b[2]);
}
///
/// \brief TODO
///
int
imstkPbdPointEdgeConstraintTest(int argc, char* argv[])
{
// Init Google Test & Mock
::testing::InitGoogleTest(&argc, argv);
// Run tests with gtest
return RUN_ALL_TESTS();
}
}
\ No newline at end of file
Source/Constraint/Testing/imstkPbdPointPointConstraintTest.cpp
View file @ e318c359
......@@ -25,44 +25,24 @@
using namespace imstk;
///
/// \brief TODO
///
class imstkPbdPointPointConstraintTest : public ::testing::Test
{
protected:
PbdPointPointConstraint m_constraint;
};
///
/// \brief Test that two points meet
///
TEST_F(imstkPbdPointPointConstraintTest, TestConvergence1)
TEST(imstkPbdPointPointConstraintTest, TestConvergence1)
{
PbdPointPointConstraint constraint;
Vec3d a = Vec3d(0.0, 0.0, 0.0);
Vec3d b = Vec3d(0.0, -1.0, 0.0);
m_constraint.initConstraint(
constraint.initConstraint(
{ &a, 1.0, nullptr },
{ &b, 1.0, nullptr },
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
m_constraint.solvePosition();
constraint.solvePosition();
}
ASSERT_EQ(a[1], b[1]);
}
///
/// \brief TODO
///
int
imstkPbdPointPointConstraintTest(int argc, char* argv[])
{
// Init Google Test & Mock
::testing::InitGoogleTest(&argc, argv);
// Run tests with gtest
return RUN_ALL_TESTS();
}
}
\ No newline at end of file
Source/Constraint/Testing/imstkPbdPointTriangleConstraintTest.cpp
View file @ e318c359
......@@ -25,20 +25,13 @@
using namespace imstk;
///
/// \brief TODO
///
class imstkPbdPointTriangleConstraintTest : public ::testing::Test
{
protected:
PbdPointTriangleConstraint m_constraint;
};
///
/// \brief Test that a point below a triangle, and the triangle, meet on y axes
///
TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence1)
TEST(imstkPbdPointTriangleConstraintTest, TestConvergence1)
{
PbdPointTriangleConstraint constraint;
Vec3d a = Vec3d(0.5, 0.0, -0.5);
Vec3d b = Vec3d(-0.5, 0.0, -0.5);
Vec3d c = Vec3d(0.0, 0.0, 0.5);
......@@ -47,7 +40,7 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence1)
x[1] -= 1.0;
Vec3d zeroVelocity = Vec3d::Zero();
m_constraint.initConstraint(
constraint.initConstraint(
{ &x, 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
......@@ -55,7 +48,7 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence1)
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
m_constraint.solvePosition();
constraint.solvePosition();
}
EXPECT_NEAR(x[1], a[1], 0.00000001);
......@@ -66,8 +59,10 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence1)
///
/// \brief Test that a point above a triangle, and the triangle, meet on y axes
///
TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence2)
TEST(imstkPbdPointTriangleConstraintTest, TestConvergence2)
{
PbdPointTriangleConstraint constraint;
Vec3d a = Vec3d(0.5, 0.0, -0.5);
Vec3d b = Vec3d(-0.5, 0.0, -0.5);
Vec3d c = Vec3d(0.0, 0.0, 0.5);
......@@ -76,7 +71,7 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence2)
x[1] += 1.0;
Vec3d zeroVelocity = Vec3d::Zero();
m_constraint.initConstraint(
constraint.initConstraint(
{ &x, 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
......@@ -84,7 +79,7 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence2)
1.0, 1.0);
for (int i = 0; i < 3; i++)
{
m_constraint.solvePosition();
constraint.solvePosition();
}
EXPECT_NEAR(x[1], a[1], 0.00000001);
......@@ -95,8 +90,10 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestConvergence2)
///
/// \brief Test that a point not within the triangle does not move at all
///
TEST_F(imstkPbdPointTriangleConstraintTest, TestNoConvergence1)
TEST(imstkPbdPointTriangleConstraintTest, TestNoConvergence1)
{
PbdPointTriangleConstraint constraint;
Vec3d a = Vec3d(0.5, 0.0, -0.5);
const Vec3d aInit = a;
Vec3d b = Vec3d(-0.5, 0.0, -0.5);
......@@ -120,7 +117,7 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestNoConvergence1)
c = cInit;
Vec3d zeroVelocity = Vec3d::Zero();
m_constraint.initConstraint(
constraint.initConstraint(
{ &testPts[i], 1.0, &zeroVelocity },
{ &a, 1.0, &zeroVelocity },
{ &b, 1.0, &zeroVelocity },
......@@ -128,7 +125,7 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestNoConvergence1)
1.0, 1.0);
for (int j = 0; j < 3; j++)
{
m_constraint.solvePosition();
constraint.solvePosition();
}
// Test they haven't moved
......@@ -148,17 +145,4 @@ TEST_F(imstkPbdPointTriangleConstraintTest, TestNoConvergence1)
EXPECT_EQ(cInit[1], c[1]);
EXPECT_EQ(cInit[2], c[2]);
}
}
///
/// \brief TODO
///
int
imstkPbdPointTriangleConstraintTest(int argc, char* argv[])
{
// Init Google Test & Mock
::testing::InitGoogleTest(&argc, argv);
// Run tests with gtest
return RUN_ALL_TESTS();
}
}
\ No newline at end of file
Source/DynamicalModels/CMakeLists.txt
View file @ e318c359
......@@ -57,6 +57,6 @@ imstk_add_library(DynamicalModels
#-----------------------------------------------------------------------------
# Testing
#-----------------------------------------------------------------------------
#if( ${PROJECT_NAME}_BUILD_TESTING )
# add_subdirectory( Testing )
#endif()
if( ${PROJECT_NAME}_BUILD_TESTING )
add_subdirectory(Testing)
endif()
Source/DynamicalModels/ObjectModels/imstkPbdConstraintFunctor.h
View file @ e318c359
......@@ -50,6 +50,7 @@ struct PbdConstraintFunctor
PbdConstraintFunctor() = default;
virtual ~PbdConstraintFunctor() = default;
public:
///
/// \brief Appends a set of constraint to the container given a geometry
///
......@@ -58,7 +59,7 @@ struct PbdConstraintFunctor
void setGeometry(std::shared_ptr<PointSet> geom) { m_geom = geom; }
public:
std::shared_ptr<PointSet> m_geom;
std::shared_ptr<PointSet> m_geom = nullptr;
};
struct PbdDistanceConstraintFunctor : public PbdConstraintFunctor
......@@ -67,10 +68,12 @@ struct PbdDistanceConstraintFunctor : public PbdConstraintFunctor
PbdDistanceConstraintFunctor() = default;
~PbdDistanceConstraintFunctor() override = default;
public:
virtual void operator()(PbdConstraintContainer& constraints) override
{
const VecDataArray<double, 3>& vertices = *m_geom->getVertexPositions();
auto addDistConstraint =
std::shared_ptr<VecDataArray<double, 3>> verticesPtr = m_geom->getVertexPositions();
const VecDataArray<double, 3>& vertices = *verticesPtr;
auto addDistConstraint =
[&](std::vector<std::vector<bool>>& E, size_t i1, size_t i2)
{
if (i1 > i2) // Make sure i1 is always smaller than i2
......@@ -134,10 +137,10 @@ struct PbdDistanceConstraintFunctor : public PbdConstraintFunctor
}
}
void setStiffness(double stiffness) { m_stiffness = stiffness; }
void setStiffness(const double stiffness) { m_stiffness = stiffness; }
protected:
double m_stiffness;
double m_stiffness = 0.0;
};
struct PbdFemConstraintFunctor : public PbdConstraintFunctor
......@@ -146,6 +149,7 @@ struct PbdFemConstraintFunctor : public PbdConstraintFunctor
PbdFemConstraintFunctor() = default;
~PbdFemConstraintFunctor() override = default;
public:
virtual void operator()(PbdConstraintContainer& constraints) override
{
// Check if constraint type matches the mesh type
......@@ -153,9 +157,11 @@ struct PbdFemConstraintFunctor : public PbdConstraintFunctor
<< "FEM Tetrahedral constraint should come with tetrahedral mesh";
// Create constraints
auto tetMesh = std::dynamic_pointer_cast<TetrahedralMesh>(m_geom);
const VecDataArray<double, 3>& vertices = *m_geom->getVertexPositions();
const VecDataArray<int, 4>& elements = *tetMesh->getTetrahedraIndices();
auto tetMesh = std::dynamic_pointer_cast<TetrahedralMesh>(m_geom);
std::shared_ptr<VecDataArray<double, 3>> verticesPtr = m_geom->getVertexPositions();
const VecDataArray<double, 3>& vertices = *verticesPtr;
std::shared_ptr<VecDataArray<int, 4>> elementsPtr = tetMesh->getTetrahedraIndices();
const VecDataArray<int, 4>& elements = *elementsPtr;
ParallelUtils::parallelFor(elements.size(),
[&](const size_t k)
......@@ -168,12 +174,12 @@ struct PbdFemConstraintFunctor : public PbdConstraintFunctor
}, elements.size() > 100);
}
void setMaterialType(PbdFEMTetConstraint::MaterialType materialType) { m_matType = materialType; }
void setMaterialType(const PbdFEMTetConstraint::MaterialType materialType) { m_matType = materialType; }
void setFemConfig(std::shared_ptr<PbdFEMConstraintConfig> femConfig) { m_femConfig = femConfig; }
protected:
PbdFEMTetConstraint::MaterialType m_matType;
std::shared_ptr<PbdFEMConstraintConfig> m_femConfig;
PbdFEMTetConstraint::MaterialType m_matType = PbdFEMTetConstraint::MaterialType::StVK;
std::shared_ptr<PbdFEMConstraintConfig> m_femConfig = nullptr;
};
struct PbdVolumeConstraintFunctor : public PbdConstraintFunctor
......@@ -182,6 +188,7 @@ struct PbdVolumeConstraintFunctor : public PbdConstraintFunctor
PbdVolumeConstraintFunctor() = default;
~PbdVolumeConstraintFunctor() override = default;
public:
virtual void operator()(PbdConstraintContainer& constraints) override
{
// Check if constraint type matches the mesh type
......@@ -189,9 +196,11 @@ struct PbdVolumeConstraintFunctor : public PbdConstraintFunctor
<< "Volume constraint should come with volumetric mesh";
// Create constraints
auto tetMesh = std::dynamic_pointer_cast<TetrahedralMesh>(m_geom);
const VecDataArray<double, 3>& vertices = *m_geom->getVertexPositions();
const VecDataArray<int, 4>& elements = *tetMesh->getTetrahedraIndices();
auto tetMesh = std::dynamic_pointer_cast<TetrahedralMesh>(m_geom);
std::shared_ptr<VecDataArray<double, 3>> verticesPtr = m_geom->getVertexPositions();
const VecDataArray<double, 3>& vertices = *verticesPtr;
std::shared_ptr<VecDataArray<int, 4>> elementsPtr = tetMesh->getTetrahedraIndices();
const VecDataArray<int, 4>& elements = *elementsPtr;
ParallelUtils::parallelFor(elements.size(),
[&](const size_t k)
......@@ -204,10 +213,10 @@ struct PbdVolumeConstraintFunctor : public PbdConstraintFunctor
});
}
void setStiffness(double stiffness) { m_stiffness = stiffness; }
void setStiffness(const double stiffness) { m_stiffness = stiffness; }
protected:
double m_stiffness;
double m_stiffness = 0.0;
};
struct PbdAreaConstraintFunctor : public PbdConstraintFunctor
......@@ -216,6 +225,7 @@ struct PbdAreaConstraintFunctor : public PbdConstraintFunctor
PbdAreaConstraintFunctor() = default;
~PbdAreaConstraintFunctor() override = default;
public:
virtual void operator()(PbdConstraintContainer& constraints) override
{
// check if constraint type matches the mesh type
......@@ -223,9 +233,11 @@ struct PbdAreaConstraintFunctor : public PbdConstraintFunctor
<< "Area constraint should come with a triangular mesh";
// ok, now create constraints
auto triMesh = std::dynamic_pointer_cast<SurfaceMesh>(m_geom);
const VecDataArray<double, 3>& vertices = *m_geom->getVertexPositions();
const VecDataArray<int, 3>& elements = *triMesh->getTriangleIndices();
auto triMesh = std::dynamic_pointer_cast<SurfaceMesh>(m_geom);
std::shared_ptr<VecDataArray<double, 3>> verticesPtr = m_geom->getVertexPositions();
const VecDataArray<double, 3>& vertices = *verticesPtr;
std::shared_ptr<VecDataArray<int, 3>> elemenstPtr = triMesh->getTriangleIndices();
const VecDataArray<int, 3>& elements = *elemenstPtr;
ParallelUtils::parallelFor(elements.size(),
[&](const size_t k)
......@@ -237,10 +249,10 @@ struct PbdAreaConstraintFunctor : public PbdConstraintFunctor
});
}
void setStiffness(double stiffness) { m_stiffness = stiffness; }
void setStiffness(const double stiffness) { m_stiffness = stiffness; }
protected:
double m_stiffness;
double m_stiffness = 0.0;
};
struct PbdBendConstraintFunctor : public PbdConstraintFunctor
......@@ -249,14 +261,17 @@ struct PbdBendConstraintFunctor : public PbdConstraintFunctor
PbdBendConstraintFunctor() = default;
~PbdBendConstraintFunctor() override = default;
public:
virtual void operator()(PbdConstraintContainer& constraints) override
{
CHECK(m_geom->getTypeName() == "LineMesh")
<< "Bend constraint should come with a line mesh";
auto lineMesh = std::dynamic_pointer_cast<LineMesh>(m_geom);
const VecDataArray<double, 3>& vertices = *m_geom->getVertexPositions();
const VecDataArray<int, 2>& elements = *lineMesh->getLinesIndices();
auto lineMesh = std::dynamic_pointer_cast<LineMesh>(m_geom);
std::shared_ptr<VecDataArray<double, 3>> verticesPtr = m_geom->getVertexPositions();
const VecDataArray<double, 3>& vertices = *verticesPtr;
/*std::shared_ptr< VecDataArray<int, 2>> indicesPtr = lineMesh->getLinesIndices();
const VecDataArray<int, 2>& indices = *indicesPtr*/
auto addBendConstraint =
[&](const double k, size_t i1, size_t i2, size_t i3)
......@@ -277,24 +292,23 @@ struct PbdBendConstraintFunctor : public PbdConstraintFunctor
constraints.addConstraint(c);
};
// Iterate sets of two segments
for (int k = 0; k < elements.size() - 1; k++)
// Iterate sets of stride # of segments
for (int k = 0; k < vertices.size() - m_stride * 2; k += m_stride)
{
auto& seg1 = elements[k];
auto& seg2 = elements[k + 1];
int i3 = seg2[0];
if (i3 == seg1[0] || i3 == seg1[1])
{
i3 = seg2[1];
}
addBendConstraint(m_stiffness, seg1[0], seg1[1], i3);
addBendConstraint(m_stiffness, k, k + m_stride, k + 2 * m_stride);
}
}
void setStiffness(double stiffness) { m_stiffness = stiffness; }
void setStiffness(const double stiffness) { m_stiffness = stiffness; }
void setStride(const int stride)
{
CHECK(m_stride > 1) << "Stride should be at least 1.";
m_stride = stride;
}
protected:
double m_stiffness;
double m_stiffness = 0.0;
int m_stride = 3;
};
struct PbdDihedralConstraintFunctor : public PbdConstraintFunctor
......@@ -303,17 +317,20 @@ struct PbdDihedralConstraintFunctor : public PbdConstraintFunctor
PbdDihedralConstraintFunctor() = default;
~PbdDihedralConstraintFunctor() override = default;
public:
virtual void operator()(PbdConstraintContainer& constraints) override