From a2119c4cbf9b2be19b4dde574d5e08b1247887ca Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller 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 &chart_point) const override; + /** + * Return the (normalized) normal vector at the point @p p. + */ + virtual Tensor<1, spacedim> + normal_vector( + const typename Triangulation::face_iterator &face, + const Point &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::push_forward_gradient( +template +Tensor<1, spacedim> +PolarManifold::normal_vector( + const typename Triangulation::face_iterator &face, + const Point & 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::vertices_per_face; + std::array distances_to_center; + std::array 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::normal_vector(face, p); +} + + + // ============================================================ // SphericalManifold // ============================================================ -- 2.39.5