Mapping::get_center()

diff --git a/doc/news/changes/minor/20190411LucaHeltai b/doc/news/changes/minor/20190411LucaHeltai
new file mode 100644 (file)
index 0000000..ad5ad6c
--- /dev/null
+++ b/doc/news/changes/minor/20190411LucaHeltai
@@ -0,0 +1,4 @@
+New: The method Mapping::get_center() allows one to retrieve the cell center of a cell when the 
+mapping object does not preserve vertex locations. 
+<br>
+(Luca Heltai, 2019/04/11)

diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h
index 21bec9ebc2a35b3d60047ae320ae6b7dda478b2c..c1c60aac8ba789f1db81d4680a92999e856f6150 100644 (file)
--- a/include/deal.II/fe/mapping.h
+++ b/include/deal.II/fe/mapping.h
@@ -336,6 +336,32 @@ public:
   get_vertices(
     const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
 
+  /**
+   * Return the mapped center 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
+   * `cell->center()`. However, there are also mappings that add displacements
+   * or choose completely different locations, e.g., MappingQEulerian,
+   * MappingQ1Eulerian, or MappingFEField.
+   *
+   * By default, this function returns the average of the vertex locations
+   * returned by the the get_vertices() method. If the parameter
+   * @p map_center_of_reference_cell is set to true, than a more expensive
+   * algorithm is used, that returns the mapped center of the reference cell,
+   * obtained by the transform_unit_to_real_cell() 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 vertex averages to
+   * transform_unit_to_real_cell() applied to the center of the reference cell
+   *
+   * @author Luca Heltai, 2019.
+   */
+  virtual Point<spacedim>
+  get_center(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+             const bool map_center_of_reference_cell = false) const;
+
   /**
    * Return whether the mapping preserves vertex locations. In other words,
    * this function returns whether the mapped location of the reference cell

diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc
index 60d57958000e9a02ad8d0abfa33e4ca2d759a2a4..97a7a38f236469fc4a9a09b60da94245125c6005 100644 (file)
--- a/source/fe/mapping.cc
+++ b/source/fe/mapping.cc
@@ -36,6 +36,31 @@ Mapping<dim, spacedim>::get_vertices(
 
 
 
+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
+{
+  if (map_center_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);
+    }
+  else
+    {
+      const auto      vertices = get_vertices(cell);
+      Point<spacedim> center;
+      for (const auto &v : vertices)
+        center += v;
+      return center / GeometryInfo<dim>::vertices_per_cell;
+    }
+}
+
+
+
 template <int dim, int spacedim>
 Point<dim - 1>
 Mapping<dim, spacedim>::project_real_point_to_unit_point_on_face(

diff --git a/tests/mappings/mapping_q_eulerian_12.cc b/tests/mappings/mapping_q_eulerian_12.cc
new file mode 100644 (file)
index 0000000..e862658
--- /dev/null
+++ b/tests/mappings/mapping_q_eulerian_12.cc
@@ -0,0 +1,96 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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 that MappingQEulerian knows how to compute centers of cells
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q_eulerian.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+
+template <int dim>
+class Displacement : public Function<dim>
+{
+public:
+  Displacement()
+    : Function<dim>(dim)
+  {}
+
+  double
+  value(const Point<dim> &p, const unsigned int component) const
+  {
+    return p[component] * p[0] / 2;
+  }
+};
+
+
+template <int dim>
+void
+test()
+{
+  deallog << "dim=" << dim << std::endl;
+
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation, -1, 1);
+
+  FESystem<dim>   fe(FE_Q<dim>(2), dim);
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(fe);
+
+  Vector<double> displacements(dof_handler.n_dofs());
+
+  VectorTools::interpolate(dof_handler, Displacement<dim>(), displacements);
+
+  MappingQEulerian<dim> euler(2, dof_handler, displacements);
+
+  // now the actual test
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      deallog << "Center: [" << cell->center() << "], with mapping ["
+              << euler.get_center(cell) << "]" << std::endl;
+      deallog << "Center: [" << cell->center() << "], with mapping + flag ["
+              << euler.get_center(cell, true) << "]" << std::endl;
+    }
+}
+
+
+
+int
+main()
+{
+  initlog(1);
+
+  test<1>();
+  test<2>();
+  test<3>();
+}

diff --git a/tests/mappings/mapping_q_eulerian_12.output b/tests/mappings/mapping_q_eulerian_12.output
new file mode 100644 (file)
index 0000000..eb04a22
--- /dev/null
+++ b/tests/mappings/mapping_q_eulerian_12.output
@@ -0,0 +1,10 @@
+
+DEAL::dim=1
+DEAL::Center: [0.00000], with mapping [0.500000]
+DEAL::Center: [0.00000], with mapping + flag [0.00000]
+DEAL::dim=2
+DEAL::Center: [0.00000 0.00000], with mapping [0.500000 0.00000]
+DEAL::Center: [0.00000 0.00000], with mapping + flag [0.00000 0.00000]
+DEAL::dim=3
+DEAL::Center: [0.00000 0.00000 0.00000], with mapping [0.500000 0.00000 0.00000]
+DEAL::Center: [0.00000 0.00000 0.00000], with mapping + flag [0.00000 0.00000 0.00000]