From: Martin Kronbichler Date: Sat, 2 Dec 2017 17:30:40 +0000 (+0100) Subject: Revert Mapping::transform_real_to_unit_cell to long double. X-Git-Tag: v9.0.0-rc1~697^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F5566%2Fhead;p=dealii.git Revert Mapping::transform_real_to_unit_cell to long double. --- diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 1534005c44..8fa6560d4b 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -95,33 +95,37 @@ namespace internal const Point &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::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) @@ -145,16 +149,16 @@ namespace internal 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) { @@ -163,10 +167,10 @@ namespace internal } 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;