From c68c758827494a753904d9928df1da2bc7e32273 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 1 Jul 2020 12:56:46 +0200 Subject: [PATCH] Intoduce ReferenceCellInfo --- include/deal.II/grid/reference_cell.h | 343 ++++++++++++++++++ include/deal.II/grid/tria_accessor.h | 7 + .../deal.II/grid/tria_accessor.templates.h | 178 ++++----- 3 files changed, 424 insertions(+), 104 deletions(-) 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 +#include + DEAL_II_NAMESPACE_OPEN @@ -42,6 +44,347 @@ namespace ReferenceCell Hex = 7, Invalid = static_cast(-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(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' + 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 + 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 + face_indices() const + { + return {0U, n_faces()}; + } + + /** + * Standard decomposition of vertex index into face and face-vertex + * index. + */ + virtual std::array + 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 + 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 + struct TensorProductBase : Base + { + unsigned int + n_vertices() const override + { + return GeometryInfo::vertices_per_cell; + } + + unsigned int + n_lines() const override + { + return GeometryInfo::lines_per_cell; + } + + unsigned int + n_faces() const override + { + return GeometryInfo::faces_per_cell; + } + }; + + + + /* + * Vertex. + */ + struct Vertex : public TensorProductBase<0> + {}; + + + + /* + * Line. + */ + struct Line : public TensorProductBase<1> + {}; + + + + /** + * Quad. + */ + struct Quad : public TensorProductBase<2> + { + std::array + 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 + 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(line_orientation) == + bool_table[line / 2][face_orientation][face_flip][face_rotation]); + } + + std::array + 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(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 unsigned int TriaAccessor::n_vertices() const { - return GeometryInfo::vertices_per_cell; + return this->reference_cell_info().n_vertices(); } @@ -2174,7 +2116,7 @@ template unsigned int TriaAccessor::n_lines() const { - return GeometryInfo::lines_per_cell; + return this->reference_cell_info().n_lines(); } @@ -2185,7 +2127,7 @@ TriaAccessor::n_faces() const { AssertDimension(structdim, dim); - return GeometryInfo::faces_per_cell; + return this->reference_cell_info().n_faces(); } @@ -2216,6 +2158,34 @@ TriaAccessor::face_indices() const } + +template +inline const ReferenceCell::internal::Info::Base & +TriaAccessor::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 -- 2.39.5