* 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
* 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>
{
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>::
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)
{
// 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,
}
else
{
+ // we are closer to one of the vertices than to the
+ // center of the face
switch (min_index)
{
case 0:
break;
}
}
- Tensor<1,spacedim> normal = cross_product_3d(t1,t2);
+
+ const Tensor<1,spacedim> normal = cross_product_3d(t1,t2);
return normal/normal.norm();
}
projects / dealii.git / commitdiff