Get mapped vertices of face

diff --git a/doc/news/changes/minor/20230125Heinz b/doc/news/changes/minor/20230125Heinz
new file mode 100644 (file)
index 0000000..c4bed5b
--- /dev/null
+++ b/doc/news/changes/minor/20230125Heinz
@@ -0,0 +1,3 @@
+Improved: Mapping can now return vertices on a face of a cell.
+<br>
+(Johannes Heinz, 2023/01/25)

diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h
index 9c2568abc6704378b0c1c213c73d911d4643d47f..99c3fc210919d89ef518305c0a7d5f675eb33d27 100644 (file)
--- a/include/deal.II/fe/mapping.h
+++ b/include/deal.II/fe/mapping.h
@@ -347,6 +347,19 @@ public:
   get_vertices(
     const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
 
+  /**
+   * Return the mapped vertices of a face.
+   *
+   * Same as above but working on a given face of a cell.
+   *
+   * @param[in] cell The cell containing the face.
+   * @param[in] face_no The number of the face within the cell.
+   */
+  boost::container::small_vector<Point<spacedim>,
+                                 GeometryInfo<dim>::vertices_per_face>
+  get_vertices(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+               const unsigned int face_no) const;
+
   /**
    * Return the mapped center of a cell.
    *

diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc
index dc87105a038de6367cef56fbdbd7923a68cf0d72..f4f7fdfec4a530854c069982592629591ed7c45b 100644 (file)
--- a/source/fe/mapping.cc
+++ b/source/fe/mapping.cc
@@ -52,6 +52,34 @@ Mapping<dim, spacedim>::get_vertices(
 
 
 
+template <int dim, int spacedim>
+boost::container::small_vector<Point<spacedim>,
+                               GeometryInfo<dim>::vertices_per_face>
+Mapping<dim, spacedim>::get_vertices(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const unsigned int                                          face_no) const
+{
+  boost::container::small_vector<Point<spacedim>,
+                                 GeometryInfo<dim>::vertices_per_face>
+    face_vertices;
+
+  const auto &cell_vertices    = get_vertices(cell);
+  const auto &reference_cell   = cell->reference_cell();
+  const auto  face_orientation = cell->combined_face_orientation(face_no);
+
+  for (const unsigned int v :
+       reference_cell.face_reference_cell(face_no).vertex_indices())
+    {
+      face_vertices.push_back(
+        cell_vertices[reference_cell.face_to_cell_vertices(
+          face_no, v, face_orientation)]);
+    }
+
+  return face_vertices;
+}
+
+
+
 template <int dim, int spacedim>
 Point<spacedim>
 Mapping<dim, spacedim>::get_center(

diff --git a/tests/mappings/mapping_get_vertices_on_face.cc b/tests/mappings/mapping_get_vertices_on_face.cc
new file mode 100644 (file)
index 0000000..53ef12b
--- /dev/null
+++ b/tests/mappings/mapping_get_vertices_on_face.cc
@@ -0,0 +1,84 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2001 - 2022 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Checks if vertices on face are given in correct oder
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+std::vector<Point<dim>>
+get_vertices(const typename Triangulation<dim>::cell_iterator &cell,
+             const typename Triangulation<dim>::face_iterator &face,
+             const Mapping<dim> &                              mapping)
+{
+  std::vector<Point<dim>> vertices(face->n_vertices());
+
+  for (unsigned int i = 0; i < face->n_vertices(); ++i)
+    {
+      vertices[i] = mapping.transform_unit_to_real_cell(
+        cell, mapping.transform_real_to_unit_cell(cell, face->vertex(i)));
+    }
+  return vertices;
+}
+
+template <int dim>
+void
+test_vertex_order()
+{
+  deallog << "dim=" << dim << std::endl;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, -1.0, 1.0);
+
+  MappingQ<dim> mapping(1);
+
+  for (const auto &cell : tria.active_cell_iterators())
+    for (unsigned int f = 0; f < cell->n_faces(); ++f)
+      {
+        auto const &vertices_detour =
+          get_vertices(cell, cell->face(f), mapping);
+        auto const &vertices = mapping.get_vertices(cell, f);
+
+        AssertDimension(vertices_detour.size(), vertices.size());
+
+        for (unsigned int i = 0; i < vertices.size(); ++i)
+          {
+            auto const diff = vertices[i] - vertices_detour[i];
+            AssertThrow(diff.norm() < 1e-12, ExcMessage("Vertices differ!"));
+          }
+      }
+
+  deallog << "OK" << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(12);
+
+  test_vertex_order<1>();
+  test_vertex_order<2>();
+  test_vertex_order<3>();
+
+  return 0;
+}

diff --git a/tests/mappings/mapping_get_vertices_on_face.output b/tests/mappings/mapping_get_vertices_on_face.output
new file mode 100644 (file)
index 0000000..4afcfeb
--- /dev/null
+++ b/tests/mappings/mapping_get_vertices_on_face.output
@@ -0,0 +1,7 @@
+
+DEAL::dim=1
+DEAL::OK
+DEAL::dim=2
+DEAL::OK
+DEAL::dim=3
+DEAL::OK