From: Peter Munch <peterrmuench@gmail.com>
Date: Wed, 1 Jul 2020 10:56:46 +0000 (+0200)
Subject: Intoduce ReferenceCellInfo
X-Git-Tag: v9.3.0-rc1~1313^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10641%2Fhead;p=dealii.git

Intoduce ReferenceCellInfo
---

diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index 6f78f6e788..49e2033ccd 100644
--- a/include/deal.II/grid/reference_cell.h
+++ b/include/deal.II/grid/reference_cell.h
@@ -18,6 +18,8 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/geometry_info.h>
+
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -42,6 +44,347 @@ namespace ReferenceCell
     Hex     = 7,
     Invalid = static_cast<std::uint8_t>(-1)
   };
+
+  namespace internal
+  {
+    /**
+     * Check if the bit at position @p n in @p number is set.
+     */
+    inline static bool
+    get_bit(const unsigned char number, const unsigned int n)
+    {
+      AssertIndexRange(n, 8);
+
+      // source:
+      // https://stackoverflow.com/questions/47981/how-do-you-set-clear-and-toggle-a-single-bit
+      // "Checking a bit"
+      return (number >> n) & 1U;
+    }
+
+
+
+    /**
+     * Set the bit at position @p n in @p number to value @p x.
+     */
+    inline static void
+    set_bit(unsigned char &number, const unsigned int n, const bool x)
+    {
+      AssertIndexRange(n, 8);
+
+      // source:
+      // https://stackoverflow.com/questions/47981/how-do-you-set-clear-and-toggle-a-single-bit
+      // "Changing the nth bit to x"
+      number ^= (-static_cast<unsigned char>(x) ^ number) & (1U << n);
+    }
+
+    /**
+     * A namespace for geometric information on reference cells.
+     */
+    namespace Info
+    {
+      /**
+       * Interface to be used in TriaAccessor/TriaCellAccessor to access
+       * sub-entities of dimension d' of geometric entities of dimension d, with
+       * 0<=d'<d<=3.
+       */
+      struct Base
+      {
+        /**
+         * Number of vertices.
+         */
+        virtual unsigned int
+        n_vertices() const
+        {
+          Assert(false, ExcNotImplemented());
+          return 0;
+        }
+
+        /**
+         * Number of lines.
+         */
+        virtual unsigned int
+        n_lines() const
+        {
+          Assert(false, ExcNotImplemented());
+          return 0;
+        }
+
+
+        /**
+         * Number of faces.
+         */
+        virtual unsigned int
+        n_faces() const
+        {
+          Assert(false, ExcNotImplemented());
+          return 0;
+        }
+
+        /**
+         * Return an object that can be thought of as an array containing all
+         * indices from zero to n_vertices().
+         */
+        inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+        vertex_indices() const
+        {
+          return {0U, n_vertices()};
+        }
+
+        /**
+         * Return an object that can be thought of as an array containing all
+         * indices from zero to n_lines().
+         */
+        inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+        line_indices() const
+        {
+          return {0U, n_lines()};
+        }
+
+        /**
+         * Return an object that can be thought of as an array containing all
+         * indices from zero to n_faces().
+         */
+        inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+        face_indices() const
+        {
+          return {0U, n_faces()};
+        }
+
+        /**
+         * Standard decomposition of vertex index into face and face-vertex
+         * index.
+         */
+        virtual std::array<unsigned int, 2>
+        standard_vertex_to_face_and_vertex_index(
+          const unsigned int vertex) const
+        {
+          Assert(false, ExcNotImplemented());
+
+          (void)vertex;
+
+          return {0u, 0u};
+        }
+
+        /**
+         * Standard decomposition of line index into face and face-line index.
+         */
+        virtual std::array<unsigned int, 2>
+        standard_line_to_face_and_line_index(const unsigned int line) const
+        {
+          Assert(false, ExcNotImplemented());
+
+          (void)line;
+
+          return {0, 0};
+        }
+
+        /**
+         * Correct vertex index depending on face orientation.
+         */
+        virtual unsigned int
+        standard_to_real_face_vertex(const unsigned int  vertex,
+                                     const unsigned int  face,
+                                     const unsigned char face_orientation) const
+        {
+          Assert(false, ExcNotImplemented());
+
+          (void)vertex;
+          (void)face;
+          (void)face_orientation;
+
+          return 0;
+        }
+
+        /**
+         * Correct line index depending on face orientation.
+         */
+        virtual unsigned int
+        standard_to_real_face_line(const unsigned int  line,
+                                   const unsigned int  face,
+                                   const unsigned char face_orientation) const
+        {
+          Assert(false, ExcNotImplemented());
+
+          (void)line;
+          (void)face;
+          (void)face_orientation;
+
+          return 0;
+        }
+
+        /**
+         * Combine face and line orientation.
+         */
+        virtual bool
+        combine_face_and_line_orientation(
+          const unsigned int  line,
+          const unsigned char face_orientation,
+          const unsigned char line_orientation) const
+        {
+          Assert(false, ExcNotImplemented());
+
+          (void)line;
+          (void)face_orientation;
+          (void)line_orientation;
+
+          return true;
+        }
+      };
+
+
+      /**
+       * Base class for tensor-product geometric entities.
+       */
+      template <int dim>
+      struct TensorProductBase : Base
+      {
+        unsigned int
+        n_vertices() const override
+        {
+          return GeometryInfo<dim>::vertices_per_cell;
+        }
+
+        unsigned int
+        n_lines() const override
+        {
+          return GeometryInfo<dim>::lines_per_cell;
+        }
+
+        unsigned int
+        n_faces() const override
+        {
+          return GeometryInfo<dim>::faces_per_cell;
+        }
+      };
+
+
+
+      /*
+       * Vertex.
+       */
+      struct Vertex : public TensorProductBase<0>
+      {};
+
+
+
+      /*
+       * Line.
+       */
+      struct Line : public TensorProductBase<1>
+      {};
+
+
+
+      /**
+       * Quad.
+       */
+      struct Quad : public TensorProductBase<2>
+      {
+        std::array<unsigned int, 2>
+        standard_vertex_to_face_and_vertex_index(
+          const unsigned int vertex) const override
+        {
+          return GeometryInfo<2>::standard_quad_vertex_to_line_vertex_index(
+            vertex);
+        }
+
+        unsigned int
+        standard_to_real_face_vertex(
+          const unsigned int  vertex,
+          const unsigned int  face,
+          const unsigned char line_orientation) const override
+        {
+          (void)face;
+
+          return GeometryInfo<2>::standard_to_real_line_vertex(
+            vertex, line_orientation);
+        }
+      };
+
+
+
+      /**
+       * Hex.
+       */
+      struct Hex : public TensorProductBase<3>
+      {
+        std::array<unsigned int, 2>
+        standard_line_to_face_and_line_index(
+          const unsigned int line) const override
+        {
+          return GeometryInfo<3>::standard_hex_line_to_quad_line_index(line);
+        }
+
+        unsigned int
+        standard_to_real_face_line(
+          const unsigned int  line,
+          const unsigned int  face,
+          const unsigned char face_orientation) const override
+        {
+          (void)face;
+
+          return GeometryInfo<3>::standard_to_real_face_line(
+            line,
+            get_bit(face_orientation, 0),
+            get_bit(face_orientation, 2),
+            get_bit(face_orientation, 1));
+        }
+
+        bool
+        combine_face_and_line_orientation(
+          const unsigned int  line,
+          const unsigned char face_orientation_raw,
+          const unsigned char line_orientation) const override
+        {
+          static const bool bool_table[2][2][2][2] = {
+            {{{true, false},    // lines 0/1, face_orientation=false,
+                                // face_flip=false, face_rotation=false and true
+              {false, true}},   // lines 0/1, face_orientation=false,
+                                // face_flip=true, face_rotation=false and true
+             {{true, true},     // lines 0/1, face_orientation=true,
+                                // face_flip=false, face_rotation=false and true
+              {false, false}}}, // lines 0/1, face_orientation=true,
+                                // face_flip=true, face_rotation=false and true
+
+            {{{true, true}, // lines 2/3 ...
+              {false, false}},
+             {{true, false}, {false, true}}}};
+
+          const bool face_orientation = get_bit(face_orientation_raw, 0);
+          const bool face_flip        = get_bit(face_orientation_raw, 2);
+          const bool face_rotation    = get_bit(face_orientation_raw, 1);
+
+          return (
+            static_cast<bool>(line_orientation) ==
+            bool_table[line / 2][face_orientation][face_flip][face_rotation]);
+        }
+
+        std::array<unsigned int, 2>
+        standard_vertex_to_face_and_vertex_index(
+          const unsigned int vertex) const override
+        {
+          return GeometryInfo<3>::standard_hex_vertex_to_quad_vertex_index(
+            vertex);
+        }
+
+        unsigned int
+        standard_to_real_face_vertex(
+          const unsigned int  vertex,
+          const unsigned int  face,
+          const unsigned char face_orientation) const override
+        {
+          (void)face;
+
+          return GeometryInfo<3>::standard_to_real_face_vertex(
+            vertex,
+            get_bit(face_orientation, 0),
+            get_bit(face_orientation, 2),
+            get_bit(face_orientation, 1));
+        }
+      };
+
+    } // namespace Info
+  }   // namespace internal
 } // namespace ReferenceCell
 
 
diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h
index 07a449b705..9022adfade 100644
--- a/include/deal.II/grid/tria_accessor.h
+++ b/include/deal.II/grid/tria_accessor.h
@@ -1659,6 +1659,13 @@ public:
    * @}
    */
 
+protected:
+  /**
+   * Return additional information related to the current geometric entity
+   * type.
+   */
+  inline const ReferenceCell::internal::Info::Base &
+  reference_cell_info() const;
 
 private:
   /**
diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h
index 4865251c52..0f9408e7ab 100644
--- a/include/deal.II/grid/tria_accessor.templates.h
+++ b/include/deal.II/grid/tria_accessor.templates.h
@@ -536,13 +536,12 @@ namespace internal
       line_index(const TriaAccessor<3, 3, 3> &accessor, const unsigned int i)
       {
         const auto pair =
-          GeometryInfo<3>::standard_hex_line_to_quad_line_index(i);
+          accessor.reference_cell_info().standard_line_to_face_and_line_index(
+            i);
         const auto quad_index = pair[0];
-        const auto line_index = GeometryInfo<3>::standard_to_real_face_line(
-          pair[1],
-          accessor.face_orientation(quad_index),
-          accessor.face_flip(quad_index),
-          accessor.face_rotation(quad_index));
+        const auto line_index =
+          accessor.reference_cell_info().standard_to_real_face_line(
+            pair[1], pair[0], face_orientation_raw(accessor, quad_index));
 
         return accessor.quad(quad_index)->line_index(line_index);
       }
@@ -605,13 +604,22 @@ namespace internal
       face_orientation(const TriaAccessor<3, 3, 3> &accessor,
                        const unsigned int           face)
       {
-        return get_bit(
+        return ReferenceCell::internal::get_bit(
           accessor.tria->levels[accessor.present_level]->face_orientations
             [accessor.present_index * GeometryInfo<3>::faces_per_cell + face],
           0 /*=orientation_bit*/);
       }
 
 
+      inline static unsigned int
+      face_orientation_raw(const TriaAccessor<3, 3, 3> &accessor,
+                           const unsigned int           face)
+      {
+        AssertIndexRange(face, accessor.n_faces());
+        return accessor.tria->levels[accessor.present_level]->face_orientations
+          [accessor.present_index * GeometryInfo<3>::faces_per_cell + face];
+      }
+
 
       /**
        * Implementation of the function of some name in the mother class.
@@ -636,37 +644,6 @@ namespace internal
       }
 
 
-      /**
-       * Check if the bit at position @p n in @p number is set.
-       */
-      inline static bool
-      get_bit(const unsigned char number, const unsigned int n)
-      {
-        AssertIndexRange(n, 8);
-
-        // source:
-        // https://stackoverflow.com/questions/47981/how-do-you-set-clear-and-toggle-a-single-bit
-        // "Checking a bit"
-        return (number >> n) & 1U;
-      }
-
-
-
-      /**
-       * Set the bit at position @p n in @p number to value @p x.
-       */
-      inline static void
-      set_bit(unsigned char &number, const unsigned int n, const bool x)
-      {
-        AssertIndexRange(n, 8);
-
-        // source:
-        // https://stackoverflow.com/questions/47981/how-do-you-set-clear-and-toggle-a-single-bit
-        // "Changing the nth bit to x"
-        number ^= (-static_cast<unsigned char>(x) ^ number) & (1U << n);
-      }
-
-
 
       inline static bool
       face_flip(const TriaAccessor<3, 3, 3> &accessor, const unsigned int face)
@@ -677,7 +654,7 @@ namespace internal
                    ->face_orientations.size(),
                ExcInternalError());
 
-        return get_bit(
+        return ReferenceCell::internal::get_bit(
           accessor.tria->levels[accessor.present_level]->face_orientations
             [accessor.present_index * GeometryInfo<3>::faces_per_cell + face],
           2 /*=flip_bit*/);
@@ -713,7 +690,7 @@ namespace internal
                    ->face_orientations.size(),
                ExcInternalError());
 
-        return get_bit(
+        return ReferenceCell::internal::get_bit(
           accessor.tria->levels[accessor.present_level]->face_orientations
             [accessor.present_index * GeometryInfo<3>::faces_per_cell + face],
           1 /*=rotation_bit*/);
@@ -772,51 +749,17 @@ namespace internal
         AssertIndexRange(line, accessor.n_lines());
 
         const auto pair =
-          GeometryInfo<3>::standard_hex_line_to_quad_line_index(line);
-        const auto quad_index     = pair[0];
-        const auto std_line_index = pair[1];
-        const auto line_index     = GeometryInfo<3>::standard_to_real_face_line(
-          std_line_index,
-          accessor.face_orientation(quad_index),
-          accessor.face_flip(quad_index),
-          accessor.face_rotation(quad_index));
-
-        // now we got to the correct line and ask the quad for its
-        // line_orientation. however, if the face is rotated, it might be
-        // possible, that a standard orientation of the line with respect to
-        // the face corresponds to a non-standard orientation for the line with
-        // respect to the cell.
-        //
-        // set up a table indicating if the two standard orientations coincide
-        //
-        // first index: two pairs of lines 0(lines 0/1) and 1(lines 2/3)
-        //
-        // second index: face_orientation; 0: opposite normal, 1: standard
-        //
-        // third index: face_flip; 0: standard, 1: face rotated by 180 degrees
-        //
-        // forth index: face_rotation: 0: standard, 1: face rotated by 90
-        // degrees
-
-        static const bool bool_table[2][2][2][2] = {
-          {{{true, false},  // lines 0/1, face_orientation=false,
-                            // face_flip=false, face_rotation=false and true
-            {false, true}}, // lines 0/1, face_orientation=false,
-                            // face_flip=true, face_rotation=false and true
-           {{true, true}, // lines 0/1, face_orientation=true, face_flip=false,
-                          // face_rotation=false and true
-            {false, false}}}, // lines 0/1, face_orientation=true,
-                              // face_flip=true, face_rotation=false and true
-
-          {{{true, true}, // lines 2/3 ...
-            {false, false}},
-           {{true, false}, {false, true}}}};
-
-
-        return (accessor.quad(quad_index)->line_orientation(line_index) ==
-                bool_table[std_line_index / 2][accessor.face_orientation(
-                  quad_index)][accessor.face_flip(quad_index)]
-                          [accessor.face_rotation(quad_index)]);
+          accessor.reference_cell_info().standard_line_to_face_and_line_index(
+            line);
+        const auto quad_index = pair[0];
+        const auto line_index =
+          accessor.reference_cell_info().standard_to_real_face_line(
+            pair[1], pair[0], face_orientation_raw(accessor, quad_index));
+
+        return accessor.reference_cell_info().combine_face_and_line_orientation(
+          pair[1],
+          face_orientation_raw(accessor, quad_index),
+          accessor.quad(quad_index)->line_orientation(line_index));
       }
 
 
@@ -845,7 +788,7 @@ namespace internal
                  accessor.tria->levels[accessor.present_level]
                    ->face_orientations.size(),
                ExcInternalError());
-        set_bit(
+        ReferenceCell::internal::set_bit(
           accessor.tria->levels[accessor.present_level]->face_orientations
             [accessor.present_index * GeometryInfo<3>::faces_per_cell + face],
           0 /*=orientation_bit*/,
@@ -878,7 +821,7 @@ namespace internal
                    ->face_orientations.size(),
                ExcInternalError());
 
-        set_bit(
+        ReferenceCell::internal::set_bit(
           accessor.tria->levels[accessor.present_level]->face_orientations
             [accessor.present_index * GeometryInfo<3>::faces_per_cell + face],
           2 /*=flip_bit*/,
@@ -911,7 +854,7 @@ namespace internal
                    ->face_orientations.size(),
                ExcInternalError());
 
-        set_bit(
+        ReferenceCell::internal::set_bit(
           accessor.tria->levels[accessor.present_level]->face_orientations
             [accessor.present_index * GeometryInfo<3>::faces_per_cell + face],
           1 /*=rotation_bit*/,
@@ -991,11 +934,12 @@ namespace internal
       vertex_index(const TriaAccessor<2, dim, spacedim> &accessor,
                    const unsigned int                    corner)
       {
-        const auto pair =
-          GeometryInfo<2>::standard_quad_vertex_to_line_vertex_index(corner);
-        const auto line_index   = pair[0];
-        const auto vertex_index = GeometryInfo<2>::standard_to_real_line_vertex(
-          pair[1], accessor.line_orientation(line_index));
+        const auto pair = accessor.reference_cell_info()
+                            .standard_vertex_to_face_and_vertex_index(corner);
+        const auto line_index = pair[0];
+        const auto vertex_index =
+          accessor.reference_cell_info().standard_to_real_face_vertex(
+            pair[1], pair[0], accessor.line_orientation(line_index));
 
         return accessor.line(line_index)->vertex_index(vertex_index);
       }
@@ -1006,14 +950,12 @@ namespace internal
       vertex_index(const TriaAccessor<3, 3, 3> &accessor,
                    const unsigned int           corner)
       {
-        const auto pair =
-          GeometryInfo<3>::standard_hex_vertex_to_quad_vertex_index(corner);
-        const auto face_index   = pair[0];
-        const auto vertex_index = GeometryInfo<3>::standard_to_real_face_vertex(
-          pair[1],
-          accessor.face_orientation(face_index),
-          accessor.face_flip(face_index),
-          accessor.face_rotation(face_index));
+        const auto pair = accessor.reference_cell_info()
+                            .standard_vertex_to_face_and_vertex_index(corner);
+        const auto face_index = pair[0];
+        const auto vertex_index =
+          accessor.reference_cell_info().standard_to_real_face_vertex(
+            pair[1], pair[0], face_orientation_raw(accessor, face_index));
 
         return accessor.quad(face_index)->vertex_index(vertex_index);
       }
@@ -2165,7 +2107,7 @@ template <int structdim, int dim, int spacedim>
 unsigned int
 TriaAccessor<structdim, dim, spacedim>::n_vertices() const
 {
-  return GeometryInfo<structdim>::vertices_per_cell;
+  return this->reference_cell_info().n_vertices();
 }
 
 
@@ -2174,7 +2116,7 @@ template <int structdim, int dim, int spacedim>
 unsigned int
 TriaAccessor<structdim, dim, spacedim>::n_lines() const
 {
-  return GeometryInfo<structdim>::lines_per_cell;
+  return this->reference_cell_info().n_lines();
 }
 
 
@@ -2185,7 +2127,7 @@ TriaAccessor<structdim, dim, spacedim>::n_faces() const
 {
   AssertDimension(structdim, dim);
 
-  return GeometryInfo<structdim>::faces_per_cell;
+  return this->reference_cell_info().n_faces();
 }
 
 
@@ -2216,6 +2158,34 @@ TriaAccessor<structdim, dim, spacedim>::face_indices() const
 }
 
 
+
+template <int structdim, int dim, int spacedim>
+inline const ReferenceCell::internal::Info::Base &
+TriaAccessor<structdim, dim, spacedim>::reference_cell_info() const
+{
+  static ReferenceCell::internal::Info::Base   gei_invalid;
+  static ReferenceCell::internal::Info::Vertex gei_vertex;
+  static ReferenceCell::internal::Info::Line   gei_line;
+  static ReferenceCell::internal::Info::Quad   gei_quad;
+  static ReferenceCell::internal::Info::Hex    gei_hex;
+
+  switch (structdim) // TODO: use ReferenceCell::Type
+    {
+      case 0:
+        return gei_vertex;
+      case 1:
+        return gei_line;
+      case 2:
+        return gei_quad;
+      case 3:
+        return gei_hex;
+      default:
+        Assert(false, StandardExceptions::ExcNotImplemented());
+        return gei_invalid;
+    }
+}
+
+
 /*----------------- Functions: TriaAccessor<0,dim,spacedim> -----------------*/
 
 template <int dim, int spacedim>