]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Revert Mapping::transform_real_to_unit_cell to long double. 5566/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 2 Dec 2017 17:30:40 +0000 (18:30 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 2 Dec 2017 17:42:18 +0000 (18:42 +0100)
source/fe/mapping_q_generic.cc

index 1534005c440cda984eb2907fe1cdf1a9a5e23e50..8fa6560d4b868c0dab8634dbce6128adb902556f 100644 (file)
@@ -95,33 +95,37 @@ namespace internal
        const Point<spacedim> &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<spacedim,spacedim>::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;

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.