]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Update the documentation of Manifold::normal_vector(). 3198/head
authorWolfgang Bangerth <bangerth@colostate.edu>
Wed, 5 Oct 2016 00:36:56 +0000 (18:36 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Wed, 5 Oct 2016 00:36:56 +0000 (18:36 -0600)
include/deal.II/grid/manifold.h
source/grid/manifold.cc

index 927ce04cd10437f0df6fd64e0cf54d26f1e08157..bee19b2757444dfc6f30eb9eba4f4c33eb479f0f 100644 (file)
@@ -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>
index b0d90b63fa09a27f0f344fa4ae29d8ac9674b10f..f36b8ca3eaf2d25db7327a793e13383cbc824723 100644 (file)
@@ -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();
 }
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.