From: Rene Gassmoeller Date: Mon, 6 Nov 2017 18:08:54 +0000 (-0700) Subject: Test new algorithm X-Git-Tag: v9.0.0-rc1~811^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=15c60b32da11be5ba75478560eca110e78a98c49;p=dealii.git Test new algorithm --- diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index 69e70fd527..8653fc26c1 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -156,49 +156,48 @@ normal_vector (const Triangulation<3, 3>::face_iterator &face, const Point<3> &p) const { const int spacedim=3; + Tensor<1,spacedim> t1,t2,normal; - // value of i used to compute t1 - unsigned int i_t1 = numbers::invalid_unsigned_int; - Tensor<1,spacedim> t1,t2; - Tensor<1,spacedim> normal; - double normal_direction = numbers::signaling_nan(); - - // Look for a combination of tangent vectors that - // are of approximately equal length and not linearly dependent - for (unsigned int i=0,j=1 ; - i::vertices_per_face && j::vertices_per_face; - ++j) + // the counter-clockwise sequence of face vertices (see + // the GeometryInfo class for more documentation). + const double face_vertices[4] = {0,1,3,2}; + + // Look for the vertex with the largest distance to the given point + unsigned int max_index=0; + double max_distance = (p-face->vertex(face_vertices[0])).norm_square(); + + for (unsigned int i=1; i::vertices_per_face; ++i) { - // if p is too close to vertex i try again with different i and j - if ((p - face->vertex(i)).norm_square() < - 1e4 * std::numeric_limits::epsilon() * (p - face->vertex(j)).norm_square()) + const double distance = (p-face->vertex(face_vertices[i])).norm_square(); + if (distance > max_distance) { - ++i; - continue; + max_index = i; + max_distance = distance; } + } - // if p is too close to vertex j try again with different j - if ((p - face->vertex(j)).norm_square() < - 1e4 * std::numeric_limits::epsilon() * (p - face->vertex(i)).norm_square()) - continue; + // Compute tangent for max_index vertex. + t1 = get_tangent_vector(p, face->vertex(face_vertices[max_index])); - // i might not have changed, and get_tangent_vector is potentially expensive. - // Only compute a new t1 if necessary - if (i != i_t1) - { - i_t1 = i; - t1 = get_tangent_vector(p, face->vertex(i)); - } + // If p is sufficiently far away from next counter-clockwise vertex, compute normal. + const unsigned int next_index = (max_index + 1) % 4; - t2 = get_tangent_vector(p, face->vertex(j)); + if ((p-face->vertex(face_vertices[next_index])).norm_square() > + 1e4 * std::numeric_limits::epsilon() * max_distance) + { + t2 = get_tangent_vector(p, face->vertex(face_vertices[next_index])); normal = cross_product_3d(t1,t2); + } - // if t1 and t2 are (nearly) linearly dependent try again with different j / t2 - if (normal.norm_square() < 1e4 * std::numeric_limits::epsilon() * - t1.norm_square() * t2.norm_square()) - continue; - - break; + // If p is too close to next_index, or the tangents are linearly dependent + // fall back to previous counter clockwise vertex. + // In this case the normal direction is flipped. + if (normal.norm_square() < 1e4 * std::numeric_limits::epsilon() * + t1.norm_square() * t2.norm_square()) + { + const unsigned int previous_index = (max_index - 1) % 4; + t2 = get_tangent_vector(p, face->vertex(face_vertices[previous_index])); + normal = cross_product_3d(t2,t1); } Assert(normal.norm_square() >= 1e4 * std::numeric_limits::epsilon() * @@ -207,17 +206,7 @@ normal_vector (const Triangulation<3, 3>::face_iterator &face, "of vertices to compute a normal on this face. Check for distorted " "faces in your triangulation.")); - // Make sure all found normal vectors on this face point in the same direction - // as the 'reference' normal vector created at the center position. - const Point center = face->center(); - const Tensor<1,spacedim> reference_t1 = get_tangent_vector(center, face->vertex(0)); - const Tensor<1,spacedim> reference_t2 = get_tangent_vector(center, face->vertex(1)); - const Tensor<1,spacedim> reference_normal = cross_product_3d(reference_t1,reference_t2); - - if (reference_normal * normal < 0.0) - normal *= -1; - - return normal/normal.norm(); + return normal / normal.norm(); }