Fix ReferenceCell::standard_vs_true_line_orientation() for simplices

doc/news/changes/minor/20231010Munch

@@ -0,0 +1,4 @@
+Fixed: The querying of line orientations in the case
+of tetrahedra has been fixed.
+<br>
+(Peter Munch, David Wells, 2023/10/10)

include/deal.II/grid/reference_cell.h

@@ -522,6 +522,7 @@ class ReferenceCell
    */
   bool
   standard_vs_true_line_orientation(const unsigned int  line,
+                                    const unsigned int  face,
                                     const unsigned char face_orientation,
                                     const bool          line_orientation) const;
 
@@ -2598,6 +2599,7 @@ ReferenceCell::n_face_orientations(const unsigned int face_no) const
 inline bool
 ReferenceCell::standard_vs_true_line_orientation(
   const unsigned int  line,
+  const unsigned int  face,
   const unsigned char combined_face_orientation,
   const bool          line_orientation) const
 {
@@ -2610,6 +2612,26 @@ ReferenceCell::standard_vs_true_line_orientation(
       return (line_orientation ==
               bool_table[line / 2][combined_face_orientation]);
     }
+  else if (*this == ReferenceCells::Tetrahedron)
+    {
+      static constexpr unsigned int X = numbers::invalid_unsigned_int;
+      static constexpr dealii::ndarray<unsigned int, 4, 3> combined_lines{
+        {{{0, 0, 0}}, {{X, 0, 1}}, {{X, 0, X}}, {{X, X, X}}}};
+
+      const auto combined_line = combined_lines[face][line];
+
+      Assert(combined_line != X,
+             ExcMessage(
+               "This function can only be called for following face-line "
+               "combinations: (0,0), (0,1), (0,2), (1,1), (1,2), (2,1),"));
+
+      static constexpr dealii::ndarray<bool, 2, 6> bool_table{
+        {{{false, true, false, true, false, true}},
+         {{true, false, true, false, true, false}}}};
+
+      return (line_orientation ==
+              bool_table[combined_line][combined_face_orientation]);
+    }
   else
     // TODO: This might actually be wrong for some of the other
     // kinds of objects. We should check this

include/deal.II/grid/tria_accessor.templates.h

@@ -854,6 +854,7 @@ namespace internal
         // Then query how that line is oriented within that face:
         return accessor.reference_cell().standard_vs_true_line_orientation(
           line_index,
+          face_index,
           combined_face_orientation(accessor, face_index),
           accessor.quad(face_index)->line_orientation(line_within_face_index));
       }
@@ -1170,13 +1171,16 @@ namespace internal
                 const auto                quad = cell.quad(f);
                 const std::array<bool, 4> my_orientations{
                   {ref_cell.standard_vs_true_line_orientation(
-                     0, orientation, quad->line_orientation(my_indices[0])),
+                     0, f, orientation, quad->line_orientation(my_indices[0])),
                    ref_cell.standard_vs_true_line_orientation(
-                     1, orientation, quad->line_orientation(my_indices[1])),
+                     1, f, orientation, quad->line_orientation(my_indices[1])),
                    ref_cell.standard_vs_true_line_orientation(
-                     2, orientation, quad->line_orientation(my_indices[2])),
+                     2, f, orientation, quad->line_orientation(my_indices[2])),
                    ref_cell.standard_vs_true_line_orientation(
-                     3, orientation, quad->line_orientation(my_indices[3]))}};
+                     3,
+                     f,
+                     orientation,
+                     quad->line_orientation(my_indices[3]))}};
                 for (unsigned int l = 0; l < 4; ++l)
                   line_orientations[4 * (f - 4) + l] = my_orientations[l];
               }
@@ -1193,9 +1197,12 @@ namespace internal
                 const auto                quad = cell.quad(f);
                 const std::array<bool, 2> my_orientations{
                   {ref_cell.standard_vs_true_line_orientation(
-                     0, orientation, quad->line_orientation(my_indices[0])),
+                     0, f, orientation, quad->line_orientation(my_indices[0])),
                    ref_cell.standard_vs_true_line_orientation(
-                     1, orientation, quad->line_orientation(my_indices[1]))}};
+                     1,
+                     f,
+                     orientation,
+                     quad->line_orientation(my_indices[1]))}};
                 line_orientations[8 + f]  = my_orientations[0];
                 line_orientations[10 + f] = my_orientations[1];
               }
@@ -1214,27 +1221,33 @@ namespace internal
                ref_cell.standard_to_real_face_line(2, 1, orientations[1]),
                ref_cell.standard_to_real_face_line(1, 2, orientations[2])}};
             line_orientations[0] = ref_cell.standard_vs_true_line_orientation(
+              0,
               0,
               orientations[0],
               cell.quad(0)->line_orientation(my_indices[0]));
             line_orientations[1] = ref_cell.standard_vs_true_line_orientation(
               1,
+              0,
               orientations[0],
               cell.quad(0)->line_orientation(my_indices[1]));
             line_orientations[2] = ref_cell.standard_vs_true_line_orientation(
               2,
+              0,
               orientations[0],
               cell.quad(0)->line_orientation(my_indices[2]));
             line_orientations[3] = ref_cell.standard_vs_true_line_orientation(
+              1,
               1,
               orientations[1],
               cell.quad(1)->line_orientation(my_indices[3]));
             line_orientations[4] = ref_cell.standard_vs_true_line_orientation(
               2,
+              1,
               orientations[1],
               cell.quad(1)->line_orientation(my_indices[4]));
             line_orientations[5] = ref_cell.standard_vs_true_line_orientation(
               1,
+              2,
               orientations[2],
               cell.quad(2)->line_orientation(my_indices[5]));
           }

tests/simplex/orientation_04.cc

@@ -0,0 +1,187 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2022 - 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.
+ *
+ * ---------------------------------------------------------------------
+ */
+
+// Test a mesh with two tetrahedra for all possible orientations. Similar to
+// orientation_02 but also checks that line orientations are correct.
+
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools_integrate_difference.h>
+#include <deal.II/numerics/vector_tools_interpolate.h>
+#include <deal.II/numerics/vector_tools_project.h>
+
+#include "../tests.h"
+
+std::tuple<unsigned int, unsigned int>
+create_triangulation(const std::vector<Point<3>>    &vertices_,
+                     const std::vector<CellData<3>> &cell_data_,
+                     const unsigned int              face_n,
+                     const unsigned int              n_permuations,
+                     Triangulation<3>               &tria)
+{
+  const ReferenceCell ref_cell  = ReferenceCells::Tetrahedron;
+  auto                vertices  = vertices_;
+  auto                cell_data = cell_data_;
+
+  Point<3> extra_vertex;
+  for (unsigned int i = 0; i < 3; ++i)
+    extra_vertex += ref_cell.template vertex<3>(ref_cell.face_to_cell_vertices(
+      face_n, i, ReferenceCell::default_combined_face_orientation()));
+
+  extra_vertex /= 3.0;
+  extra_vertex += ref_cell.template unit_normal_vectors<3>(face_n);
+
+  vertices.push_back(extra_vertex);
+
+  cell_data.emplace_back();
+  cell_data.back().vertices.resize(0);
+  for (unsigned int i = 0; i < 3; ++i)
+    cell_data.back().vertices.push_back(ref_cell.face_to_cell_vertices(
+      face_n, i, ref_cell.default_combined_face_orientation()));
+  cell_data.back().vertices.push_back(ref_cell.n_vertices());
+  std::sort(cell_data.back().vertices.begin(), cell_data.back().vertices.end());
+
+  unsigned int permutation_n = 0;
+  do
+    {
+      tria.clear();
+      tria.create_triangulation(vertices, cell_data, SubCellData());
+      ++permutation_n;
+    }
+  while ((permutation_n < n_permuations) &&
+         std::next_permutation(cell_data.back().vertices.begin(),
+                               cell_data.back().vertices.end()));
+
+  const auto cell = tria.begin();
+
+  const auto face = cell->face(face_n);
+
+  auto ncell = tria.begin();
+  ncell++;
+  ncell->face(face_n);
+
+  unsigned int nf = 0;
+  for (; nf < ref_cell.n_faces(); ++nf)
+    if (ncell->face(nf) == face)
+      break;
+
+  return {nf, ncell->combined_face_orientation(nf)};
+}
+
+
+
+void
+test()
+{
+  const unsigned int dim = 3;
+
+  double previous_error = 1.0;
+
+  for (unsigned int f = 0; f < 4; ++f)
+    {
+      for (unsigned int r = 0; r < 24; ++r)
+        {
+          unsigned int orientation = r;
+          unsigned int face_no     = f;
+
+          Triangulation<dim> tria;
+
+          Triangulation<3> dummy;
+          GridGenerator::reference_cell(dummy, ReferenceCells::Tetrahedron);
+
+          auto vertices = dummy.get_vertices();
+
+          std::vector<CellData<3>> cells;
+
+          {
+            CellData<3> cell;
+            cell.vertices    = {0, 1, 2, 3};
+            cell.material_id = 0;
+            cells.push_back(cell);
+          }
+
+          std::tie(face_no, orientation) =
+            create_triangulation(vertices, cells, face_no, r, tria);
+
+          bool success = true;
+
+          auto cell = tria.begin();
+          cell++;
+
+          std::vector<unsigned int> verticess;
+
+          for (const auto v : cell->vertex_indices())
+            verticess.emplace_back(cell->vertex_index(v));
+
+          for (unsigned int ll = 0; ll < 3; ++ll)
+            {
+              const unsigned int l =
+                cell->reference_cell().face_to_cell_lines(face_no, ll, 1);
+
+              const auto orientation_exp = cell->line_orientation(l);
+
+              std::pair<unsigned int, unsigned int> p0;
+              p0.first =
+                verticess[cell->reference_cell().line_to_cell_vertices(l, 0)];
+              p0.second =
+                verticess[cell->reference_cell().line_to_cell_vertices(l, 1)];
+
+              std::pair<unsigned int, unsigned int> p1;
+              p1.first  = cell->line(l)->vertex_index(0);
+              p1.second = cell->line(l)->vertex_index(1);
+
+              if (orientation_exp == false)
+                std::swap(p1.first, p1.second);
+
+              success &= (p0 == p1);
+            }
+
+          if (success)
+            deallog << "x ";
+          else
+            deallog << "o ";
+        }
+    }
+
+  deallog << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+
+  test();
+}

tests/simplex/orientation_04.output

@@ -0,0 +1,2 @@
+
+DEAL::x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x