Further improve test to avoid square roots.

diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc
index 09008f51f8a122aef920c095498eaba246e7db35..8eb5fbae4c84d4987160504b5cc37bf9bed77760 100644 (file)
--- a/source/grid/tria_accessor.cc
+++ b/source/grid/tria_accessor.cc
@@ -983,9 +983,17 @@ namespace
     Point<3> normal;
     cross_product(normal, v01, v02);
 
-    const Point<3> v03 = accessor.vertex(3) - accessor.vertex(0); //1st diagonal
-
-    if (std::abs((v03 * normal) / (v03.norm() * v01.norm() * v02.norm())) >= 1e-12)
+    const Point<3> v03 = accessor.vertex(3) - accessor.vertex(0);
+
+    // check whether v03 does not lie in the plane of v01 and v02
+    // (i.e., whether the face is not planar). we do so by checking
+    // whether the triple product (v01 x v02) * v03 forms a positive
+    // volume relative to |v01|*|v02|*|v03|. the test checks the
+    // squares of these to avoid taking norms/square roots:
+    if (std::abs((v03 * normal) * (v03 * normal) /
+                (v03.square() * v01.square() * v02.square()))
+       >=
+       1e-24)
       {
        Assert (false,
                ExcMessage("Computing the measure of a nonplanar face is not implemented!"));