From c3a71712e6de22a24804a356aacbd0c1ec43aea0 Mon Sep 17 00:00:00 2001 From: David Wells Date: Tue, 1 Sep 2015 17:51:37 -0400 Subject: [PATCH] Use exceptions to mark Q1 transform failures. The previous implementation returned magical numbers in some instances. --- source/fe/mapping_q1.cc | 53 +++++++++++++++++++++++++---------------- 1 file changed, 32 insertions(+), 21 deletions(-) diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc index 98287bf2a9..faefe2be6c 100644 --- a/source/fe/mapping_q1.cc +++ b/source/fe/mapping_q1.cc @@ -103,11 +103,12 @@ namespace internal const double c = (x0 - x1)*y - (x - x1)*y0 + (x - x0)*y1; const double discriminant = b*b - 4*a*c; - // fast exit if the point is not in the cell (this is the only case - // where the discriminant is negative) + // exit if the point is not in the cell (this is the only case where the + // discriminant is negative) if (discriminant < 0.0) { - return Point<2>(2, 2); + AssertThrow (false, + (typename Mapping::ExcTransformationFailed())); } double eta1; @@ -160,9 +161,15 @@ namespace internal } else // give up and try Newton iteration { - return Point<2>(2, 2); + AssertThrow (false, + (typename Mapping::ExcTransformationFailed())); } } + // bogus return to placate compiler. It should not be possible to get + // here. + Assert(false, ExcInternalError()); + return Point<2>(std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN()); } @@ -421,29 +428,33 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it // algorithm is used. const std_cxx11::array, GeometryInfo::vertices_per_cell> vertices = this->get_vertices(cell); - // These internal routines do not throw exceptions when the point does - // not lie in the cell because manually checking (see the tolerances - // before the return statements) is about 10% faster than a throw/catch. - Point point = internal::MappingQ1::transform_real_to_unit_cell(vertices, p); - - if (dim == 1) - { - // formula not subject to any issues - return point; - } - else if (dim == 2) + try { - // formula not guaranteed to work for points outside of the cell - const double eps = 1e-15; - if (-eps <= point(1) && point(1) <= 1 + eps && - -eps <= point(0) && point(0) <= 1 + eps) + Point point = internal::MappingQ1::transform_real_to_unit_cell(vertices, p); + + if (dim == 1) { + // formula not subject to any issues return point; } + else if (dim == 2) + { + // formula not guaranteed to work for points outside of the cell + const double eps = 1e-15; + if (-eps <= point(1) && point(1) <= 1 + eps && + -eps <= point(0) && point(0) <= 1 + eps) + { + return point; + } + } + else + { + Assert(false, ExcInternalError()); + } } - else + catch (const typename Mapping::ExcTransformationFailed &) { - Assert(false, ExcInternalError()); + // continue on to the standard Newton code } } // Find the initial value for the -- 2.39.5