From bb16819dc905abcd6257d7ba6eb44ae6d4c5e99f Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 4 Oct 2016 18:36:56 -0600 Subject: [PATCH] Update the documentation of Manifold::normal_vector(). --- include/deal.II/grid/manifold.h | 38 +++++++++++++++++++++++++++---- source/grid/manifold.cc | 40 +++++++++++++++++++++------------ 2 files changed, 60 insertions(+), 18 deletions(-) diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h index 927ce04cd1..bee19b2757 100644 --- a/include/deal.II/grid/manifold.h +++ b/include/deal.II/grid/manifold.h @@ -548,6 +548,24 @@ public: * mapped surface, however, will not usually coincide with the * actual surface.) * + * This function only makes sense if dim==spacedim because + * otherwise there is no unique normal vector but in fact a + * (spacedim-dim+1)-dimensional tangent space of vectors that + * are all both normal to the face and normal to the dim-dimensional + * surface that lives in spacedim-dimensional space. For example, + * think of a two-dimensional mesh that covers a two-dimensional + * surface in three-dimensional space. In that case, each + * face (edge) is one-dimensional, and there are two linearly independent + * vectors that are both normal to the edge: one is normal to the + * edge and tangent to the surface (intuitively, that would be the + * one that points from the current cell to the neighboring one, + * if the surface was locally flat), and the other one is rooted + * in the edge but points perpendicular to the surface (which is + * also perpendicular to the edge that lives within the surface). + * Thus, because there are no obviously correct semantics for this function + * if spacedim is greater than dim, the function will simply throw + * an error in that situation. + * * The face iterator gives an indication which face this function is * supposed to compute the normal vector for. This is useful if the * boundary of the domain is composed of different nondifferential @@ -556,10 +574,22 @@ public: * mesh, with piecewise (bi-)linear components between the vertices, * but where the boundary may have a kink at the vertices itself). * - * @note The default implementation of this function computes the - * normal vector by taking the cross product between the tangent - * vectors from p to the most orthogonal and further non consecutive - * vertices of the face. + * @note In 2d, the default implementation of this function computes the + * normal vector by taking the tangent direction from p to + * the further one of the two vertices that make up an edge, and then + * rotates it outward (with respect to the coordinate system of the edge) + * by 90 degrees. In 3d, the default implementation is more + * complicated, aiming at avoiding problems with numerical round-off + * for points close to one of the vertices. If the point p is closer + * to the center of the face than to any of the vertices, the + * normal vector is computed by the cross product of the tangent + * vectors from p to either vertex zero and one of the face (if + * the closest vertex is either vertex two or three), or of the tangent + * vectors from p to vertices two and three (if the closest vertex is + * either vertex zero or one). On the other hand, if the point p + * is closer to one of the vertices than to the center of the face, + * then we take the cross product of the tangent vectors from p + * to the two vertices that are adjacent to the closest one. */ virtual Tensor<1,spacedim> diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index b0d90b63fa..f36b8ca3ea 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -139,13 +139,20 @@ normal_vector (const Triangulation<2, 2>::face_iterator &face, { const int spacedim=2; - Tensor<1,spacedim> tangent = ((p-face->vertex(0)).norm_square() > (p-face->vertex(1)).norm_square() ? - -get_tangent_vector(p, face->vertex(0)) : - get_tangent_vector(p, face->vertex(1))); - Tensor<1,spacedim> normal = cross_product_2d(tangent); + // get the tangent vector from the point 'p' in the direction of the further + // one of the two vertices that make up the face of this 2d cell + const Tensor<1,spacedim> tangent + = ((p-face->vertex(0)).norm_square() > (p-face->vertex(1)).norm_square() ? + -get_tangent_vector(p, face->vertex(0)) : + get_tangent_vector(p, face->vertex(1))); + + // then rotate it by 90 degrees + const Tensor<1,spacedim> normal = cross_product_2d(tangent); return normal/normal.norm(); } + + template<> Tensor<1,3> Manifold<3, 3>:: @@ -156,13 +163,12 @@ normal_vector (const Triangulation<3, 3>::face_iterator &face, Tensor<1,spacedim> t1,t2; // Take the difference between p and all four vertices - int min_index=0; - Tensor<1,spacedim> dp = p-face->vertex(0); - double min_distance = dp.norm_square(); + unsigned int min_index=0; + double min_distance = (p-face->vertex(0)).norm_square(); for (unsigned int i=1; i<4; ++i) { - dp = p-face->vertex(i); + const Tensor<1,spacedim> dp = p-face->vertex(i); double distance = dp.norm_square(); if (distance < min_distance) { @@ -173,11 +179,14 @@ normal_vector (const Triangulation<3, 3>::face_iterator &face, // Verify we have a valid vertex index AssertIndexRange(min_index, 4); - // Now figure out which vertices are better to compute tangent vectors - // we split the cell in 4 quadrants, and use the most orthogonal vertices - // to the closest vertex if we ar far from the center, othewise we use - // the two consecutive vertices, on the opposite side with respect to - // the face center. + // Now figure out which vertices are best to compute tangent vectors. + // We split the cell in a central diamond of points closer to the + // center than to any of the vertices, and the 4 triangles in the + // corner. The central diamond is split into its upper and lower + // half. For each of these 6 cases, the following encodes a list + // of two vertices each to which we compute the tangent vectors, + // and then take the cross product. See the documentation of this + // function for exact details. if ((p-face->center()).norm_square() < min_distance) { // we are close to the face center: pick two consecutive vertices, @@ -196,6 +205,8 @@ normal_vector (const Triangulation<3, 3>::face_iterator &face, } else { + // we are closer to one of the vertices than to the + // center of the face switch (min_index) { case 0: @@ -227,7 +238,8 @@ normal_vector (const Triangulation<3, 3>::face_iterator &face, break; } } - Tensor<1,spacedim> normal = cross_product_3d(t1,t2); + + const Tensor<1,spacedim> normal = cross_product_3d(t1,t2); return normal/normal.norm(); } -- 2.39.5