Make Mapping::get_center() work with simplices.

include/deal.II/fe/mapping.h

@@ -348,7 +348,7 @@ class Mapping : public Subscriptor
     const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
 
   /**
-   * Return the mapped center of a cell.
+   * Return one of two possible mapped centers of a cell.
    *
    * If you are using a (bi-,tri-)linear mapping that preserves vertex
    * locations, this function simply returns the value also produced by
@@ -358,21 +358,21 @@ class Mapping : public Subscriptor
    * polynomials, for which the center may not coincide with the average of
    * the vertex locations.
    *
-   * By default, this function returns the push forward of the center of the
+   * By default, this function returns the push forward of the barycenter of the
    * reference cell. If the parameter
-   * @p map_center_of_reference_cell is set to false, than the return value
-   * will be the average of the vertex locations, as returned by the
+   * @p map_barycenter_of_reference_cell is set to false, than the returned
+   * value will be the average of the vertex locations, as returned by the
    * get_vertices() method.
    *
    * @param[in] cell The cell for which you want to compute the center
-   * @param[in] map_center_of_reference_cell A flag that switches the algorithm
-   * for the computation of the cell center from
+   * @param[in] map_barycenter_of_reference_cell A flag that switches the
+   * algorithm for the computation of the cell center from
    * transform_unit_to_real_cell() applied to the center of the reference cell
    * to computing the vertex averages.
    */
   virtual Point<spacedim>
   get_center(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-             const bool map_center_of_reference_cell = true) const;
+             const bool map_barycenter_of_reference_cell = true) const;
 
   /**
    * Return the bounding box of a mapped cell.

source/fe/mapping.cc

@@ -54,22 +54,20 @@ template <int dim, int spacedim>
 Point<spacedim>
 Mapping<dim, spacedim>::get_center(
   const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-  const bool map_center_of_reference_cell) const
+  const bool map_barycenter_of_reference_cell) const
 {
-  if (map_center_of_reference_cell)
+  if (map_barycenter_of_reference_cell)
     {
-      Point<dim> reference_center;
-      for (unsigned int d = 0; d < dim; ++d)
-        reference_center[d] = .5;
-      return transform_unit_to_real_cell(cell, reference_center);
+      return transform_unit_to_real_cell(
+        cell, cell->reference_cell().template barycenter<dim>());
     }
   else
     {
       const auto      vertices = get_vertices(cell);
       Point<spacedim> center;
       for (const auto &v : vertices)
         center += v;
-      return center / GeometryInfo<dim>::vertices_per_cell;
+      return center / cell->n_vertices();
     }
 }
 

tests/simplex/mapping_01.cc

@@ -0,0 +1,105 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 - 2023 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 mapping centers with simplices
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_fe_field.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+class Position : public Function<dim>
+{
+public:
+  Position()
+    : Function<dim>(dim)
+  {}
+
+  double
+  value(const Point<dim> &point, const unsigned int component) const
+  {
+    if (component == 0)
+      return point[component] + std::sin(point[dim - 1] * numbers::PI);
+    return point[component];
+  }
+};
+
+void
+test(const unsigned int mapping_degree)
+{
+  const int dim = 3;
+
+  Triangulation<dim> tria;
+  GridGenerator::subdivided_hyper_cube_with_simplices(tria, 2);
+
+  FE_SimplexP<dim> fe(mapping_degree);
+  FESystem<dim>    position_fe(fe, dim);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  DoFHandler<dim> position_dof_handler(tria);
+  position_dof_handler.distribute_dofs(position_fe);
+
+  Vector<double> position_vector(position_dof_handler.n_dofs());
+  MappingFE<dim> mapping_interpolation(FE_SimplexP<dim>(1));
+  VectorTools::interpolate(mapping_interpolation,
+                           position_dof_handler,
+                           Position<dim>(),
+                           position_vector);
+
+  MappingFEField<dim> mapping(position_dof_handler, position_vector);
+
+  for (const auto &cell : tria.active_cell_iterators())
+    {
+      deallog << "cell = " << cell << std::endl;
+      deallog << "mapped barycenter       = " << mapping.get_center(cell, true)
+              << std::endl;
+      Point<dim> manually_mapped_center;
+      for (unsigned int d = 0; d < dim; ++d)
+        manually_mapped_center[d] =
+          Position<dim>().value(cell->barycenter(), d);
+      deallog << "exact mapped barycenter = " << manually_mapped_center
+              << std::endl;
+      deallog << "mapped vertex mean      = " << mapping.get_center(cell, false)
+              << std::endl;
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  test(1); // linear mapping
+
+  deallog << std::endl << "quadratic mapping" << std::endl << std::endl;
+
+  test(2); // quadratic mapping
+}