Skip to content

Commit

Permalink
Merge pull request #15755 from drwells/fix-planar-measurement
Browse files Browse the repository at this point in the history 
Redo the quadrilateral measure function.
  • Loading branch information
@tamiko
tamiko committed Jul 18, 2023
2 parents 34ae3d3 + cc40931 commit 0265947
Show file tree
Hide file tree
Showing 2 changed files with 63 additions and 133 deletions.
5 changes: 5 additions & 0 deletions include/deal.II/grid/tria_accessor.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1697,6 +1697,11 @@ class TriaAccessor : public TriaAccessorBase<structdim, dim, spacedim>
* 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;
Expand Down
191 changes: 58 additions & 133 deletions source/grid/tria_accessor.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
Expand Down

0 comments on commit 0265947

Please sign in to comment.