From 100405c5bf01be3000ac1e50db4ff97c1a1054c9 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 24 Feb 2020 10:52:05 -0700 Subject: [PATCH] Update documentation of SphericalManifold::normal_vector. --- source/grid/manifold_lib.cc | 95 +++++++++++++++++++++++-------------- 1 file changed, 60 insertions(+), 35 deletions(-) diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index 9b190eab23..1b8657e790 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -419,35 +419,49 @@ SphericalManifold::normal_vector( const typename Triangulation::face_iterator &face, const Point & p) const { - // if the maximum deviation for the distance from the vertices to the center - // is less than 1.e-5 of the minimum distance to the first vertex, assume we - // can simply return p-center. otherwise, we compute the normal using - // get_normal_vector + // Let us first test whether we are on a "horizontal" face + // (tangential to the sphere). In this 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(); + std::array sqr_distances_to_center; + std::array sqr_distances_to_first_vertex; + sqr_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] = + sqr_distances_to_center[i] = (face->vertex(i) - center).norm_square(); + sqr_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()); - - if (*minmax_distance.second - *minmax_distance.first < - 1.e-10 * *min_distance_to_first_vertex) + const auto minmax_sqr_distance = + std::minmax_element(sqr_distances_to_center.begin(), + sqr_distances_to_center.end()); + const auto min_sqr_distance_to_first_vertex = + std::min_element(sqr_distances_to_first_vertex.begin(), + sqr_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_sqr_distance.second - *minmax_sqr_distance.first < + 1.e-10 * *min_sqr_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; } - return Manifold::normal_vector(face, p); + else + // If it is not a horizontal face, just use the machinery of the + // base class. + return Manifold::normal_vector(face, p); } @@ -481,28 +495,39 @@ SphericalManifold::get_normals_at_vertices( typename Manifold::FaceVertexNormals &face_vertex_normals) const { - // if the maximum deviation for the distance from the vertices to the center - // is less than 1.e-5 of the minimum distance to the first vertex, assume we - // can simply return vertex-center. otherwise, we compute the normal using - // get_normal_vector + // Let us first test whether we are on a "horizontal" face + // (tangential to the sphere). In this 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(); + std::array sqr_distances_to_center; + std::array sqr_distances_to_first_vertex; + sqr_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] = + sqr_distances_to_center[i] = (face->vertex(i) - center).norm_square(); + sqr_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()); - - if (*minmax_distance.second - *minmax_distance.first < - 1.e-10 * *min_distance_to_first_vertex) + const auto minmax_sqr_distance = + std::minmax_element(sqr_distances_to_center.begin(), + sqr_distances_to_center.end()); + const auto min_sqr_distance_to_first_vertex = + std::min_element(sqr_distances_to_first_vertex.begin(), + sqr_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_sqr_distance.second - *minmax_sqr_distance.first < + 1.e-10 * *min_sqr_distance_to_first_vertex) { for (unsigned int vertex = 0; vertex < n_vertices; ++vertex) face_vertex_normals[vertex] = face->vertex(vertex) - center; -- 2.39.5