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.

diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc
index 5c3feeed41c6855c69325017cdd3f6790daeb2a6..94699c45e2905e388319d6e63c8c5e9f0911e6aa 100644 (file)
--- 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;