Commit fda17f2d authored by Stefan Pfeifer's avatar Stefan Pfeifer

Fix constant orientation subset method and add tests

parent acbe0763
Pipeline #40253930 passed with stage
in 9 minutes and 4 seconds
......@@ -195,6 +195,7 @@ add_executable(
tests/fem/HarmonicOscillator.cpp
tests/fem/LargeDeformationBeams.cpp
tests/fem/TangentStiffness.cpp
tests/numerics/Geometry.cpp
)
target_compile_definitions(bow-simulator-test PRIVATE _USE_MATH_DEFINES) # Defines math constants on MSVC
......
......@@ -195,7 +195,7 @@ void BowModel::init_string(const Callback& callback)
});
}
points = one_sided_orientation_subset(points, true);
points = constant_orientation_subset<Orientation::RightHanded>(points);
points = equipartition(points, input.settings.n_string_elements + 1);
// Create string nodes
......
......@@ -103,8 +103,10 @@ ContactElement::State ContactElement::get_state() const
system.get_u(dofs[4]) - h1*cos(system.get_u(dofs[5]))};
// If no contact, set kinematic expressions to zero and return
if(!is_right_handed(P2, P0, Q0) || !is_right_handed(P2, P1, P0) ||
!is_right_handed(P2, Q0, Q1) || !is_right_handed(P2, Q1, P1))
if(get_orientation(P2, P0, Q0) == Orientation::LeftHanded ||
get_orientation(P2, P1, P0) == Orientation::LeftHanded ||
get_orientation(P2, Q0, Q1) == Orientation::LeftHanded ||
get_orientation(P2, Q1, P1) == Orientation::LeftHanded)
{
return {0.0, Vector<8>::Zero(), Matrix<8>::Zero()};
}
......
......@@ -2,19 +2,32 @@
#include "RootFinding.hpp"
#include <algorithm>
// Checks if the triangle {a, b, c} is right-handed
static bool is_right_handed(const Vector<2>& a, const Vector<2>& b, const Vector<2>& c)
enum class Orientation
{
return (b[0] - a[0])*(c[1] - a[1]) - (b[1] - a[1])*(c[0] - a[0]) > 0;
LeftHanded,
RightHanded,
Collinear
};
// Calculates the orientation of the triangle {a, b, c}.
static Orientation get_orientation(const Vector<2>& a, const Vector<2>& b, const Vector<2>& c)
{
double product = (b[0] - a[0])*(c[1] - a[1]) - (b[1] - a[1])*(c[0] - a[0]);
if(product < 0.0)
return Orientation::LeftHanded;
else if(product > 0.0)
return Orientation::RightHanded;
else
return Orientation::Collinear;
}
// Traverses a curve given by a list of input points and returns a subset of points with
// strictly right handed or left handed orientation.
// strictly right handed or left handed orientation that contains the start- and endpoints.
// Implementation inspired by Graham scan for finding the convex hull of a point set
// (https://de.wikipedia.org/wiki/Graham_Scan)
//
// Todo: Find a better name. What's a curve with strictly positive/negative curvature called?
static std::vector<Vector<2>> one_sided_orientation_subset(const std::vector<Vector<2>>& input, bool right_handed)
// (https://en.wikipedia.org/wiki/Graham_scan)
template <Orientation orientation>
static std::vector<Vector<2>> constant_orientation_subset(const std::vector<Vector<2>>& input)
{
assert(input.size() >= 2);
......@@ -24,22 +37,23 @@ static std::vector<Vector<2>> one_sided_orientation_subset(const std::vector<Vec
output.push_back(input[0]);
output.push_back(input[1]);
size_t i = 2;
while(i < input.size())
for(size_t i = 2; i < input.size(); ++i)
{
auto p1 = output[output.size()-2];
auto p2 = output[output.size()-1];
auto p3 = input[i];
if(is_right_handed(p1, p2, p3) == right_handed || output.size() == 1)
{
output.push_back(p3);
++i;
}
else
auto fits_as_next = [&](Vector<2> p3)
{
if(output.size() == 1)
return true;
auto p1 = output[output.size()-2];
auto p2 = output[output.size()-1];
return get_orientation(p1, p2, p3) == orientation;
};
while(!fits_as_next(input[i]))
output.pop_back();
}
output.push_back(input[i]);
}
return output;
......
#include "numerics/Geometry.hpp"
#include <catch.hpp>
#include <iostream>
TEST_CASE("is-right-handed")
{
REQUIRE(get_orientation({0, 0}, {1, 0}, {0, 1}) == Orientation::RightHanded);
REQUIRE(get_orientation({0, 0}, {0, 1}, {1, 0}) == Orientation::LeftHanded);
REQUIRE(get_orientation({0, 0}, {1, 1}, {2, 2}) == Orientation::Collinear);
}
TEST_CASE("constant-orientation-subset")
{
std::vector<Vector<2>> input, output;
input = {{0, 0}, {1, 0}};
output = {{0, 0}, {1, 0}};
REQUIRE(constant_orientation_subset<Orientation::RightHanded>(input) == output);
input = {{0, 0}, {1, 0}, {2, 1}};
output = {{0, 0}, {1, 0}, {2, 1}};
REQUIRE(constant_orientation_subset<Orientation::RightHanded>(input) == output);
input = {{0, 0}, {1, 0}, {2, -1}};
output = {{0, 0}, {2, -1}};
REQUIRE(constant_orientation_subset<Orientation::RightHanded>(input) == output);
input = {{0, 0}, {1, 0}, {2, 1}, {3, 3}};
output = {{0, 0}, {1, 0}, {2, 1}, {3, 3}};
REQUIRE(constant_orientation_subset<Orientation::RightHanded>(input) == output);
input = {{0, 0}, {1, 0}, {2, 1}, {3, 3}, {4, 3}};
output = {{0, 0}, {1, 0}, {4, 3}};
REQUIRE(constant_orientation_subset<Orientation::RightHanded>(input) == output);
}
Markdown is supported
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