From 03e99f05657a935c118489f612c063e48f511f1c Mon Sep 17 00:00:00 2001
From: Jiaqi Zhang <jiaqi2@clemson.edu>
Date: Mon, 19 Sep 2022 21:46:17 +0000
Subject: [PATCH] fix normal on a coarse mesh

---
 source/grid/manifold.cc                     |  42 +++++++-
 tests/manifold/spherical_manifold_14.cc     | 103 ++++++++++++++++++++
 tests/manifold/spherical_manifold_14.output |   2 +
 3 files changed, 144 insertions(+), 3 deletions(-)
 create mode 100644 tests/manifold/spherical_manifold_14.cc
 create mode 100644 tests/manifold/spherical_manifold_14.output

diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc
index cf833fdf7a..5a36d71341 100644
--- a/source/grid/manifold.cc
+++ b/source/grid/manifold.cc
@@ -206,9 +206,45 @@ Manifold<3, 3>::normal_vector(const Triangulation<3, 3>::face_iterator &face,
          ExcMessage("The search for possible directions did not succeed."));
 
   // Compute tangents and normal for selected vertices
-  Tensor<1, spacedim> t1     = get_tangent_vector(p, vertices[first_index]);
-  Tensor<1, spacedim> t2     = get_tangent_vector(p, vertices[second_index]);
-  Tensor<1, spacedim> normal = cross_product_3d(t1, t2);
+  Tensor<1, spacedim> t1;
+  Tensor<1, spacedim> t2;
+  Tensor<1, spacedim> normal;
+
+  bool              done = false;
+  std::vector<bool> tested_vertices(vertices.size(), false);
+  tested_vertices[first_index]  = true;
+  tested_vertices[second_index] = true;
+
+  do
+    {
+      // Compute tangents and normal for selected vertices
+      t1     = get_tangent_vector(p, vertices[first_index]);
+      t2     = get_tangent_vector(p, vertices[second_index]);
+      normal = cross_product_3d(t1, t2);
+
+      // if normal is zero, try some other combination of veritices
+      if (normal.norm_square() < 1e4 * std::numeric_limits<double>::epsilon() *
+                                   t1.norm_square() * t2.norm_square())
+        {
+          // See if we have tested already all vertices
+          auto first_false =
+            std::find(tested_vertices.begin(), tested_vertices.end(), false);
+          if (first_false == tested_vertices.end())
+            {
+              done = true;
+            }
+          else
+            {
+              *first_false = true;
+              second_index = first_false - tested_vertices.begin();
+            }
+        }
+      else
+        {
+          done = true;
+        }
+    }
+  while (!done);
 
   Assert(
     normal.norm_square() > 1e4 * std::numeric_limits<double>::epsilon() *
diff --git a/tests/manifold/spherical_manifold_14.cc b/tests/manifold/spherical_manifold_14.cc
new file mode 100644
index 0000000000..105dec6617
--- /dev/null
+++ b/tests/manifold/spherical_manifold_14.cc
@@ -0,0 +1,103 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2012 - 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// A test for Issue #14183. On a very coarse mesh, the
+// vertices chosen to compute the normal vector may
+// lead to two parallel tangent vectors, resulting
+// in a zero normal vector. This can be fixed by
+// choosing other combinations of vertices when the
+// normal vector is too close to a zero vector.
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/vector_tools_constraints.h>
+
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+  deallog.get_file_stream().precision(7);
+  deallog.get_file_stream().setf(std::ios::fixed);
+
+  {
+    const unsigned int dim = 3;
+    Triangulation<dim> triangulation(
+      Triangulation<dim>::limit_level_difference_at_vertices);
+    GridGenerator::half_hyper_shell(triangulation, Point<dim>(), 0.5, 1.);
+
+    std::set<types::boundary_id> no_flux_boundary{0};
+    MappingQ<dim>                mapping(2);
+
+    FESystem<dim>   fe(FE_Q<dim>(2), dim);
+    DoFHandler<dim> dofh(triangulation);
+
+    dofh.distribute_dofs(fe);
+
+    AffineConstraints<double> constraints;
+
+
+    const unsigned int                 face_no = 0;
+    const std::vector<Point<dim - 1>> &unit_support_points =
+      fe.get_unit_face_support_points(face_no);
+
+    Assert(unit_support_points.size() == fe.n_dofs_per_face(),
+           ExcInternalError());
+
+
+    Quadrature<dim - 1> face_quadrature(unit_support_points);
+
+    FEFaceValues<dim> fe_face_values(mapping,
+                                     fe,
+                                     face_quadrature,
+                                     update_quadrature_points |
+                                       update_normal_vectors);
+
+
+    std::set<types::boundary_id>::iterator b_id;
+
+    for (const auto &cell : dofh.active_cell_iterators())
+      if (!cell->is_artificial())
+        for (const unsigned int face_no : cell->face_indices())
+          if ((b_id = no_flux_boundary.find(
+                 cell->face(face_no)->boundary_id())) != no_flux_boundary.end())
+            {
+              fe_face_values.reinit(cell, face_no);
+              for (unsigned int i = 0; i < face_quadrature.size(); ++i)
+                Tensor<1, dim> normal_vector =
+                  (cell->face(face_no)->get_manifold().normal_vector(
+                    cell->face(face_no), fe_face_values.quadrature_point(i)));
+            }
+
+    deallog << " OK! " << std::endl;
+  }
+}
diff --git a/tests/manifold/spherical_manifold_14.output b/tests/manifold/spherical_manifold_14.output
new file mode 100644
index 0000000000..73ff79ea95
--- /dev/null
+++ b/tests/manifold/spherical_manifold_14.output
@@ -0,0 +1,2 @@
+
+DEAL:: OK! 
-- 
2.39.5