From: David Wells Date: Mon, 31 Mar 2025 13:14:02 +0000 (-0400) Subject: Implement ReferenceCell::face_measure(). X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d7c94b571fa5c287791d24bcfaab1e24c9686d01;p=dealii.git Implement ReferenceCell::face_measure(). We need these values in QProjector. --- diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index 046c926870..17a1414aac 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -681,6 +681,24 @@ public: double volume() const; + /** + * Return the $d - 1$-dimensional measure of the face of a reference cell that + * corresponds to the current object, where $d$ is the dimension of the space + * it lives in. + * + * In this context the measure of a face depends on both the type of reference + * cell as well as the face number (i.e., `*this` and @p face_no). For + * example, the measures of the sides of a ReferenceCells::Triangle are $1$, + * $\sqrt{2}$, and $1$, respectively, whereas the measures of the sides of a + * ReferenceCells::Quadrilateral are all $1$. + * + * @note Like ReferenceCell::volume(), this function defines vertices to have + * a measure of $1$ to enable the definition of one-dimensional face + * quadratures. + */ + double + face_measure(const unsigned int face_no) const; + /** * Return the barycenter (i.e., the center of mass) of the reference * cell that corresponds to the current object. The function is not @@ -3022,6 +3040,43 @@ ReferenceCell::volume() const +inline double +ReferenceCell::face_measure(const unsigned int face_no) const +{ + AssertIndexRange(face_no, n_faces()); + constexpr double X = std::numeric_limits::signaling_NaN(); + + static const std::array tri_faces{ + {1.0, std::sqrt(2.0), 1.0, X, X}}; + static const std::array tet_faces{ + {0.5, 0.5, 0.5, std::sqrt(3.0) / 2.0, X}}; + static const std::array pyramid_faces{ + {4, std::sqrt(2.0), std::sqrt(2.0), std::sqrt(2.0), std::sqrt(2.0)}}; + static const std::array wedge_faces{ + {0.5, 0.5, 1.0, std::sqrt(2.0), 1.0}}; + + switch (this->kind) + { + case ReferenceCells::Vertex: + case ReferenceCells::Line: + case ReferenceCells::Quadrilateral: + case ReferenceCells::Hexahedron: + return 1.0; + case ReferenceCells::Triangle: + return tri_faces[face_no]; + case ReferenceCells::Tetrahedron: + return tet_faces[face_no]; + case ReferenceCells::Pyramid: + return pyramid_faces[face_no]; + case ReferenceCells::Wedge: + return wedge_faces[face_no]; + } + DEAL_II_ASSERT_UNREACHABLE(); + return X; +} + + + template inline Point ReferenceCell::barycenter() const diff --git a/tests/grid/reference_cell_type_06.cc b/tests/grid/reference_cell_type_06.cc new file mode 100644 index 0000000000..6f7fd591b6 --- /dev/null +++ b/tests/grid/reference_cell_type_06.cc @@ -0,0 +1,108 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2022 - 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + +// Test ReferenceCell::face_measure() + +#include + +#include +#include + +#include +#include + +#include + +#include "../tests.h" + + +template +void +test(const ReferenceCell &reference_cell) +{ + Triangulation triangulation; + GridGenerator::reference_cell(triangulation, reference_cell); + + FE_Nothing fe(reference_cell); + + hp::QCollection quadratures; + // TODO: MappingQ asserts that face quadrature collections have exactly one + // entry. For now just special-case those. + if (reference_cell.is_hyper_cube()) + quadratures.push_back(reference_cell.face_reference_cell(0) + .template get_gauss_type_quadrature(2)); + else + for (unsigned int face_no = 0; face_no < reference_cell.n_faces(); + ++face_no) + quadratures.push_back(reference_cell.face_reference_cell(face_no) + .template get_gauss_type_quadrature(2)); + + // Set up the objects to compute an integral on the reference cell + FEFaceValues fe_face_values(fe, quadratures, update_JxW_values); + + for (unsigned int face_no = 0; face_no < reference_cell.n_faces(); ++face_no) + { + fe_face_values.reinit(triangulation.begin_active(), face_no); + + double measure = 0; + for (unsigned int i = 0; i < fe_face_values.n_quadrature_points; ++i) + measure += fe_face_values.JxW(i); + + deallog << "ReferenceCell: " << reference_cell.to_string() << std::endl; + deallog << " computed face measure = " << measure << std::endl; + deallog << " self-reported face measure = " + << reference_cell.face_measure(face_no) << std::endl; + + // TriaAccessor<0, 1, 1>::measure() does not yet exist + if constexpr (dim > 1) + { + Assert(std::fabs( + measure - + triangulation.begin_active()->face(face_no)->measure()) < + 1e-12, + ExcInternalError()); + } + Assert(std::fabs(measure - reference_cell.face_measure(face_no)) < 1e-12, + ExcInternalError()); + } +} + +int +main() +{ + initlog(); + + { + deallog.push("1D"); + test<1>(ReferenceCells::Line); + deallog.pop(); + } + + { + deallog.push("2D"); + test<2>(ReferenceCells::Quadrilateral); + test<2>(ReferenceCells::Triangle); + deallog.pop(); + } + + { + deallog.push("3D"); + test<3>(ReferenceCells::Tetrahedron); + test<3>(ReferenceCells::Pyramid); + test<3>(ReferenceCells::Wedge); + test<3>(ReferenceCells::Hexahedron); + deallog.pop(); + } +} diff --git a/tests/grid/reference_cell_type_06.output b/tests/grid/reference_cell_type_06.output new file mode 100644 index 0000000000..ff07dd1f3b --- /dev/null +++ b/tests/grid/reference_cell_type_06.output @@ -0,0 +1,88 @@ + +DEAL:1D::ReferenceCell: Line +DEAL:1D:: computed face measure = 1.00000 +DEAL:1D:: self-reported face measure = 1.00000 +DEAL:1D::ReferenceCell: Line +DEAL:1D:: computed face measure = 1.00000 +DEAL:1D:: self-reported face measure = 1.00000 +DEAL:2D::ReferenceCell: Quad +DEAL:2D:: computed face measure = 1.00000 +DEAL:2D:: self-reported face measure = 1.00000 +DEAL:2D::ReferenceCell: Quad +DEAL:2D:: computed face measure = 1.00000 +DEAL:2D:: self-reported face measure = 1.00000 +DEAL:2D::ReferenceCell: Quad +DEAL:2D:: computed face measure = 1.00000 +DEAL:2D:: self-reported face measure = 1.00000 +DEAL:2D::ReferenceCell: Quad +DEAL:2D:: computed face measure = 1.00000 +DEAL:2D:: self-reported face measure = 1.00000 +DEAL:2D::ReferenceCell: Tri +DEAL:2D:: computed face measure = 1.00000 +DEAL:2D:: self-reported face measure = 1.00000 +DEAL:2D::ReferenceCell: Tri +DEAL:2D:: computed face measure = 1.41421 +DEAL:2D:: self-reported face measure = 1.41421 +DEAL:2D::ReferenceCell: Tri +DEAL:2D:: computed face measure = 1.00000 +DEAL:2D:: self-reported face measure = 1.00000 +DEAL:3D::ReferenceCell: Tet +DEAL:3D:: computed face measure = 0.500000 +DEAL:3D:: self-reported face measure = 0.500000 +DEAL:3D::ReferenceCell: Tet +DEAL:3D:: computed face measure = 0.500000 +DEAL:3D:: self-reported face measure = 0.500000 +DEAL:3D::ReferenceCell: Tet +DEAL:3D:: computed face measure = 0.500000 +DEAL:3D:: self-reported face measure = 0.500000 +DEAL:3D::ReferenceCell: Tet +DEAL:3D:: computed face measure = 0.866025 +DEAL:3D:: self-reported face measure = 0.866025 +DEAL:3D::ReferenceCell: Pyramid +DEAL:3D:: computed face measure = 4.00000 +DEAL:3D:: self-reported face measure = 4.00000 +DEAL:3D::ReferenceCell: Pyramid +DEAL:3D:: computed face measure = 1.41421 +DEAL:3D:: self-reported face measure = 1.41421 +DEAL:3D::ReferenceCell: Pyramid +DEAL:3D:: computed face measure = 1.41421 +DEAL:3D:: self-reported face measure = 1.41421 +DEAL:3D::ReferenceCell: Pyramid +DEAL:3D:: computed face measure = 1.41421 +DEAL:3D:: self-reported face measure = 1.41421 +DEAL:3D::ReferenceCell: Pyramid +DEAL:3D:: computed face measure = 1.41421 +DEAL:3D:: self-reported face measure = 1.41421 +DEAL:3D::ReferenceCell: Wedge +DEAL:3D:: computed face measure = 0.500000 +DEAL:3D:: self-reported face measure = 0.500000 +DEAL:3D::ReferenceCell: Wedge +DEAL:3D:: computed face measure = 0.500000 +DEAL:3D:: self-reported face measure = 0.500000 +DEAL:3D::ReferenceCell: Wedge +DEAL:3D:: computed face measure = 1.00000 +DEAL:3D:: self-reported face measure = 1.00000 +DEAL:3D::ReferenceCell: Wedge +DEAL:3D:: computed face measure = 1.41421 +DEAL:3D:: self-reported face measure = 1.41421 +DEAL:3D::ReferenceCell: Wedge +DEAL:3D:: computed face measure = 1.00000 +DEAL:3D:: self-reported face measure = 1.00000 +DEAL:3D::ReferenceCell: Hex +DEAL:3D:: computed face measure = 1.00000 +DEAL:3D:: self-reported face measure = 1.00000 +DEAL:3D::ReferenceCell: Hex +DEAL:3D:: computed face measure = 1.00000 +DEAL:3D:: self-reported face measure = 1.00000 +DEAL:3D::ReferenceCell: Hex +DEAL:3D:: computed face measure = 1.00000 +DEAL:3D:: self-reported face measure = 1.00000 +DEAL:3D::ReferenceCell: Hex +DEAL:3D:: computed face measure = 1.00000 +DEAL:3D:: self-reported face measure = 1.00000 +DEAL:3D::ReferenceCell: Hex +DEAL:3D:: computed face measure = 1.00000 +DEAL:3D:: self-reported face measure = 1.00000 +DEAL:3D::ReferenceCell: Hex +DEAL:3D:: computed face measure = 1.00000 +DEAL:3D:: self-reported face measure = 1.00000