const Point<spacedim> &p)
{
Assert(spacedim == 2, ExcInternalError());
- const double x = p(0);
- const double y = p(1);
- const double x0 = vertices[0](0);
- const double x1 = vertices[1](0);
- const double x2 = vertices[2](0);
- const double x3 = vertices[3](0);
+ // For accuracy reasons, we do all arithmetics in extended precision
+ // (long double). This has a noticable effect on the hit rate for
+ // borderline cases and thus makes the algorithm more robust.
+ const long double x = p(0);
+ const long double y = p(1);
- const double y0 = vertices[0](1);
- const double y1 = vertices[1](1);
- const double y2 = vertices[2](1);
- const double y3 = vertices[3](1);
+ const long double x0 = vertices[0](0);
+ const long double x1 = vertices[1](0);
+ const long double x2 = vertices[2](0);
+ const long double x3 = vertices[3](0);
- const double a = (x1 - x3)*(y0 - y2) - (x0 - x2)*(y1 - y3);
- const double b = -(x0 - x1 - x2 + x3)*y + (x - 2*x1 + x3)*y0 - (x - 2*x0 + x2)*y1
- - (x - x1)*y2 + (x - x0)*y3;
- const double c = (x0 - x1)*y - (x - x1)*y0 + (x - x0)*y1;
+ const long double y0 = vertices[0](1);
+ const long double y1 = vertices[1](1);
+ const long double y2 = vertices[2](1);
+ const long double y3 = vertices[3](1);
- const double discriminant = b*b - 4*a*c;
+ const long double a = (x1 - x3)*(y0 - y2) - (x0 - x2)*(y1 - y3);
+ const long double b = -(x0 - x1 - x2 + x3)*y + (x - 2*x1 + x3)*y0 - (x - 2*x0 + x2)*y1
+ - (x - x1)*y2 + (x - x0)*y3;
+ const long double c = (x0 - x1)*y - (x - x1)*y0 + (x - x0)*y1;
+
+ const long double discriminant = b*b - 4*a*c;
// exit if the point is not in the cell (this is the only case where the
// discriminant is negative)
AssertThrow (discriminant > 0.0,
(typename Mapping<spacedim,spacedim>::ExcTransformationFailed()));
- double eta1;
- double eta2;
- const double sqrt_discriminant = std::sqrt(discriminant);
+ long double eta1;
+ long double eta2;
+ const long double sqrt_discriminant = std::sqrt(discriminant);
// special case #1: if a is near-zero to make the discriminant exactly
// equal b, then use the linear formula
if (b != 0.0 && std::abs(b) == sqrt_discriminant)
eta2 = (-b + sqrt_discriminant) / (2*a);
}
// pick the one closer to the center of the cell.
- const double eta = (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2;
+ const long double eta = (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2;
/*
* There are two ways to compute xi from eta, but either one may have a
* zero denominator.
*/
- const double subexpr0 = -eta*x2 + x0*(eta - 1);
- const double xi_denominator0 = eta*x3 - x1*(eta - 1) + subexpr0;
- const double max_x = std::max(std::max(std::abs(x0), std::abs(x1)),
- std::max(std::abs(x2), std::abs(x3)));
+ const long double subexpr0 = -eta*x2 + x0*(eta - 1);
+ const long double xi_denominator0 = eta*x3 - x1*(eta - 1) + subexpr0;
+ const long double max_x = std::max(std::max(std::abs(x0), std::abs(x1)),
+ std::max(std::abs(x2), std::abs(x3)));
if (std::abs(xi_denominator0) > 1e-10*max_x)
{
}
else
{
- const double max_y = std::max(std::max(std::abs(y0), std::abs(y1)),
- std::max(std::abs(y2), std::abs(y3)));
- const double subexpr1 = -eta*y2 + y0*(eta - 1);
- const double xi_denominator1 = eta*y3 - y1*(eta - 1) + subexpr1;
+ const long double max_y = std::max(std::max(std::abs(y0), std::abs(y1)),
+ std::max(std::abs(y2), std::abs(y3)));
+ const long double subexpr1 = -eta*y2 + y0*(eta - 1);
+ const long double xi_denominator1 = eta*y3 - y1*(eta - 1) + subexpr1;
if (std::abs(xi_denominator1) > 1e-10*max_y)
{
const double xi = (subexpr1 + y)/xi_denominator1;