Redo the quadrilateral measure function.

This is a lot simpler and somewhat more efficient:
master, best of 10, summing hyper ball face measures: 2434 ms
feature, same setup: 2124 ms

Lets not use complex checks for planarity - either we have a parallelogram (up
to the last few significant digits) or we do not. If people are reading in CAD
data then they probably don't have uniform meshes anyway.

P.S.

benchmark source:

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>

#include <chrono>
#include <iostream>

int
main()
{
  using namespace dealii;

  Triangulation<3> tria;
  GridGenerator::hyper_ball(tria);

  tria.refine_global(3);

  double area = 0.0;
  auto   t0   = std::chrono::high_resolution_clock::now();
  for (unsigned int i = 0; i < 1000; ++i)
    for (const auto &face : tria.active_face_iterators())
      area += face->measure();
  auto t1 = std::chrono::high_resolution_clock::now();
  std::cout
    << "time elapsed = "
    << std::chrono::duration_cast<std::chrono::milliseconds>(t1 - t0).count()
    << std::endl;
}

diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h
index e21afd1182987858111d12f1649ee21cfe8cbfed..65e1a020473b1c719f00b3f18fe224f86c3acd97 100644 (file)
--- a/include/deal.II/grid/tria_accessor.h
+++ b/include/deal.II/grid/tria_accessor.h
@@ -1697,6 +1697,11 @@ public:
    * after applying quadrature equals the summing the JxW values returned by
    * the FEValues or FEFaceValues object you will want to use for the
    * integral.
+   *
+   * @note There is no analytic formula for the area of a bilinear face (i.e.,
+   * something with a quadrilateral reference cell) in 3D. This function uses
+   * the quadrature defined by QGauss<2>(4) to approximate the measure in this
+   * case.
    */
   double
   measure() const;

diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc
index 5319a70503852adfa0b34dec549d542958462620..98602affe26a40f393187caeef8a2924ca3288fb 100644 (file)
--- a/source/grid/tria_accessor.cc
+++ b/source/grid/tria_accessor.cc
@@ -1329,139 +1329,64 @@ namespace
   {
     if (accessor.reference_cell() == ReferenceCells::Quadrilateral)
       {
-        // If the face is planar, the diagonal from vertex 0 to vertex 3,
-        // v_03, should be in the plane P_012 of vertices 0, 1 and 2.  Get
-        // the normal vector of P_012 and test if v_03 is orthogonal to
-        // that. If so, the face is planar and computing its area is simple.
-        const Tensor<1, 3> v01 = accessor.vertex(1) - accessor.vertex(0);
-        const Tensor<1, 3> v02 = accessor.vertex(2) - accessor.vertex(0);
-
-        const Tensor<1, 3> normal = cross_product_3d(v01, v02);
-
-        const Tensor<1, 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 * v03) * (v01 * v01) * (v02 * v02))) >= 1e-24)
-          {
-            // If the vectors are non planar we integrate the norm of the normal
-            // vector using a numerical Gauss scheme of order 4. In particular
-            // we consider a bilinear quad x(u,v) = (1-v)((1-u)v_0 + u v_1) +
-            // v((1-u)v_2 + u v_3), consequently we compute the normal vector as
-            // n(u,v) = t_u x t_v = w_1 + u w_2 + v w_3. The integrand function
-            // is
-            // || n(u,v) || = sqrt(a + b u^2 + c v^2 + d u + e v + f uv).
-            // We integrate it using a QGauss<2> (4) computed explicitly.
-            const Tensor<1, 3> w_1 =
-              cross_product_3d(accessor.vertex(1) - accessor.vertex(0),
-                               accessor.vertex(2) - accessor.vertex(0));
-            const Tensor<1, 3> w_2 =
-              cross_product_3d(accessor.vertex(1) - accessor.vertex(0),
-                               accessor.vertex(3) - accessor.vertex(2) -
-                                 accessor.vertex(1) + accessor.vertex(0));
-            const Tensor<1, 3> w_3 =
-              cross_product_3d(accessor.vertex(3) - accessor.vertex(2) -
-                                 accessor.vertex(1) + accessor.vertex(0),
-                               accessor.vertex(2) - accessor.vertex(0));
-
-            double a = scalar_product(w_1, w_1);
-            double b = scalar_product(w_2, w_2);
-            double c = scalar_product(w_3, w_3);
-            double d = scalar_product(w_1, w_2);
-            double e = scalar_product(w_1, w_3);
-            double f = scalar_product(w_2, w_3);
-
-            return 0.03025074832140047 *
-                     std::sqrt(
-                       a + 0.0048207809894260144 * b +
-                       0.0048207809894260144 * c + 0.13886368840594743 * d +
-                       0.13886368840594743 * e + 0.0096415619788520288 * f) +
-                   0.056712962962962937 *
-                     std::sqrt(
-                       a + 0.0048207809894260144 * b + 0.10890625570683385 * c +
-                       0.13886368840594743 * d + 0.66001895641514374 * e +
-                       0.045826333352825557 * f) +
-                   0.056712962962962937 *
-                     std::sqrt(
-                       a + 0.0048207809894260144 * b + 0.44888729929169013 * c +
-                       0.13886368840594743 * d + 1.3399810435848563 * e +
-                       0.09303735505312187 * f) +
-                   0.03025074832140047 *
-                     std::sqrt(
-                       a + 0.0048207809894260144 * b + 0.86595709258347853 * c +
-                       0.13886368840594743 * d + 1.8611363115940525 * e +
-                       0.12922212642709538 * f) +
-                   0.056712962962962937 *
-                     std::sqrt(
-                       a + 0.10890625570683385 * b + 0.0048207809894260144 * c +
-                       0.66001895641514374 * d + 0.13886368840594743 * e +
-                       0.045826333352825557 * f) +
-                   0.10632332575267359 * std::sqrt(a + 0.10890625570683385 * b +
-                                                   0.10890625570683385 * c +
-                                                   0.66001895641514374 * d +
-                                                   0.66001895641514374 * e +
-                                                   0.2178125114136677 * f) +
-                   0.10632332575267359 * std::sqrt(a + 0.10890625570683385 * b +
-                                                   0.44888729929169013 * c +
-                                                   0.66001895641514374 * d +
-                                                   1.3399810435848563 * e +
-                                                   0.44220644500147605 * f) +
-                   0.056712962962962937 *
-                     std::sqrt(
-                       a + 0.10890625570683385 * b + 0.86595709258347853 * c +
-                       0.66001895641514374 * d + 1.8611363115940525 * e +
-                       0.61419262306231814 * f) +
-                   0.056712962962962937 *
-                     std::sqrt(
-                       a + 0.44888729929169013 * b + 0.0048207809894260144 * c +
-                       1.3399810435848563 * d + 0.13886368840594743 * e +
-                       0.09303735505312187 * f) +
-                   0.10632332575267359 * std::sqrt(a + 0.44888729929169013 * b +
-                                                   0.10890625570683385 * c +
-                                                   1.3399810435848563 * d +
-                                                   0.66001895641514374 * e +
-                                                   0.44220644500147605 * f) +
-                   0.10632332575267359 *
-                     std::sqrt(a + 0.44888729929169013 * b +
-                               0.44888729929169013 * c +
-                               1.3399810435848563 * d + 1.3399810435848563 * e +
-                               0.89777459858338027 * f) +
-                   0.056712962962962937 *
-                     std::sqrt(a + 0.44888729929169013 * b +
-                               0.86595709258347853 * c +
-                               1.3399810435848563 * d + 1.8611363115940525 * e +
-                               1.2469436885317342 * f) +
-                   0.03025074832140047 * std::sqrt(a + 0.86595709258347853 * b +
-                                                   0.0048207809894260144 * c +
-                                                   1.8611363115940525 * d +
-                                                   0.13886368840594743 * e +
-                                                   0.12922212642709538 * f) +
-                   0.056712962962962937 *
-                     std::sqrt(
-                       a + 0.86595709258347853 * b + 0.10890625570683385 * c +
-                       1.8611363115940525 * d + 0.66001895641514374 * e +
-                       0.61419262306231814 * f) +
-                   0.056712962962962937 *
-                     std::sqrt(a + 0.86595709258347853 * b +
-                               0.44888729929169013 * c +
-                               1.8611363115940525 * d + 1.3399810435848563 * e +
-                               1.2469436885317342 * f) +
-                   0.03025074832140047 *
-                     std::sqrt(a + 0.86595709258347853 * b +
-                               0.86595709258347853 * c +
-                               1.8611363115940525 * d + 1.8611363115940525 * e +
-                               1.7319141851669571 * f);
-          }
-
-        // the face is planar. then its area is 1/2 of the norm of the
-        // cross product of the two diagonals
-        const Tensor<1, 3> v12        = accessor.vertex(2) - accessor.vertex(1);
-        const Tensor<1, 3> twice_area = cross_product_3d(v03, v12);
-        return 0.5 * twice_area.norm();
+        const Point<3> x0 = accessor.vertex(0);
+        const Point<3> x1 = accessor.vertex(1);
+        const Point<3> x2 = accessor.vertex(2);
+        const Point<3> x3 = accessor.vertex(3);
+
+        // This is based on the approach used in libMesh (see face_quad4.C): the
+        // primary differences are the vertex numbering and quadrature order.
+        //
+        // The area of a surface is the integral of the magnitude of its normal
+        // vector, which may be computed via the cross product of two tangent
+        // vectors. We can easily get tangent vectors from the surface
+        // parameterization. Hence, given a bilinear surface
+        //
+        //     X(chi, eta) = x0 + (x1 - x0) chi + (x2 - x0) eta
+        //                 + (x3 + x0 - x1 - x2) chi eta
+        //
+        // the tangent vectors are
+        //
+        //     t1 = (x1 - x0) + (x3 + x0 - x1 - x2) eta
+        //     t2 = (x2 - x0) + (x3 + x0 - x1 - x2) xi
+        const Tensor<1, 3> b0 = x1 - x0;
+        const Tensor<1, 3> b1 = x2 - x0;
+        const Tensor<1, 3> a  = x3 - x2 - b0;
+
+        // The diameter is the maximum distance between any pair of vertices and
+        // we can use it as a length scale for the cell. If all components of a
+        // (the vector connecting x3 and the last vertex of the parallelogram
+        // defined by the first three vertices) are zero within some tolerance,
+        // then we have a parallelogram and can use a much simpler formula.
+        double a_max = 0.0;
+        for (unsigned int d = 0; d < 3; ++d)
+          a_max = std::max(std::abs(a[d]), a_max);
+        if (a_max < 1e-14 * accessor.diameter())
+          return cross_product_3d(b0, b1).norm();
+
+        // Otherwise, use a 4x4 quadrature to approximate the surface area.
+        // Hard-code this in to prevent the extra overhead of always creating
+        // the same QGauss rule.
+        constexpr unsigned int n_qp = 4;
+        const double           c1   = 2.0 / 7.0 * std::sqrt(6.0 / 5.0);
+        const double           w0   = (18.0 - std::sqrt(30)) / 72.0;
+        const double           w1   = (18.0 + std::sqrt(30)) / 72.0;
+
+        const std::array<double, n_qp> q{{
+          0.5 - std::sqrt(3.0 / 7.0 + c1) / 2.0,
+          0.5 - std::sqrt(3.0 / 7.0 - c1) / 2.0,
+          0.5 + std::sqrt(3.0 / 7.0 - c1) / 2.0,
+          0.5 + std::sqrt(3.0 / 7.0 + c1) / 2.0,
+        }};
+        const std::array<double, n_qp> w{{w0, w1, w1, w0}};
+
+        double area = 0.;
+        for (unsigned int i = 0; i < n_qp; ++i)
+          for (unsigned int j = 0; j < n_qp; ++j)
+            area += cross_product_3d(q[i] * a + b0, q[j] * a + b1).norm() *
+                    w[i] * w[j];
+
+        return area;
       }
     else if (accessor.reference_cell() == ReferenceCells::Triangle)
       {