]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix singularity (FPE) in MappingQ1::real_to_unit_cell.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 1 Dec 2017 08:41:59 +0000 (09:41 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 1 Dec 2017 08:41:59 +0000 (09:41 +0100)
source/fe/mapping_q_generic.cc

index 678348928dd55f0112ecaeaac0852d5eb9e420b3..1534005c440cda984eb2907fe1cdf1a9a5e23e50 100644 (file)
@@ -95,8 +95,8 @@ namespace internal
        const Point<spacedim> &p)
       {
         Assert(spacedim == 2, ExcInternalError());
-        const long double x = p(0);
-        const long double y = p(1);
+        const double x = p(0);
+        const double y = p(1);
 
         const double x0 = vertices[0](0);
         const double x1 = vertices[1](0);
@@ -108,21 +108,23 @@ namespace internal
         const double y2 = vertices[2](1);
         const double y3 = vertices[3](1);
 
-        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 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 discriminant = b*b - 4*a*c;
+        const 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()));
 
-        long double eta1;
-        long double eta2;
-        // special case #1: if a is zero, then use the linear formula
-        if (a == 0.0 && b != 0.0)
+        double eta1;
+        double eta2;
+        const 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)
           {
             eta1 = -c/b;
             eta2 = -c/b;
@@ -133,41 +135,41 @@ namespace internal
           {
             // 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));
+            eta1 = 2*c / (-b - sqrt_discriminant);
+            eta2 = 2*c / (-b + sqrt_discriminant);
           }
         // finally, use the plain version:
         else
           {
-            eta1 = (-b - std::sqrt(discriminant)) / (2*a);
-            eta2 = (-b + std::sqrt(discriminant)) / (2*a);
+            eta1 = (-b - sqrt_discriminant) / (2*a);
+            eta2 = (-b + sqrt_discriminant) / (2*a);
           }
         // pick the one closer to the center of the cell.
-        const long double eta = (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2;
+        const 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 long double subexpr0 = -eta*x2 + x0*(eta - 1);
-        const long double xi_denominator0 = eta*x3 - x1*(eta - 1) + subexpr0;
+        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)));
 
         if (std::abs(xi_denominator0) > 1e-10*max_x)
           {
-            const long double xi = (x + subexpr0)/xi_denominator0;
+            const double xi = (x + subexpr0)/xi_denominator0;
             return Point<2>(xi, eta);
           }
         else
           {
             const 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;
+            const double subexpr1 = -eta*y2 + y0*(eta - 1);
+            const double xi_denominator1 = eta*y3 - y1*(eta - 1) + subexpr1;
             if (std::abs(xi_denominator1) > 1e-10*max_y)
               {
-                const long double xi = (subexpr1 + y)/xi_denominator1;
+                const double xi = (subexpr1 + y)/xi_denominator1;
                 return Point<2>(xi, eta);
               }
             else // give up and try Newton iteration

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.