Merge pull request #13678 from bangerth/ref-cell

Add ReferenceCell::volume() and ::barycenter().

doc/news/changes/minor/20220504Bangerth

@@ -0,0 +1,5 @@
+New: There are now functions ReferenceCell::volume() and
+ReferenceCells::barycenter() that compute what their names suggest for
+the various kinds of reference cells deal.II supports.
+<br>
+(Wolfgang Bangerth, 2022/05/04)

include/deal.II/grid/reference_cell.h

@@ -457,6 +457,37 @@ class ReferenceCell
    * @{
    */
 
+  /**
+   * Return the $d$-dimensional volume of the reference cell that corresponds
+   * to the current object, where $d$ is the dimension of the space it lives
+   * in. For example, since the quadrilateral reference cell is $[0,1]^2$,
+   * its volume is one, whereas the volume of the reference triangle is
+   * 0.5 because it occupies the area $\{0 \le x,y \le 1, x+y\le 1\}$.
+   *
+   * For ReferenceCells::Vertex, the reference cell is a zero-dimensional
+   * point in a zero-dimensional space. As a consequence, one cannot
+   * meaningfully define a volume for it. The function returns one for
+   * this case, because this makes it possible to define useful quadrature
+   * rules based on the center of a reference cell and its volume.
+   */
+  double
+  volume() const;
+
+  /**
+   * Return the barycenter of the reference cell that corresponds
+   * to the current object. The function is not called `center()` because
+   * one can define the center of an object in a number of different ways
+   * whereas the barycenter of a reference cell $K$ is unambiguously defined as
+   * @f[
+   *   \mathbf x_K = \frac{1}{V} \int_K \mathbf x \; dx
+   * @f]
+   * where $V$ is the volume of the reference cell (see also the volume()
+   * function).
+   */
+  template <int dim>
+  Point<dim>
+  barycenter() const;
+
   /**
    * Return true if the given point is inside the reference cell of the present
    * space dimension up to some tolerance. This function accepts an additional
@@ -1922,6 +1953,61 @@ ReferenceCell::d_linear_shape_function_gradient(const Point<dim> & xi,
 
 
 
+inline double
+ReferenceCell::volume() const
+{
+  if (*this == ReferenceCells::Vertex)
+    return 0;
+  else if (*this == ReferenceCells::Line)
+    return 1;
+  else if (*this == ReferenceCells::Triangle)
+    return 1. / 2.;
+  else if (*this == ReferenceCells::Quadrilateral)
+    return 1;
+  else if (*this == ReferenceCells::Tetrahedron)
+    return 1. / 6.;
+  else if (*this == ReferenceCells::Wedge)
+    return 1. / 2.;
+  else if (*this == ReferenceCells::Pyramid)
+    return 4. / 3.;
+  else if (*this == ReferenceCells::Hexahedron)
+    return 1;
+
+  Assert(false, ExcNotImplemented());
+  return 0.0;
+}
+
+
+
+template <int dim>
+inline Point<dim>
+ReferenceCell::barycenter() const
+{
+  AssertDimension(dim, get_dimension());
+
+  if (*this == ReferenceCells::Vertex)
+    return Point<dim>();
+  else if (*this == ReferenceCells::Line)
+    return Point<dim>(1. / 2.);
+  else if (*this == ReferenceCells::Triangle)
+    return Point<dim>(1. / 3., 1. / 3.);
+  else if (*this == ReferenceCells::Quadrilateral)
+    return Point<dim>(1. / 2., 1. / 2.);
+  else if (*this == ReferenceCells::Tetrahedron)
+    return Point<dim>(1. / 4., 1. / 4., 1. / 4.);
+  else if (*this == ReferenceCells::Wedge)
+    return Point<dim>(1. / 3, 1. / 3, 1. / 2.);
+  else if (*this == ReferenceCells::Pyramid)
+    return Point<dim>(0, 0, 1. / 4.);
+  else if (*this == ReferenceCells::Hexahedron)
+    return Point<dim>(1. / 2., 1. / 2., 1. / 2.);
+
+  Assert(false, ExcNotImplemented());
+  return Point<dim>();
+}
+
+
+
 template <int dim>
 inline bool
 ReferenceCell::contains_point(const Point<dim> &p, const double tolerance) const

tests/grid/reference_cell_type_02.cc

@@ -0,0 +1,99 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Test ReferenceCell::volume()
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test(const ReferenceCell &reference_cell)
+{
+  Triangulation<dim> triangulation;
+  GridGenerator::reference_cell(triangulation, reference_cell);
+
+  std::unique_ptr<Quadrature<dim>> q;
+  if (reference_cell == ReferenceCells::Line)
+    q = std::make_unique<QGauss<dim>>(2);
+  else if (reference_cell == ReferenceCells::Triangle)
+    q = std::make_unique<QGaussSimplex<dim>>(2);
+  else if (reference_cell == ReferenceCells::Quadrilateral)
+    q = std::make_unique<QGauss<dim>>(2);
+  else if (reference_cell == ReferenceCells::Tetrahedron)
+    q = std::make_unique<QGaussSimplex<dim>>(2);
+  else if (reference_cell == ReferenceCells::Wedge)
+    q = std::make_unique<QGaussWedge<dim>>(2);
+  else if (reference_cell == ReferenceCells::Pyramid)
+    q = std::make_unique<QGaussPyramid<dim>>(2);
+  else if (reference_cell == ReferenceCells::Hexahedron)
+    q = std::make_unique<QGauss<dim>>(2);
+
+  FE_Nothing<dim> fe(reference_cell);
+
+  // Set up the objects to compute an integral on the reference cell
+  FEValues<dim> fe_values(fe, *q, update_JxW_values);
+  fe_values.reinit(triangulation.begin_active());
+
+  double volume = 0;
+  for (unsigned int i = 0; i < q->size(); ++i)
+    volume += fe_values.JxW(i);
+
+  deallog << "ReferenceCell: " << reference_cell.to_string() << std::endl;
+  deallog << "  computed volume = " << volume << std::endl;
+  deallog << "  self-reported volume = " << reference_cell.volume()
+          << std::endl;
+
+  Assert(std::fabs(volume - reference_cell.volume()) < 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();
+  }
+}

tests/grid/reference_cell_type_02.output

@@ -0,0 +1,22 @@
+
+DEAL:1D::ReferenceCell: Line
+DEAL:1D::  computed volume = 1.00000
+DEAL:1D::  self-reported volume = 1.00000
+DEAL:2D::ReferenceCell: Quad
+DEAL:2D::  computed volume = 1.00000
+DEAL:2D::  self-reported volume = 1.00000
+DEAL:2D::ReferenceCell: Tri
+DEAL:2D::  computed volume = 0.500000
+DEAL:2D::  self-reported volume = 0.500000
+DEAL:3D::ReferenceCell: Tet
+DEAL:3D::  computed volume = 0.166667
+DEAL:3D::  self-reported volume = 0.166667
+DEAL:3D::ReferenceCell: Pyramid
+DEAL:3D::  computed volume = 1.33333
+DEAL:3D::  self-reported volume = 1.33333
+DEAL:3D::ReferenceCell: Wedge
+DEAL:3D::  computed volume = 0.500000
+DEAL:3D::  self-reported volume = 0.500000
+DEAL:3D::ReferenceCell: Hex
+DEAL:3D::  computed volume = 1.00000
+DEAL:3D::  self-reported volume = 1.00000

tests/grid/reference_cell_type_03.cc

@@ -0,0 +1,104 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Test ReferenceCell::barycenter()
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test(const ReferenceCell &reference_cell)
+{
+  Triangulation<dim> triangulation;
+  GridGenerator::reference_cell(triangulation, reference_cell);
+
+  std::unique_ptr<Quadrature<dim>> q;
+  if (reference_cell == ReferenceCells::Line)
+    q = std::make_unique<QGauss<dim>>(2);
+  else if (reference_cell == ReferenceCells::Triangle)
+    q = std::make_unique<QGaussSimplex<dim>>(2);
+  else if (reference_cell == ReferenceCells::Quadrilateral)
+    q = std::make_unique<QGauss<dim>>(2);
+  else if (reference_cell == ReferenceCells::Tetrahedron)
+    q = std::make_unique<QGaussSimplex<dim>>(2);
+  else if (reference_cell == ReferenceCells::Wedge)
+    q = std::make_unique<QGaussWedge<dim>>(2);
+  else if (reference_cell == ReferenceCells::Pyramid)
+    q = std::make_unique<QGaussPyramid<dim>>(2);
+  else if (reference_cell == ReferenceCells::Hexahedron)
+    q = std::make_unique<QGauss<dim>>(2);
+
+  FE_Nothing<dim> fe(reference_cell);
+
+  // Set up the objects to compute an integral on the reference cell
+  FEValues<dim> fe_values(fe, *q, update_JxW_values | update_quadrature_points);
+  fe_values.reinit(triangulation.begin_active());
+
+  double     volume = 0;
+  Point<dim> barycenter;
+  for (unsigned int i = 0; i < q->size(); ++i)
+    {
+      volume += fe_values.JxW(i);
+      barycenter += fe_values.quadrature_point(i) * fe_values.JxW(i);
+    }
+  barycenter /= volume;
+
+  deallog << "ReferenceCell: " << reference_cell.to_string() << std::endl;
+  deallog << "  computed barycenter = " << barycenter << std::endl;
+  deallog << "  self-reported barycenter = " << reference_cell.barycenter<dim>()
+          << std::endl;
+
+  Assert((barycenter - reference_cell.barycenter<dim>()).norm() <= 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();
+  }
+}

tests/grid/reference_cell_type_03.output

@@ -0,0 +1,22 @@
+
+DEAL:1D::ReferenceCell: Line
+DEAL:1D::  computed barycenter = 0.500000
+DEAL:1D::  self-reported barycenter = 0.500000
+DEAL:2D::ReferenceCell: Quad
+DEAL:2D::  computed barycenter = 0.500000 0.500000
+DEAL:2D::  self-reported barycenter = 0.500000 0.500000
+DEAL:2D::ReferenceCell: Tri
+DEAL:2D::  computed barycenter = 0.333333 0.333333
+DEAL:2D::  self-reported barycenter = 0.333333 0.333333
+DEAL:3D::ReferenceCell: Tet
+DEAL:3D::  computed barycenter = 0.250000 0.250000 0.250000
+DEAL:3D::  self-reported barycenter = 0.250000 0.250000 0.250000
+DEAL:3D::ReferenceCell: Pyramid
+DEAL:3D::  computed barycenter = 0 0 0.250000
+DEAL:3D::  self-reported barycenter = 0.00000 0.00000 0.250000
+DEAL:3D::ReferenceCell: Wedge
+DEAL:3D::  computed barycenter = 0.333333 0.333333 0.500000
+DEAL:3D::  self-reported barycenter = 0.333333 0.333333 0.500000
+DEAL:3D::ReferenceCell: Hex
+DEAL:3D::  computed barycenter = 0.500000 0.500000 0.500000
+DEAL:3D::  self-reported barycenter = 0.500000 0.500000 0.500000