From a2119c4cbf9b2be19b4dde574d5e08b1247887ca Mon Sep 17 00:00:00 2001
From: Rene Gassmoeller <rene.gassmoeller@mailbox.org>
Date: Mon, 24 Feb 2020 15:06:41 -0500
Subject: [PATCH] Add PolarManifold::normal_vector

---
 include/deal.II/grid/manifold_lib.h |  8 +++++
 source/grid/manifold_lib.cc         | 52 +++++++++++++++++++++++++++++
 2 files changed, 60 insertions(+)

diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h
index 0a4a3ec307..1211f9c912 100644
--- a/include/deal.II/grid/manifold_lib.h
+++ b/include/deal.II/grid/manifold_lib.h
@@ -110,6 +110,14 @@ public:
   virtual DerivativeForm<1, spacedim, spacedim>
   push_forward_gradient(const Point<spacedim> &chart_point) const override;
 
+  /**
+   * Return the (normalized) normal vector at the point @p p.
+   */
+  virtual Tensor<1, spacedim>
+  normal_vector(
+    const typename Triangulation<dim, spacedim>::face_iterator &face,
+    const Point<spacedim> &p) const override;
+
   /**
    * The center of the spherical coordinate system.
    */
diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index 1b8657e790..c9bf5ba433 100644
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -274,6 +274,58 @@ PolarManifold<dim, spacedim>::push_forward_gradient(
 
 
 
+template <int dim, int spacedim>
+Tensor<1, spacedim>
+PolarManifold<dim, spacedim>::normal_vector(
+  const typename Triangulation<dim, spacedim>::face_iterator &face,
+  const Point<spacedim> &                                     p) const
+{
+  // Let us first test whether we are on a "horizontal" face.
+  // In that case, the normal vector is easy to compute
+  // since it is proportional to the vector from the center to the
+  // point 'p'.
+  //
+  // We test whether that is the case by checking that the vertices
+  // all have roughly the same distance from the center: If the
+  // maximum deviation for the distances from the vertices to the
+  // center is less than 1.e-5 of the distance between vertices (as
+  // measured by the minimum distance from any of the other vertices
+  // to the first vertex), then we call this a horizontal face.
+  constexpr unsigned int n_vertices = GeometryInfo<spacedim>::vertices_per_face;
+  std::array<double, n_vertices>     distances_to_center;
+  std::array<double, n_vertices - 1> distances_to_first_vertex;
+  distances_to_center[0] = (face->vertex(0) - center).norm_square();
+  for (unsigned int i = 1; i < n_vertices; ++i)
+    {
+      distances_to_center[i] = (face->vertex(i) - center).norm_square();
+      distances_to_first_vertex[i - 1] =
+        (face->vertex(i) - face->vertex(0)).norm_square();
+    }
+  const auto minmax_distance =
+    std::minmax_element(distances_to_center.begin(), distances_to_center.end());
+  const auto min_distance_to_first_vertex =
+    std::min_element(distances_to_first_vertex.begin(),
+                     distances_to_first_vertex.end());
+
+  // So, if this is a horizontal face, then just compute the normal
+  // vector as the one from the center to the point 'p', adequately
+  // scaled.
+  if (*minmax_distance.second - *minmax_distance.first <
+      1.e-10 * *min_distance_to_first_vertex)
+    {
+      const Tensor<1, spacedim> unnormalized_spherical_normal = p - center;
+      const Tensor<1, spacedim> normalized_spherical_normal =
+        unnormalized_spherical_normal / unnormalized_spherical_normal.norm();
+      return normalized_spherical_normal;
+    }
+
+  // If it is not a horizontal face, just use the machinery of the
+  // base class.
+  return Manifold<dim, spacedim>::normal_vector(face, p);
+}
+
+
+
 // ============================================================
 // SphericalManifold
 // ============================================================
-- 
2.39.5