Skip to content

Commit

Permalink
#123 Added functionality to shift ellipse
Browse files Browse the repository at this point in the history
  • Loading branch information
@pou036
pou036 committed Mar 21, 2018
1 parent f438eb3 commit 34460e5
Show file tree
Hide file tree
Showing 4 changed files with 52 additions and 55 deletions.
3 changes: 2 additions & 1 deletion include/utils/Ellipse.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,9 +101,10 @@ class Ellipse
* @param y: coordinates of point to project
* @param x: coordinates of projection point on the ellipse
* @param s: arc-length from point to yield point
* @param shift: horizontal shift of the ellipse
* @return d: distance to ellipse
*/ static void
getYieldPointCC(Real const m, Real const p_c, Real const y0, Real const y1, Real & x0, Real & x1, Real & s);
getYieldPointCC(Real const m, Real const p_c, Real const y0, Real const y1, Real & x0, Real & x1, Real & s, Real shift);

/**
* Function to compute the distance from a query point to an ellipse
Expand Down
2 changes: 1 addition & 1 deletion src/materials/RedbackMechMaterialCC.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,7 +179,7 @@ RedbackMechMaterialCC::get_py_qy(Real p, Real q, Real & p_y, Real & q_y, Real yi

// get yield point in any case (even if elastic)
//Ellipse::distanceCC(_slope_yield_surface, -yield_stress, p, q, p_y, q_y, _shift_ellipse);
Ellipse::getYieldPointCC(_slope_yield_surface, -yield_stress, p, q, p_y, q_y, s);
Ellipse::getYieldPointCC(_slope_yield_surface, -yield_stress, p, q, p_y, q_y, s, _shift_ellipse);
}

void
Expand Down
39 changes: 17 additions & 22 deletions src/utils/Ellipse.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,64 +171,59 @@ Ellipse::distanceCC(Real const m, Real const p_c, Real const y0, Real const y1,
}

void
Ellipse::getYieldPointCC(Real const m, Real const p_c, Real const y0, Real const y1, Real & x0, Real & x1, Real & s)
Ellipse::getYieldPointCC(Real const m, Real const p_c, Real const y0, Real const y1, Real & x0, Real & x1, Real & s, Real shift)
{
if (y1 < 0)
{
// algorithm below only works for y1 > 0, so use symmetry
Real neg_x1;
getYieldPointCC(m, p_c, y0, -y1, x0, neg_x1, s);
getYieldPointCC(m, p_c, y0, -y1, x0, neg_x1, s, shift);
x1 = -neg_x1;
return;
}
// Deal with easy cases first
if (y0 == p_c / 2.0)
if (y0 == p_c / 2.0 + shift)
{
// Real t = -std::pow(m, 2)*std::log(-m*p_c/(2.0*y1))/2.0;
x0 = p_c / 2.0;
x0 = y0;
x1 = -m * p_c / 2.0;
s = std::abs(y1 - x1);
return;
}
else if (y1 == 0)
{
// Real t = -std::log(std::abs(p_c/(2*y0-p_c)))/2.0;
if (y0 > p_c / 2.0)
{
x0 = 0.0;
s = std::abs(y0);
}
if (y0 > p_c / 2.0 + shift)
x0 = shift;
else
{
x0 = p_c;
s = std::abs(y0 - x0);
}
x0 = p_c + shift;
s = std::abs(y0 - x0);
x1 = 0.0;
return;
}
// For all other cases, we use Newton Raphson
Real tol = 1e-15; // tolerance on ellipse potential to be close to 0
int nb_iter_max = 1000; // number of dichotomies to avoid infinite loops
int nb_iter_max = 1000; // max number of iterations to avoid infinite loops
Real m2 = m * m;
Real t = 0;
Real phi = std::pow(y1 / m, 2) + std::pow(y0 - p_c / 2.0, 2) - p_c * p_c / 4.0;
Real phi_prime = -4 * std::pow(y1 / m2, 2) - 4 * std::pow(y0 - p_c / 2.0, 2);
Real phi = std::pow(y1 / m, 2) + std::pow(y0 - shift - p_c / 2.0, 2) - p_c * p_c / 4.0;
Real phi_prime = -4 * std::pow(y1 / m2, 2) - 4 * std::pow(y0 - shift - p_c / 2.0, 2);
int nb_iter = 0;
while (std::abs(phi) > tol && nb_iter < nb_iter_max)
{
nb_iter += 1;
t -= phi / phi_prime;
phi =
std::pow(y1 / m, 2) * std::exp(-4 * t / m2) + std::pow(y0 - p_c / 2.0, 2) * std::exp(-4 * t) - p_c * p_c / 4.0;
phi_prime = -4 * std::pow(y1 / m2, 2) * std::exp(-4 * t / m2) - 4 * std::pow(y0 - p_c / 2.0, 2) * std::exp(-4 * t);
std::pow(y1 / m, 2) * std::exp(-4 * t / m2) + std::pow(y0 - shift - p_c / 2.0, 2) * std::exp(-4 * t) - p_c * p_c / 4.0;
phi_prime = -4 * std::pow(y1 / m2, 2) * std::exp(-4 * t / m2) - 4 * std::pow(y0 - shift - p_c / 2.0, 2) * std::exp(-4 * t);
}
if (nb_iter == nb_iter_max)
mooseError("Newton Raphson failed to converge after ", nb_iter, " iterations.");
x0 = p_c / 2.0 + (y0 - p_c / 2.0) * std::exp(-2 * t);
mooseError("Newton Raphson (getYieldPointCC) failed to converge after ", nb_iter, " iterations.");
x0 = p_c / 2.0 + shift + (y0 - shift - p_c / 2.0) * std::exp(-2 * t);
x1 = y1 * std::exp(-2 * t / m2);

// Compute arc-length from point to yield point
/*// If we have access to hypergeometric function 2F1:
/*// If we have access to hypergeometric function 2F1 (not accounting for shift!)
gamma = p - pc/2.0
alpha = q*q/(M*M*M*M*gamma*gamma)
beta = 2*(1 - 1/(M*M))
Expand All @@ -243,7 +238,7 @@ Ellipse::getYieldPointCC(Real const m, Real const p_c, Real const y0, Real const
for (int i=1; i<n_iter+1; i++)
{
t_new = i*t/n_iter;
s += std::sqrt( std::pow((y0 - p_c/2.0)*(std::exp(-2*t_new) - std::exp(-2*t_old)), 2) \
s += std::sqrt( std::pow((y0 - shift - p_c/2.0)*(std::exp(-2*t_new) - std::exp(-2*t_old)), 2) \
+ std::pow(y1*(std::exp(-2*t_new/m2) - std::exp(-2*t_old/m2)), 2) );
t_old = t_new;
}
Expand Down
63 changes: 32 additions & 31 deletions unit/src/EllipseTest.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -589,18 +589,19 @@ EllipseTest::isPointOutsideOfEllipseTestMajorAxisVertical()
void
EllipseTest::distanceCCTestShift()
{
Real d, x2, y2;//, x3, y3;
Real d, x2, y2, x3, y3, s;
Real shift = 0.5;

// Case: circle, p_c < 0
// top right quadrant
d = Ellipse::distanceCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/2.0, /*y1=*/3.0, x2, y2, shift);
//Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/2.0, /*y1=*/3.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/2.0, /*y1=*/3.0, x3, y3, s, shift);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3 * std::sqrt(2) - 1.5, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.5 / std::sqrt(2) - 1.0, x2, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.5 / std::sqrt(2), y2, 1e-5);
//CPPUNIT_ASSERT_DOUBLES_EQUAL(1.5 / std::sqrt(2) - 1.5, x3, 1e-5);
//CPPUNIT_ASSERT_DOUBLES_EQUAL(1.5 / std::sqrt(2), y3, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3 * std::sqrt(2) - 1.5, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.5 / std::sqrt(2) - 1.0, x2, 1e-14);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.5 / std::sqrt(2), y2, 1e-14);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.5 / std::sqrt(2) - 1.0, x3, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.5 / std::sqrt(2), y3, 1e-5);
}

/**
Expand All @@ -615,7 +616,7 @@ EllipseTest::distanceCCTestCircle()
// Case: circle, p_c < 0
// top right quadrant
d = Ellipse::distanceCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/1.5, /*y1=*/3.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/1.5, /*y1=*/3.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/1.5, /*y1=*/3.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3 * std::sqrt(2) - 1.5, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3 * std::sqrt(2) - 1.5, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.5 / std::sqrt(2) - 1.5, x2, 1e-5);
Expand All @@ -624,7 +625,7 @@ EllipseTest::distanceCCTestCircle()
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.5 / std::sqrt(2), y3, 1e-5);
// top left quadrant
d = Ellipse::distanceCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-4.5, /*y1=*/3.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-4.5, /*y1=*/3.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-4.5, /*y1=*/3.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3 * std::sqrt(2) - 1.5, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3 * std::sqrt(2) - 1.5, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.5 / std::sqrt(2) - 1.5, x2, 1e-5);
Expand All @@ -633,7 +634,7 @@ EllipseTest::distanceCCTestCircle()
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.5 / std::sqrt(2), y3, 1e-5);
// bottom right quadrant
d = Ellipse::distanceCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/1.5, /*y1=*/-3.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/1.5, /*y1=*/-3.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/1.5, /*y1=*/-3.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3 * std::sqrt(2) - 1.5, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3 * std::sqrt(2) - 1.5, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.5 / std::sqrt(2) - 1.5, x2, 1e-5);
Expand All @@ -642,7 +643,7 @@ EllipseTest::distanceCCTestCircle()
CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.5 / std::sqrt(2), y3, 1e-5);
// bottom left quadrant
d = Ellipse::distanceCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-4.5, /*y1=*/-3.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-4.5, /*y1=*/-3.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-4.5, /*y1=*/-3.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3 * std::sqrt(2) - 1.5, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3 * std::sqrt(2) - 1.5, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.5 / std::sqrt(2) - 1.5, x2, 1e-5);
Expand All @@ -651,7 +652,7 @@ EllipseTest::distanceCCTestCircle()
CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.5 / std::sqrt(2), y3, 1e-5);
// on horizontal axis, left hand side of ellipse
d = Ellipse::distanceCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-4.0, /*y1=*/0.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-4.0, /*y1=*/0.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-4.0, /*y1=*/0.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-3.0, x2, 1e-5);
Expand All @@ -660,7 +661,7 @@ EllipseTest::distanceCCTestCircle()
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, y3, 1e-5);
// on horizontal axis, right hand side of ellipse
d = Ellipse::distanceCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/1.0, /*y1=*/0.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/1.0, /*y1=*/0.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/1.0, /*y1=*/0.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, x2, 1e-5);
Expand All @@ -669,7 +670,7 @@ EllipseTest::distanceCCTestCircle()
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, y3, 1e-5);
// on vertical axis, below ellipse
d = Ellipse::distanceCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/-2.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/-2.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/-2.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.5, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.5, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.5, x2, 1e-5);
Expand All @@ -678,7 +679,7 @@ EllipseTest::distanceCCTestCircle()
CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.5, y3, 1e-5);
// on vertical axis, above ellipse
d = Ellipse::distanceCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/2.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/2.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.0, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/2.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.5, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.5, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.5, x2, 1e-5);
Expand All @@ -699,7 +700,7 @@ EllipseTest::distanceCCTestMajorAxisHorizontal()
// Case: ellipse, major axis horizontal, p_c < 0
// top right quadrant
d = Ellipse::distanceCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/1.5, /*y1=*/3.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/1.5, /*y1=*/3.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/1.5, /*y1=*/3.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3.1210069, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3.31752131594, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.282058, x2, 1e-5);
Expand All @@ -708,7 +709,7 @@ EllipseTest::distanceCCTestMajorAxisHorizontal()
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.168952819574, y3, 1e-10);
// top left quadrant
d = Ellipse::distanceCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-4.5, /*y1=*/3.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-4.5, /*y1=*/3.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-4.5, /*y1=*/3.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3.1210069, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3.31752131594, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-2.717941, x2, 1e-5);
Expand All @@ -717,7 +718,7 @@ EllipseTest::distanceCCTestMajorAxisHorizontal()
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.168952819574, y3, 1e-10);
// bottom right quadrant
d = Ellipse::distanceCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/1.5, /*y1=*/-3.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/1.5, /*y1=*/-3.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/1.5, /*y1=*/-3.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3.1210069, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3.31752131594, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.282058, x2, 1e-5);
Expand All @@ -726,7 +727,7 @@ EllipseTest::distanceCCTestMajorAxisHorizontal()
CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.168952819574, y3, 1e-10);
// bottom left quadrant
d = Ellipse::distanceCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-4.5, /*y1=*/-3.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-4.5, /*y1=*/-3.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-4.5, /*y1=*/-3.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3.1210069, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(3.31752131594, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-2.717941, x2, 1e-5);
Expand All @@ -735,7 +736,7 @@ EllipseTest::distanceCCTestMajorAxisHorizontal()
CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.168952819574, y3, 1e-10);
// on horizontal axis, left hand side of ellipse
d = Ellipse::distanceCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-4.0, /*y1=*/0.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-4.0, /*y1=*/0.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-4.0, /*y1=*/0.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-3.0, x2, 1e-5);
Expand All @@ -744,7 +745,7 @@ EllipseTest::distanceCCTestMajorAxisHorizontal()
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, y3, 1e-10);
// on horizontal axis, right hand side of ellipse
d = Ellipse::distanceCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/1.0, /*y1=*/0.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/1.0, /*y1=*/0.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/1.0, /*y1=*/0.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, x2, 1e-5);
Expand All @@ -753,7 +754,7 @@ EllipseTest::distanceCCTestMajorAxisHorizontal()
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, y3, 1e-10);
// on vertical axis, below ellipse
d = Ellipse::distanceCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/-1.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/-1.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/-1.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.5, x2, 1e-5);
Expand All @@ -762,7 +763,7 @@ EllipseTest::distanceCCTestMajorAxisHorizontal()
CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.75, y3, 1e-10);
// on vertical axis, above ellipse
d = Ellipse::distanceCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/1.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/1.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/0.5, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/1.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.5, x2, 1e-5);
Expand All @@ -778,12 +779,12 @@ EllipseTest::distanceCCTestMajorAxisHorizontal()
void
EllipseTest::distanceCCTestMajorAxisVertical()
{
Real d, x2, y2, x3, y3;
Real d, x2, y2, x3, y3, s;

// Case: ellipse, major axis vertical, p_c < 0
// top right quadrant
d = Ellipse::distanceCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-0.5, /*y1=*/2.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-0.5, /*y1=*/2.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-0.5, /*y1=*/2.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.2094051, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.209475269026, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.650317, x2, 1e-5);
Expand All @@ -792,7 +793,7 @@ EllipseTest::distanceCCTestMajorAxisVertical()
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.85734897716, y3, 1e-10);
// top left quadrant
d = Ellipse::distanceCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-2.5, /*y1=*/2.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-2.5, /*y1=*/2.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-2.5, /*y1=*/2.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.2094051, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.209475269026, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-2.349682, x2, 1e-5);
Expand All @@ -801,7 +802,7 @@ EllipseTest::distanceCCTestMajorAxisVertical()
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.85734897716, y3, 1e-10);
// bottom right quadrant
d = Ellipse::distanceCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-0.5, /*y1=*/-2.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-0.5, /*y1=*/-2.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-0.5, /*y1=*/-2.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.2094051, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.209475269026, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.650317, x2, 1e-5);
Expand All @@ -810,7 +811,7 @@ EllipseTest::distanceCCTestMajorAxisVertical()
CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.85734897716, y3, 1e-10);
// bottom left quadrant
d = Ellipse::distanceCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-2.5, /*y1=*/-2.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-2.5, /*y1=*/-2.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-2.5, /*y1=*/-2.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.2094051, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.209475269026, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-2.349682, x2, 1e-5);
Expand All @@ -820,7 +821,7 @@ EllipseTest::distanceCCTestMajorAxisVertical()

// on horizontal axis, left hand side of ellipse
d = Ellipse::distanceCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-4.0, /*y1=*/0.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-4.0, /*y1=*/0.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-4.0, /*y1=*/0.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-3.0, x2, 1e-5);
Expand All @@ -829,7 +830,7 @@ EllipseTest::distanceCCTestMajorAxisVertical()
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, y3, 1e-10);
// on horizontal axis, right hand side of ellipse
d = Ellipse::distanceCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/1.0, /*y1=*/0.0, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/1.0, /*y1=*/0.0, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/1.0, /*y1=*/0.0, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, x2, 1e-5);
Expand All @@ -838,7 +839,7 @@ EllipseTest::distanceCCTestMajorAxisVertical()
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, y3, 1e-10);
// on vertical axis, below ellipse
d = Ellipse::distanceCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/-2.5, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/-2.5, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/-2.5, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.5, x2, 1e-5);
Expand All @@ -847,7 +848,7 @@ EllipseTest::distanceCCTestMajorAxisVertical()
CPPUNIT_ASSERT_DOUBLES_EQUAL(-2.25, y3, 1e-10);
// on vertical axis, above ellipse
d = Ellipse::distanceCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/2.5, x2, y2, 0.);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/2.5, x3, y3, s);
Ellipse::getYieldPointCC(/*m=*/1.5, /*p_c=*/-3.0, /*x1=*/-1.5, /*y1=*/2.5, x3, y3, s, 0.);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25, d, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25, s, 1e-5);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-1.5, x2, 1e-5);
Expand Down

0 comments on commit 34460e5

Please sign in to comment.