From cc409317c50b34a9af08e8abb22a1b8da0076cb3 Mon Sep 17 00:00:00 2001 From: David Wells Date: Thu, 6 Jul 2023 09:49:02 -0400 Subject: [PATCH] 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 #include #include #include 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(t1 - t0).count() << std::endl; } --- include/deal.II/grid/tria_accessor.h | 5 + source/grid/tria_accessor.cc | 191 ++++++++------------------- 2 files changed, 63 insertions(+), 133 deletions(-) diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h index e21afd1182..65e1a02047 100644 --- 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 5319a70503..98602affe2 100644 --- 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 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 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) { -- 2.39.5