From 3279cf50e0891d3fce7fbe84a054b6132d25d06b Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Fri, 3 Nov 2017 14:42:02 -0600 Subject: [PATCH] Fix algorithm in Manifold::normal_vector --- .../changes/minor/20171103ReneGassmoeller | 5 + include/deal.II/grid/manifold.h | 12 +- source/grid/manifold.cc | 115 +++++++----------- tests/manifold/flat_manifold_08.cc | 64 ++++++++++ tests/manifold/flat_manifold_08.output | 3 + 5 files changed, 116 insertions(+), 83 deletions(-) create mode 100644 doc/news/changes/minor/20171103ReneGassmoeller create mode 100644 tests/manifold/flat_manifold_08.cc create mode 100644 tests/manifold/flat_manifold_08.output diff --git a/doc/news/changes/minor/20171103ReneGassmoeller b/doc/news/changes/minor/20171103ReneGassmoeller new file mode 100644 index 0000000000..6a94b70a2a --- /dev/null +++ b/doc/news/changes/minor/20171103ReneGassmoeller @@ -0,0 +1,5 @@ +Fixed: The algorithm in Manifold::normal_vector could +cause division-by-zero errors for special cases. This +was fixed. +
+(Rene Gassmoeller, 2017/11/03) diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h index 02dd9a4214..4926b09bb2 100644 --- a/include/deal.II/grid/manifold.h +++ b/include/deal.II/grid/manifold.h @@ -613,16 +613,8 @@ public: * 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. + * for points close to one of the vertices, and avoiding tangent directions + * that are linearly dependent. */ virtual Tensor<1,spacedim> diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index 5940478ad5..54168ea40f 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -157,85 +157,54 @@ normal_vector (const Triangulation<3, 3>::face_iterator &face, { const int spacedim=3; Tensor<1,spacedim> t1,t2; + Tensor<1,spacedim> normal; - // Take the difference between p and all four vertices - unsigned int min_index=0; - double min_distance = (p-face->vertex(0)).norm_square(); - - for (unsigned int i=1; i<4; ++i) - { - const Tensor<1,spacedim> dp = p-face->vertex(i); - double distance = dp.norm_square(); - if (distance < min_distance) - { - min_index = i; - min_distance = distance; - } - } - // Verify we have a valid vertex index - AssertIndexRange(min_index, 4); - - // 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) + // 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<4 && j<4; ++j) { - // we are close to the face center: pick two consecutive vertices, - // but not the closest one. We make sure the direction is always - // the same. - if (min_index < 2) - { - t1 = get_tangent_vector(p, face->vertex(3)); - t2 = get_tangent_vector(p, face->vertex(2)); - } - else + // if p is too close to vertex i try again with different i and j + if ((p - face->vertex(i)).norm_square() < + std::numeric_limits::epsilon() * (p - face->vertex(j)).norm_square()) { - t1 = get_tangent_vector(p, face->vertex(0)); - t2 = get_tangent_vector(p, face->vertex(1)); - } - } - else - { - // we are closer to one of the vertices than to the - // center of the face - switch (min_index) - { - case 0: - { - t1 = get_tangent_vector(p, face->vertex(1)); - t2 = get_tangent_vector(p, face->vertex(2)); - break; - } - case 1: - { - t1 = get_tangent_vector(p, face->vertex(3)); - t2 = get_tangent_vector(p, face->vertex(0)); - break; - } - case 2: - { - t1 = get_tangent_vector(p, face->vertex(0)); - t2 = get_tangent_vector(p, face->vertex(3)); - break; - } - case 3: - { - t1 = get_tangent_vector(p, face->vertex(2)); - t2 = get_tangent_vector(p, face->vertex(1)); - break; - } - default: - Assert(false, ExcInternalError()); - break; + ++i; + continue; } + + // if p is too close to vertex j try again with different j + if ((p - face->vertex(j)).norm_square() < + std::numeric_limits::epsilon() * (p - face->vertex(i)).norm_square()) + continue; + + t1 = get_tangent_vector(p, face->vertex(i)); + t2 = get_tangent_vector(p, face->vertex(j)); + + normal = cross_product_3d(t1,t2); + + // if t1 and t2 are (nearly) linearly dependent try again with different j / t2 + if (normal.norm_square() < std::numeric_limits::epsilon() * + t1.norm_square() * t2.norm_square()) + continue; + + break; } - const Tensor<1,spacedim> normal = cross_product_3d(t1,t2); + Assert(normal.norm_square() >= std::numeric_limits::epsilon() * + t1.norm_square() * t2.norm_square(), + ExcMessage("Manifold::normal_vector was unable to find a suitable combination " + "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(); } diff --git a/tests/manifold/flat_manifold_08.cc b/tests/manifold/flat_manifold_08.cc new file mode 100644 index 0000000000..34a750e7fb --- /dev/null +++ b/tests/manifold/flat_manifold_08.cc @@ -0,0 +1,64 @@ +//---------------------------- manifold_id_01.cc --------------------------- +// Copyright (C) 2011 - 2015 by the mathLab team. +// +// This file is subject to LGPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- flat_manifold_01.cc --------------------------- + + +// Test that the flat manifold does what it should + +#include "../tests.h" + + +// all include files you need here +#include +#include +#include +#include +#include +#include + +// Helper function +template +void test(unsigned int ref=1) +{ + std::vector > vertices (GeometryInfo::vertices_per_cell); + std::vector > cells (1); + for (unsigned int i=0; i::vertices_per_cell; ++i) + cells[0].vertices[i] = i; + cells[0].material_id = 0; + + vertices[0] = Point(0,0,0); + vertices[1] = Point(1,0,0); + vertices[2] = Point(0.5,0.4,0); + vertices[3] = Point(1.5,0.4,0); + + vertices[4] = Point(0,0,1); + vertices[5] = Point(1,0,1); + vertices[6] = Point(0.5,0.4,1); + vertices[7] = Point(1.5,0.4,1); + + Triangulation tria; + tria.create_triangulation (vertices, cells, SubCellData()); + + typename Triangulation::active_cell_iterator + cell = tria.begin_active(); + + Point p1 (0.5,0,0); + deallog << "Normal vector of face 4: " << cell->get_manifold().normal_vector(cell->face(4),p1) << std::endl; + deallog << "Center of face 4: " << cell->face(4)->center() << std::endl; +} + +int main () +{ + initlog(); + + test<3,3>(); + + return 0; +} + diff --git a/tests/manifold/flat_manifold_08.output b/tests/manifold/flat_manifold_08.output new file mode 100644 index 0000000000..4f5536f658 --- /dev/null +++ b/tests/manifold/flat_manifold_08.output @@ -0,0 +1,3 @@ + +DEAL::Normal vector of face 4: 0.00000 0.00000 1.00000 +DEAL::Center of face 4: 0.750000 0.200000 0.00000 -- 2.39.5