From f39ba236493968e4718023cdb7c6fdfac32cd388 Mon Sep 17 00:00:00 2001 From: David Wells Date: Tue, 25 Jul 2017 17:44:25 -0400 Subject: [PATCH] Rearrange the MappingQ1 quadratic equation solver. This gets around an issue where, for parallelograms, a should be zero: the previous version did not handle this correctly. --- source/fe/mapping_q_generic.cc | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 5c3feeed41..94699c45e2 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -123,19 +123,20 @@ namespace internal eta1 = -c/b; eta2 = -c/b; } - // special case #2: if c is very small or the square root of the - // discriminant is nearly b. - else if (std::abs(c) < 1e-12*std::abs(b) - || std::abs(std::sqrt(discriminant) - b) <= 1e-14*std::abs(b)) + // special case #2: a is zero for parallelograms and very small for + // near-parallelograms: + else if (std::abs(a) < 1e-8*std::abs(b)) { - eta1 = (-b - std::sqrt(discriminant)) / (2*a); - eta2 = (-b + std::sqrt(discriminant)) / (2*a); + // if both a and c are very small then the root should be near + // zero: this first case will capture that + eta1 = 2*c / (-b - std::sqrt(discriminant)); + eta2 = 2*c / (-b + std::sqrt(discriminant)); } - // finally, use the numerically stable version of the quadratic formula: + // finally, use the plain version: else { - eta1 = 2*c / (-b - std::sqrt(discriminant)); - eta2 = 2*c / (-b + std::sqrt(discriminant)); + eta1 = (-b - std::sqrt(discriminant)) / (2*a); + eta2 = (-b + std::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; -- 2.39.5