diff --git a/include/Makefile.in b/include/Makefile.in index 00edd929c8..140ebb4df7 100644 --- a/include/Makefile.in +++ b/include/Makefile.in @@ -715,6 +715,7 @@ include_HEADERS = \ fe/fe_lagrange_shape_1D.h \ fe/fe_macro.h \ fe/fe_map.h \ + fe/fe_reference_element_traits.h \ fe/fe_transformation_base.h \ fe/fe_type.h \ fe/fe_xyz_map.h \ diff --git a/include/fe/fe_reference_element_traits.h b/include/fe/fe_reference_element_traits.h new file mode 100644 index 0000000000..1af7be5bc9 --- /dev/null +++ b/include/fe/fe_reference_element_traits.h @@ -0,0 +1,170 @@ +// The libMesh Finite Element Library. +// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// This library is free software; you can 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. + +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. + +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +#ifndef LIBMESH_FE_REFERENCE_ELEMENT_TRAITS_H +#define LIBMESH_FE_REFERENCE_ELEMENT_TRAITS_H + +#include "libmesh/libmesh.h" // invalid_uint + +namespace libMesh +{ + +/** + * The second-order elements' node maps are not independent facts. + * Their edge_nodes_map is the first-order equivalent's edge_nodes_map + * with the mid-edge node appended, and mid-edge nodes are numbered + * after the vertices in edge order. Their side_nodes_map row is the + * first-order equivalent's vertex row, then the mid-edge node of each + * consecutive vertex pair, then the face node if the side has one. + * + * The functions here compute those tables at compile time from the + * first-order tables, so a second-order element class stores neither: + * it holds the derived tables as static constexpr members and exposes + * them through its usual side_nodes_map / edge_nodes_map names as + * references to the underlying arrays. Only the first-order tables + * and the face-node rules are written by hand. + */ + +/** + * The value the element classes pad short side_nodes_map rows with, + * e.g. the triangular sides of a Prism. + */ +static constexpr unsigned int unused_side_node = 99; + +/** + * A fixed-size table that constexpr functions can return. Element + * classes bind their side_nodes_map / edge_nodes_map references to + * \p values. + */ +template +struct NodeMapTable +{ + unsigned int values[Rows][Cols]; +}; + +/** + * The face-node rule for elements with no face nodes. + */ +constexpr unsigned int no_face_node (const unsigned int) +{ + return invalid_uint; +} + +/** + * \returns The edge_nodes_map of the second-order element whose + * first-order equivalent is \p FirstOrder: each edge's two vertices, then + * its mid-edge node, numbered after the vertices in edge order. + */ +template +constexpr NodeMapTable +derived_edge_nodes () +{ + NodeMapTable t {}; + for (unsigned int e = 0; e != Edges; ++e) + { + t.values[e][0] = FirstOrder::edge_nodes_map[e][0]; + t.values[e][1] = FirstOrder::edge_nodes_map[e][1]; + t.values[e][2] = FirstOrder::num_nodes + e; + } + return t; +} + +/** + * \returns The number of vertices on side \p s of a first-order + * element, i.e. the entries of its side_nodes_map row that aren't + * padding. + */ +template +constexpr unsigned int n_side_vertices (const unsigned int s) +{ + unsigned int n = 0; + for (unsigned int k = 0; k != FirstOrder::nodes_per_side; ++k) + if (FirstOrder::side_nodes_map[s][k] != unused_side_node) + ++n; + return n; +} + +/** + * \returns The third entry of the row of \p edges joining vertices \p a + * and \p b, i.e. that edge's mid-edge node. + */ +template +constexpr unsigned int mid_edge_node (const unsigned int (&edges)[Edges][3], + const unsigned int a, + const unsigned int b) +{ + for (unsigned int e = 0; e != Edges; ++e) + if ((edges[e][0] == a && edges[e][1] == b) || + (edges[e][0] == b && edges[e][1] == a)) + return edges[e][2]; + return invalid_uint; +} + +/** + * \returns The side_nodes_map of a 3D second-order element with + * \p Sides sides of up to \p Cols nodes each, whose first-order + * equivalent is \p FirstOrder, whose edge_nodes_map is \p edges, and whose + * \p face_node(s) is the node at the center of side \p s (or + * \p invalid_uint if there is none). + */ +template +constexpr NodeMapTable +derived_side_nodes (const unsigned int (&edges)[Edges][3], + FaceNode face_node) +{ + NodeMapTable t {}; + for (unsigned int s = 0; s != Sides; ++s) + { + const unsigned int nv = n_side_vertices(s); + unsigned int n = 0; + for (unsigned int k = 0; k != nv; ++k) + t.values[s][n++] = FirstOrder::side_nodes_map[s][k]; + for (unsigned int k = 0; k != nv; ++k) + t.values[s][n++] = mid_edge_node(edges, + FirstOrder::side_nodes_map[s][k], + FirstOrder::side_nodes_map[s][(k+1) % nv]); + if (face_node(s) != invalid_uint) + t.values[s][n++] = face_node(s); + for (; n != Cols; ++n) + t.values[s][n] = unused_side_node; + } + return t; +} + +/** + * \returns The side_nodes_map of a 2D second-order element whose + * first-order equivalent is \p FirstOrder: each side's two vertices, then + * its mid-side node, numbered after the vertices in side order. + */ +template +constexpr NodeMapTable +derived_side_nodes () +{ + NodeMapTable t {}; + for (unsigned int s = 0; s != Sides; ++s) + { + t.values[s][0] = FirstOrder::side_nodes_map[s][0]; + t.values[s][1] = FirstOrder::side_nodes_map[s][1]; + t.values[s][2] = FirstOrder::num_nodes + s; + } + return t; +} + +} // namespace libMesh + +#endif // LIBMESH_FE_REFERENCE_ELEMENT_TRAITS_H diff --git a/include/geom/cell_hex20.h b/include/geom/cell_hex20.h index 70c37c23f7..f857d3e018 100644 --- a/include/geom/cell_hex20.h +++ b/include/geom/cell_hex20.h @@ -22,6 +22,8 @@ // Local includes #include "libmesh/cell_hex.h" +#include "libmesh/cell_hex8.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -217,16 +219,18 @@ class Hex20 final : public Hex static const int nodes_per_edge = 3; /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or + * side to element node numbers. They are derived from the + * first-order Hex8 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; - - /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to - * element node numbers. - */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr NodeMapTable + _edge_nodes = derived_edge_nodes(); + static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values; + + static constexpr NodeMapTable + _side_nodes = derived_side_nodes + (_edge_nodes.values, no_face_node); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; /** * A specialization for computing the volume of a Hex20. diff --git a/include/geom/cell_hex27.h b/include/geom/cell_hex27.h index 0777540a30..1dd50c3079 100644 --- a/include/geom/cell_hex27.h +++ b/include/geom/cell_hex27.h @@ -22,6 +22,8 @@ // Local includes #include "libmesh/cell_hex.h" +#include "libmesh/cell_hex8.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -232,16 +234,18 @@ class Hex27 final : public Hex static const int nodes_per_edge = 3; /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or + * side to element node numbers. They are derived from the + * first-order Hex8 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; - - /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to - * element node numbers. - */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr NodeMapTable + _edge_nodes = derived_edge_nodes(); + static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values; + + static constexpr NodeMapTable + _side_nodes = derived_side_nodes + (_edge_nodes.values, [](unsigned int s) { return 20 + s; }); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; /** * A specialization for computing the volume of a Hex27. diff --git a/include/geom/cell_hex8.h b/include/geom/cell_hex8.h index 9abee93171..cf77eb4594 100644 --- a/include/geom/cell_hex8.h +++ b/include/geom/cell_hex8.h @@ -170,13 +170,35 @@ class Hex8 final : public Hex * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to * element node numbers. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; + static constexpr unsigned int side_nodes_map[num_sides][nodes_per_side] = + { + {0, 3, 2, 1}, // Side 0 + {0, 1, 5, 4}, // Side 1 + {1, 2, 6, 5}, // Side 2 + {2, 3, 7, 6}, // Side 3 + {3, 0, 4, 7}, // Side 4 + {4, 5, 6, 7} // Side 5 + }; /** * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to * element node numbers. */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr unsigned int edge_nodes_map[num_edges][nodes_per_edge] = + { + {0, 1}, // Edge 0 + {1, 2}, // Edge 1 + {2, 3}, // Edge 2 + {0, 3}, // Edge 3 + {0, 4}, // Edge 4 + {1, 5}, // Edge 5 + {2, 6}, // Edge 6 + {3, 7}, // Edge 7 + {4, 5}, // Edge 8 + {5, 6}, // Edge 9 + {6, 7}, // Edge 10 + {4, 7} // Edge 11 + }; /** * Class static helper function that computes the centroid of a diff --git a/include/geom/cell_prism15.h b/include/geom/cell_prism15.h index 68374b06f6..8227ca142c 100644 --- a/include/geom/cell_prism15.h +++ b/include/geom/cell_prism15.h @@ -22,6 +22,8 @@ // Local includes #include "libmesh/cell_prism.h" +#include "libmesh/cell_prism6.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -222,16 +224,18 @@ class Prism15 final : public Prism static const int nodes_per_edge = 3; /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or + * side to element node numbers. They are derived from the + * first-order Prism6 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; - - /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to - * element node numbers. - */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr NodeMapTable + _edge_nodes = derived_edge_nodes(); + static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values; + + static constexpr NodeMapTable + _side_nodes = derived_side_nodes + (_edge_nodes.values, no_face_node); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; /** * A specialization for computing the volume of a Prism15. diff --git a/include/geom/cell_prism18.h b/include/geom/cell_prism18.h index 530f6f3797..646a76675a 100644 --- a/include/geom/cell_prism18.h +++ b/include/geom/cell_prism18.h @@ -22,6 +22,8 @@ // Local includes #include "libmesh/cell_prism.h" +#include "libmesh/cell_prism6.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -237,16 +239,18 @@ class Prism18 final : public Prism static const int nodes_per_edge = 3; /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or + * side to element node numbers. They are derived from the + * first-order Prism6 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; - - /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to - * element node numbers. - */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr NodeMapTable + _edge_nodes = derived_edge_nodes(); + static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values; + + static constexpr NodeMapTable + _side_nodes = derived_side_nodes + (_edge_nodes.values, [](unsigned int s) { return (s >= 1 && s <= 3) ? 14 + s : invalid_uint; }); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; /** * A specialization for computing the volume of a Prism18. diff --git a/include/geom/cell_prism20.h b/include/geom/cell_prism20.h index af1bedf263..b9a2eb335d 100644 --- a/include/geom/cell_prism20.h +++ b/include/geom/cell_prism20.h @@ -22,6 +22,8 @@ // Local includes #include "libmesh/cell_prism.h" +#include "libmesh/cell_prism6.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -242,16 +244,18 @@ class Prism20 final : public Prism static const int nodes_per_edge = 3; /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or + * side to element node numbers. They are derived from the + * first-order Prism6 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; - - /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to - * element node numbers. - */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr NodeMapTable + _edge_nodes = derived_edge_nodes(); + static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values; + + static constexpr NodeMapTable + _side_nodes = derived_side_nodes + (_edge_nodes.values, [](unsigned int s) { return s == 0 ? 18u : s == 4 ? 19u : 14 + s; }); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; virtual void permute(unsigned int perm_num) override final; diff --git a/include/geom/cell_prism21.h b/include/geom/cell_prism21.h index 894f86789f..2ceb256775 100644 --- a/include/geom/cell_prism21.h +++ b/include/geom/cell_prism21.h @@ -22,6 +22,8 @@ // Local includes #include "libmesh/cell_prism.h" +#include "libmesh/cell_prism6.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -245,16 +247,18 @@ class Prism21 final : public Prism static const int nodes_per_edge = 3; /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or + * side to element node numbers. They are derived from the + * first-order Prism6 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; - - /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to - * element node numbers. - */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr NodeMapTable + _edge_nodes = derived_edge_nodes(); + static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values; + + static constexpr NodeMapTable + _side_nodes = derived_side_nodes + (_edge_nodes.values, [](unsigned int s) { return s == 0 ? 18u : s == 4 ? 19u : 14 + s; }); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; virtual void permute(unsigned int perm_num) override final; diff --git a/include/geom/cell_prism6.h b/include/geom/cell_prism6.h index 27cdd174dc..c172c12ce7 100644 --- a/include/geom/cell_prism6.h +++ b/include/geom/cell_prism6.h @@ -169,7 +169,14 @@ class Prism6 final : public Prism * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to * element node numbers. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; + static constexpr unsigned int side_nodes_map[num_sides][nodes_per_side] = + { + {0, 2, 1, 99}, // Side 0 + {0, 1, 4, 3}, // Side 1 + {1, 2, 5, 4}, // Side 2 + {2, 0, 3, 5}, // Side 3 + {3, 4, 5, 99} // Side 4 + }; /** * This maps the child elements with the associated side of the parent element @@ -180,7 +187,18 @@ class Prism6 final : public Prism * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to * element node numbers. */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr unsigned int edge_nodes_map[num_edges][nodes_per_edge] = + { + {0, 1}, // Edge 0 + {1, 2}, // Edge 1 + {0, 2}, // Edge 2 + {0, 3}, // Edge 3 + {1, 4}, // Edge 4 + {2, 5}, // Edge 5 + {3, 4}, // Edge 6 + {4, 5}, // Edge 7 + {3, 5} // Edge 8 + }; /** * An Optimized numerical quadrature approach for computing the diff --git a/include/geom/cell_pyramid13.h b/include/geom/cell_pyramid13.h index f0d2819fb2..8208f1f84c 100644 --- a/include/geom/cell_pyramid13.h +++ b/include/geom/cell_pyramid13.h @@ -22,6 +22,8 @@ // Local includes #include "libmesh/cell_pyramid.h" +#include "libmesh/cell_pyramid5.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -210,16 +212,18 @@ class Pyramid13 final : public Pyramid static const int nodes_per_edge = 3; /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or + * side to element node numbers. They are derived from the + * first-order Pyramid5 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; - - /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to - * element node numbers. - */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr NodeMapTable + _edge_nodes = derived_edge_nodes(); + static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values; + + static constexpr NodeMapTable + _side_nodes = derived_side_nodes + (_edge_nodes.values, no_face_node); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; /** * Specialization for computing the volume of a Pyramid13. diff --git a/include/geom/cell_pyramid14.h b/include/geom/cell_pyramid14.h index 87547dc63e..4329d6a896 100644 --- a/include/geom/cell_pyramid14.h +++ b/include/geom/cell_pyramid14.h @@ -22,6 +22,8 @@ // Local includes #include "libmesh/cell_pyramid.h" +#include "libmesh/cell_pyramid5.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -228,16 +230,18 @@ class Pyramid14 final : public Pyramid static const int nodes_per_edge = 3; /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or + * side to element node numbers. They are derived from the + * first-order Pyramid5 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; - - /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to - * element node numbers. - */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr NodeMapTable + _edge_nodes = derived_edge_nodes(); + static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values; + + static constexpr NodeMapTable + _side_nodes = derived_side_nodes + (_edge_nodes.values, [](unsigned int s) { return s == 4 ? 13 : invalid_uint; }); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; /** * Specialization for computing the volume of a Pyramid14. diff --git a/include/geom/cell_pyramid18.h b/include/geom/cell_pyramid18.h index 33f5c21e70..0df26624e6 100644 --- a/include/geom/cell_pyramid18.h +++ b/include/geom/cell_pyramid18.h @@ -22,6 +22,8 @@ // Local includes #include "libmesh/cell_pyramid.h" +#include "libmesh/cell_pyramid5.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -235,16 +237,18 @@ class Pyramid18 final : public Pyramid static const int nodes_per_edge = 3; /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or + * side to element node numbers. They are derived from the + * first-order Pyramid5 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; - - /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to - * element node numbers. - */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr NodeMapTable + _edge_nodes = derived_edge_nodes(); + static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values; + + static constexpr NodeMapTable + _side_nodes = derived_side_nodes + (_edge_nodes.values, [](unsigned int s) { return s == 4 ? 13 : 14 + s; }); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; virtual void permute(unsigned int perm_num) override final; diff --git a/include/geom/cell_pyramid5.h b/include/geom/cell_pyramid5.h index 835a09b850..4edd008351 100644 --- a/include/geom/cell_pyramid5.h +++ b/include/geom/cell_pyramid5.h @@ -167,13 +167,30 @@ class Pyramid5 final : public Pyramid * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to * element node numbers. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; + static constexpr unsigned int side_nodes_map[num_sides][nodes_per_side] = + { + {0, 1, 4, 99}, // Side 0 + {1, 2, 4, 99}, // Side 1 + {2, 3, 4, 99}, // Side 2 + {3, 0, 4, 99}, // Side 3 + {0, 3, 2, 1} // Side 4 + }; /** * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to * element node numbers. */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr unsigned int edge_nodes_map[num_edges][nodes_per_edge] = + { + {0, 1}, // Edge 0 + {1, 2}, // Edge 1 + {2, 3}, // Edge 2 + {0, 3}, // Edge 3 + {0, 4}, // Edge 4 + {1, 4}, // Edge 5 + {2, 4}, // Edge 6 + {3, 4} // Edge 7 + }; /** * We compute the centroid of the Pyramid by treating it as a diff --git a/include/geom/cell_tet10.h b/include/geom/cell_tet10.h index 5f454fe755..12b5b7c20b 100644 --- a/include/geom/cell_tet10.h +++ b/include/geom/cell_tet10.h @@ -22,6 +22,8 @@ // Local includes #include "libmesh/cell_tet.h" +#include "libmesh/cell_tet4.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -218,16 +220,18 @@ class Tet10 final : public Tet static const int nodes_per_edge = 3; /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or + * side to element node numbers. They are derived from the + * first-order Tet4 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; - - /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to - * element node numbers. - */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr NodeMapTable + _edge_nodes = derived_edge_nodes(); + static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values; + + static constexpr NodeMapTable + _side_nodes = derived_side_nodes + (_edge_nodes.values, no_face_node); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; /** * A specialization for computing the volume of a Tet10. diff --git a/include/geom/cell_tet14.h b/include/geom/cell_tet14.h index 43245751ee..42362fc8f2 100644 --- a/include/geom/cell_tet14.h +++ b/include/geom/cell_tet14.h @@ -22,6 +22,8 @@ // Local includes #include "libmesh/cell_tet.h" +#include "libmesh/cell_tet4.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -223,16 +225,18 @@ class Tet14 final : public Tet static const int nodes_per_edge = 3; /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * These map the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge or + * side to element node numbers. They are derived from the + * first-order Tet4 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; - - /** - * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to - * element node numbers. - */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr NodeMapTable + _edge_nodes = derived_edge_nodes(); + static constexpr const unsigned int (&edge_nodes_map)[num_edges][nodes_per_edge] = _edge_nodes.values; + + static constexpr NodeMapTable + _side_nodes = derived_side_nodes + (_edge_nodes.values, [](unsigned int s) { return 10 + s; }); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; virtual void permute(unsigned int perm_num) override final; diff --git a/include/geom/cell_tet4.h b/include/geom/cell_tet4.h index c6efd46cf3..f31e259a7c 100644 --- a/include/geom/cell_tet4.h +++ b/include/geom/cell_tet4.h @@ -192,13 +192,27 @@ class Tet4 final : public Tet * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to * element node numbers. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; + static constexpr unsigned int side_nodes_map[num_sides][nodes_per_side] = + { + {0, 2, 1}, // Side 0 + {0, 1, 3}, // Side 1 + {1, 2, 3}, // Side 2 + {2, 0, 3} // Side 3 + }; /** * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to * element node numbers. */ - static const unsigned int edge_nodes_map[num_edges][nodes_per_edge]; + static constexpr unsigned int edge_nodes_map[num_edges][nodes_per_edge] = + { + {0, 1}, // Edge 0 + {1, 2}, // Edge 1 + {0, 2}, // Edge 2 + {0, 3}, // Edge 3 + {1, 3}, // Edge 4 + {2, 3} // Edge 5 + }; /** * The centroid of a 4-node tetrahedron is simply given by the diff --git a/include/geom/elem.h b/include/geom/elem.h index 75611008ce..43ec15d5cc 100644 --- a/include/geom/elem.h +++ b/include/geom/elem.h @@ -468,6 +468,27 @@ class Elem : public ReferenceCountedObject, virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const = 0; + /** + * \returns The local node id for node \p side_node on side \p side of + * an element of type \p t, without needing an instantiated Elem. + * The Polygon and Polyhedron subclasses have no such map, and the + * infinite elements' maps are not read here, so those must be queried + * through an actual Elem. + */ + static unsigned int local_side_node(ElemType t, + unsigned int side, + unsigned int side_node); + + /** + * \returns The local node id for node \p edge_node on edge \p edge of + * an element of type \p t, without needing an instantiated Elem. For + * 2D types this is local_side_node(); 1D types have no edges. The + * same types are unsupported here as in local_side_node(). + */ + static unsigned int local_edge_node(ElemType t, + unsigned int edge, + unsigned int edge_node); + /** * \returns \p true if a vertex of \p e is contained * in this element. If \p mesh_connection is true, looks @@ -640,7 +661,67 @@ class Elem : public ReferenceCountedObject, * is fixed; for more general types like Polygon subclasses an actual * instantiated Elem must be queried. */ - static const unsigned int type_to_n_nodes_map[INVALID_ELEM]; + static constexpr unsigned int type_to_n_nodes_map[INVALID_ELEM] = + { + 2, // EDGE2 + 3, // EDGE3 + 4, // EDGE4 + + 3, // TRI3 + 6, // TRI6 + + 4, // QUAD4 + 8, // QUAD8 + 9, // QUAD9 + + 4, // TET4 + 10, // TET10 + + 8, // HEX8 + 20, // HEX20 + 27, // HEX27 + + 6, // PRISM6 + 15, // PRISM15 + 18, // PRISM18 + + 5, // PYRAMID5 + 13, // PYRAMID13 + 14, // PYRAMID14 + + 2, // INFEDGE2 + + 4, // INFQUAD4 + 6, // INFQUAD6 + + 8, // INFHEX8 + 16, // INFHEX16 + 18, // INFHEX18 + + 6, // INFPRISM6 + 12, // INFPRISM12 + + 1, // NODEELEM + + 0, // REMOTEELEM + + 3, // TRI3SUBDIVISION + 3, // TRISHELL3 + 4, // QUADSHELL4 + 8, // QUADSHELL8 + + 7, // TRI7 + 14, // TET14 + 20, // PRISM20 + 21, // PRISM21 + 18, // PYRAMID18 + + 9, // QUADSHELL9 + + invalid_uint, // C0POLYGON + invalid_uint, // C0POLYHEDRON + + }; /** * \returns The number of nodes this element contains. @@ -689,6 +770,24 @@ class Elem : public ReferenceCountedObject, */ virtual ElemType side_type (const unsigned int s) const = 0; + /** + * \returns The type of side \p s of an element of type \p t, without + * needing an instantiated Elem. The Polygon and Polyhedron subclasses + * have one side type but no fixed number of sides, so \p s goes + * unchecked for them; query an actual Elem when you have one. + */ + static ElemType side_type (const ElemType t, + const unsigned int s); + + /** + * \returns The type of every edge of an element of type \p t, or + * \p INVALID_ELEM for the 1D types, which have no edges. Unlike its + * sides, a finite element's edges all have the same type, so no edge + * index is needed; the infinite elements, whose finite and infinite + * edges differ, are not answered here. + */ + static ElemType edge_type (const ElemType t); + /** * \returns the normal (outwards-facing) of the side of the element at the vertex-average of the side * @param s the side of interest diff --git a/include/geom/face_quad4.h b/include/geom/face_quad4.h index 2bc20fa073..0fa0ae6a19 100644 --- a/include/geom/face_quad4.h +++ b/include/geom/face_quad4.h @@ -152,7 +152,13 @@ class Quad4 : public Quad * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to * element node numbers. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; + static constexpr unsigned int side_nodes_map[num_sides][nodes_per_side] = + { + {0, 1}, // Side 0 + {1, 2}, // Side 1 + {2, 3}, // Side 2 + {3, 0} // Side 3 + }; /** * An optimized method for computing the centroid of a diff --git a/include/geom/face_quad8.h b/include/geom/face_quad8.h index 4bcdaa3943..6e5ad55600 100644 --- a/include/geom/face_quad8.h +++ b/include/geom/face_quad8.h @@ -23,6 +23,8 @@ // Local includes #include "libmesh/libmesh_common.h" #include "libmesh/face_quad.h" +#include "libmesh/face_quad4.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -192,9 +194,12 @@ class Quad8 : public Quad /** * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * element node numbers. It is derived from the first-order + * Quad4 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; + static constexpr NodeMapTable + _side_nodes = derived_side_nodes(); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; /** * An optimized method for approximating the area of a diff --git a/include/geom/face_quad9.h b/include/geom/face_quad9.h index e065d8bafa..28b3aff481 100644 --- a/include/geom/face_quad9.h +++ b/include/geom/face_quad9.h @@ -23,6 +23,8 @@ // Local includes #include "libmesh/libmesh_common.h" #include "libmesh/face_quad.h" +#include "libmesh/face_quad4.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -199,9 +201,12 @@ class Quad9 : public Quad /** * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * element node numbers. It is derived from the first-order + * Quad4 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; + static constexpr NodeMapTable + _side_nodes = derived_side_nodes(); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; /** * An optimized method for approximating the area of a diff --git a/include/geom/face_tri3.h b/include/geom/face_tri3.h index 583b80c8ab..71f993a8ee 100644 --- a/include/geom/face_tri3.h +++ b/include/geom/face_tri3.h @@ -166,7 +166,12 @@ class Tri3 : public Tri * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to * element node numbers. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; + static constexpr unsigned int side_nodes_map[num_sides][nodes_per_side] = + { + {0, 1}, // Side 0 + {1, 2}, // Side 1 + {2, 0} // Side 2 + }; /** * The centroid of a 3-node triangle is simply given by the diff --git a/include/geom/face_tri6.h b/include/geom/face_tri6.h index 6417999e9f..38494ce86f 100644 --- a/include/geom/face_tri6.h +++ b/include/geom/face_tri6.h @@ -23,6 +23,8 @@ // Local includes #include "libmesh/libmesh_common.h" #include "libmesh/face_tri.h" +#include "libmesh/face_tri3.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -202,9 +204,12 @@ class Tri6 : public Tri /** * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * element node numbers. It is derived from the first-order + * Tri3 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; + static constexpr NodeMapTable + _side_nodes = derived_side_nodes(); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; /** * An optimized method for approximating the area of a diff --git a/include/geom/face_tri7.h b/include/geom/face_tri7.h index 833bb3ac3b..70fb352277 100644 --- a/include/geom/face_tri7.h +++ b/include/geom/face_tri7.h @@ -23,6 +23,8 @@ // Local includes #include "libmesh/libmesh_common.h" #include "libmesh/face_tri.h" +#include "libmesh/face_tri3.h" +#include "libmesh/fe_reference_element_traits.h" namespace libMesh { @@ -206,9 +208,12 @@ class Tri7 : public Tri /** * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to - * element node numbers. + * element node numbers. It is derived from the first-order + * Tri3 tables; see fe_reference_element_traits.h. */ - static const unsigned int side_nodes_map[num_sides][nodes_per_side]; + static constexpr NodeMapTable + _side_nodes = derived_side_nodes(); + static constexpr const unsigned int (&side_nodes_map)[num_sides][nodes_per_side] = _side_nodes.values; /** * \returns A bounding box (not necessarily the minimal bounding box) diff --git a/include/include_HEADERS b/include/include_HEADERS index a7ddd59ad0..e19ebb70d4 100644 --- a/include/include_HEADERS +++ b/include/include_HEADERS @@ -90,6 +90,7 @@ include_HEADERS = \ fe/fe_lagrange_shape_1D.h \ fe/fe_macro.h \ fe/fe_map.h \ + fe/fe_reference_element_traits.h \ fe/fe_transformation_base.h \ fe/fe_type.h \ fe/fe_xyz_map.h \ diff --git a/include/libmesh/Makefile.am b/include/libmesh/Makefile.am index b3eb689f91..0083c34c37 100644 --- a/include/libmesh/Makefile.am +++ b/include/libmesh/Makefile.am @@ -81,6 +81,7 @@ BUILT_SOURCES = \ fe_lagrange_shape_1D.h \ fe_macro.h \ fe_map.h \ + fe_reference_element_traits.h \ fe_transformation_base.h \ fe_type.h \ fe_xyz_map.h \ @@ -852,6 +853,9 @@ fe_macro.h: $(top_srcdir)/include/fe/fe_macro.h fe_map.h: $(top_srcdir)/include/fe/fe_map.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ +fe_reference_element_traits.h: $(top_srcdir)/include/fe/fe_reference_element_traits.h + $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ + fe_transformation_base.h: $(top_srcdir)/include/fe/fe_transformation_base.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ diff --git a/include/libmesh/Makefile.in b/include/libmesh/Makefile.in index b2f773c1b6..98e72774ce 100644 --- a/include/libmesh/Makefile.in +++ b/include/libmesh/Makefile.in @@ -574,7 +574,8 @@ BUILT_SOURCES = dirichlet_boundaries.h dof_map.h dof_map_base.h \ weighted_patch_recovery_error_estimator.h fe.h fe_abstract.h \ fe_base.h fe_compute_data.h fe_interface.h \ fe_interface_macros.h fe_lagrange_shape_1D.h fe_macro.h \ - fe_map.h fe_transformation_base.h fe_type.h fe_xyz_map.h \ + fe_map.h fe_reference_element_traits.h \ + fe_transformation_base.h fe_type.h fe_xyz_map.h \ h1_fe_transformation.h hcurl_fe_transformation.h \ hdiv_fe_transformation.h inf_fe.h inf_fe_instantiate_1D.h \ inf_fe_instantiate_2D.h inf_fe_instantiate_3D.h inf_fe_macro.h \ @@ -1200,6 +1201,9 @@ fe_macro.h: $(top_srcdir)/include/fe/fe_macro.h fe_map.h: $(top_srcdir)/include/fe/fe_map.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ +fe_reference_element_traits.h: $(top_srcdir)/include/fe/fe_reference_element_traits.h + $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ + fe_transformation_base.h: $(top_srcdir)/include/fe/fe_transformation_base.h $(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@ diff --git a/src/geom/cell_hex20.C b/src/geom/cell_hex20.C index f8351801d7..077d01b5f1 100644 --- a/src/geom/cell_hex20.C +++ b/src/geom/cell_hex20.C @@ -34,32 +34,6 @@ const int Hex20::num_nodes; const int Hex20::nodes_per_side; const int Hex20::nodes_per_edge; -const unsigned int Hex20::side_nodes_map[Hex20::num_sides][Hex20::nodes_per_side] = - { - {0, 3, 2, 1, 11, 10, 9, 8}, // Side 0 - {0, 1, 5, 4, 8, 13, 16, 12}, // Side 1 - {1, 2, 6, 5, 9, 14, 17, 13}, // Side 2 - {2, 3, 7, 6, 10, 15, 18, 14}, // Side 3 - {3, 0, 4, 7, 11, 12, 19, 15}, // Side 4 - {4, 5, 6, 7, 16, 17, 18, 19} // Side 5 - }; - -const unsigned int Hex20::edge_nodes_map[Hex20::num_edges][Hex20::nodes_per_edge] = - { - {0, 1, 8}, // Edge 0 - {1, 2, 9}, // Edge 1 - {2, 3, 10}, // Edge 2 - {0, 3, 11}, // Edge 3 - {0, 4, 12}, // Edge 4 - {1, 5, 13}, // Edge 5 - {2, 6, 14}, // Edge 6 - {3, 7, 15}, // Edge 7 - {4, 5, 16}, // Edge 8 - {5, 6, 17}, // Edge 9 - {6, 7, 18}, // Edge 10 - {4, 7, 19} // Edge 11 - }; - // ------------------------------------------------------------ // Hex20 class member functions diff --git a/src/geom/cell_hex27.C b/src/geom/cell_hex27.C index 9352873889..4877a1fce6 100644 --- a/src/geom/cell_hex27.C +++ b/src/geom/cell_hex27.C @@ -34,32 +34,6 @@ const int Hex27::num_nodes; const int Hex27::nodes_per_side; const int Hex27::nodes_per_edge; -const unsigned int Hex27::side_nodes_map[Hex27::num_sides][Hex27::nodes_per_side] = - { - {0, 3, 2, 1, 11, 10, 9, 8, 20}, // Side 0 - {0, 1, 5, 4, 8, 13, 16, 12, 21}, // Side 1 - {1, 2, 6, 5, 9, 14, 17, 13, 22}, // Side 2 - {2, 3, 7, 6, 10, 15, 18, 14, 23}, // Side 3 - {3, 0, 4, 7, 11, 12, 19, 15, 24}, // Side 4 - {4, 5, 6, 7, 16, 17, 18, 19, 25} // Side 5 - }; - -const unsigned int Hex27::edge_nodes_map[Hex27::num_edges][Hex27::nodes_per_edge] = - { - {0, 1, 8}, // Edge 0 - {1, 2, 9}, // Edge 1 - {2, 3, 10}, // Edge 2 - {0, 3, 11}, // Edge 3 - {0, 4, 12}, // Edge 4 - {1, 5, 13}, // Edge 5 - {2, 6, 14}, // Edge 6 - {3, 7, 15}, // Edge 7 - {4, 5, 16}, // Edge 8 - {5, 6, 17}, // Edge 9 - {6, 7, 18}, // Edge 10 - {4, 7, 19} // Edge 11 - }; - // ------------------------------------------------------------ // Hex27 class member functions diff --git a/src/geom/cell_hex8.C b/src/geom/cell_hex8.C index e3fb308739..f157c7b955 100644 --- a/src/geom/cell_hex8.C +++ b/src/geom/cell_hex8.C @@ -39,32 +39,6 @@ const int Hex8::num_nodes; const int Hex8::nodes_per_side; const int Hex8::nodes_per_edge; -const unsigned int Hex8::side_nodes_map[Hex8::num_sides][Hex8::nodes_per_side] = - { - {0, 3, 2, 1}, // Side 0 - {0, 1, 5, 4}, // Side 1 - {1, 2, 6, 5}, // Side 2 - {2, 3, 7, 6}, // Side 3 - {3, 0, 4, 7}, // Side 4 - {4, 5, 6, 7} // Side 5 - }; - -const unsigned int Hex8::edge_nodes_map[Hex8::num_edges][Hex8::nodes_per_edge] = - { - {0, 1}, // Edge 0 - {1, 2}, // Edge 1 - {2, 3}, // Edge 2 - {0, 3}, // Edge 3 - {0, 4}, // Edge 4 - {1, 5}, // Edge 5 - {2, 6}, // Edge 6 - {3, 7}, // Edge 7 - {4, 5}, // Edge 8 - {5, 6}, // Edge 9 - {6, 7}, // Edge 10 - {4, 7} // Edge 11 - }; - // ------------------------------------------------------------ // Hex8 class member functions diff --git a/src/geom/cell_prism15.C b/src/geom/cell_prism15.C index 1fdeded2e3..a89e280263 100644 --- a/src/geom/cell_prism15.C +++ b/src/geom/cell_prism15.C @@ -35,28 +35,6 @@ const int Prism15::num_nodes; const int Prism15::nodes_per_side; const int Prism15::nodes_per_edge; -const unsigned int Prism15::side_nodes_map[Prism15::num_sides][Prism15::nodes_per_side] = - { - {0, 2, 1, 8, 7, 6, 99, 99}, // Side 0 - {0, 1, 4, 3, 6, 10, 12, 9}, // Side 1 - {1, 2, 5, 4, 7, 11, 13, 10}, // Side 2 - {2, 0, 3, 5, 8, 9, 14, 11}, // Side 3 - {3, 4, 5, 12, 13, 14, 99, 99} // Side 4 - }; - -const unsigned int Prism15::edge_nodes_map[Prism15::num_edges][Prism15::nodes_per_edge] = - { - {0, 1, 6}, // Edge 0 - {1, 2, 7}, // Edge 1 - {0, 2, 8}, // Edge 2 - {0, 3, 9}, // Edge 3 - {1, 4, 10}, // Edge 4 - {2, 5, 11}, // Edge 5 - {3, 4, 12}, // Edge 6 - {4, 5, 13}, // Edge 7 - {3, 5, 14} // Edge 8 - }; - // ------------------------------------------------------------ // Prism15 class member functions diff --git a/src/geom/cell_prism18.C b/src/geom/cell_prism18.C index 8b708d6d05..373ec1e19e 100644 --- a/src/geom/cell_prism18.C +++ b/src/geom/cell_prism18.C @@ -36,28 +36,6 @@ const int Prism18::num_nodes; const int Prism18::nodes_per_side; const int Prism18::nodes_per_edge; -const unsigned int Prism18::side_nodes_map[Prism18::num_sides][Prism18::nodes_per_side] = - { - {0, 2, 1, 8, 7, 6, 99, 99, 99}, // Side 0 - {0, 1, 4, 3, 6, 10, 12, 9, 15}, // Side 1 - {1, 2, 5, 4, 7, 11, 13, 10, 16}, // Side 2 - {2, 0, 3, 5, 8, 9, 14, 11, 17}, // Side 3 - {3, 4, 5, 12, 13, 14, 99, 99, 99} // Side 4 - }; - -const unsigned int Prism18::edge_nodes_map[Prism18::num_edges][Prism18::nodes_per_edge] = - { - {0, 1, 6}, // Edge 0 - {1, 2, 7}, // Edge 1 - {0, 2, 8}, // Edge 2 - {0, 3, 9}, // Edge 3 - {1, 4, 10}, // Edge 4 - {2, 5, 11}, // Edge 5 - {3, 4, 12}, // Edge 6 - {4, 5, 13}, // Edge 7 - {3, 5, 14} // Edge 8 - }; - // ------------------------------------------------------------ // Prism18 class member functions diff --git a/src/geom/cell_prism20.C b/src/geom/cell_prism20.C index c1cab40856..8cd5de7306 100644 --- a/src/geom/cell_prism20.C +++ b/src/geom/cell_prism20.C @@ -36,28 +36,6 @@ const int Prism20::num_nodes; const int Prism20::nodes_per_side; const int Prism20::nodes_per_edge; -const unsigned int Prism20::side_nodes_map[Prism20::num_sides][Prism20::nodes_per_side] = - { - {0, 2, 1, 8, 7, 6, 18, 99, 99}, // Side 0 - {0, 1, 4, 3, 6, 10, 12, 9, 15}, // Side 1 - {1, 2, 5, 4, 7, 11, 13, 10, 16}, // Side 2 - {2, 0, 3, 5, 8, 9, 14, 11, 17}, // Side 3 - {3, 4, 5, 12, 13, 14, 19, 99, 99} // Side 4 - }; - -const unsigned int Prism20::edge_nodes_map[Prism20::num_edges][Prism20::nodes_per_edge] = - { - {0, 1, 6}, // Edge 0 - {1, 2, 7}, // Edge 1 - {0, 2, 8}, // Edge 2 - {0, 3, 9}, // Edge 3 - {1, 4, 10}, // Edge 4 - {2, 5, 11}, // Edge 5 - {3, 4, 12}, // Edge 6 - {4, 5, 13}, // Edge 7 - {3, 5, 14} // Edge 8 - }; - // ------------------------------------------------------------ // Prism20 class member functions diff --git a/src/geom/cell_prism21.C b/src/geom/cell_prism21.C index 2e6a577784..807f74cf40 100644 --- a/src/geom/cell_prism21.C +++ b/src/geom/cell_prism21.C @@ -51,28 +51,6 @@ const int Prism21::num_nodes; const int Prism21::nodes_per_side; const int Prism21::nodes_per_edge; -const unsigned int Prism21::side_nodes_map[Prism21::num_sides][Prism21::nodes_per_side] = - { - {0, 2, 1, 8, 7, 6, 18, 99, 99}, // Side 0 - {0, 1, 4, 3, 6, 10, 12, 9, 15}, // Side 1 - {1, 2, 5, 4, 7, 11, 13, 10, 16}, // Side 2 - {2, 0, 3, 5, 8, 9, 14, 11, 17}, // Side 3 - {3, 4, 5, 12, 13, 14, 19, 99, 99} // Side 4 - }; - -const unsigned int Prism21::edge_nodes_map[Prism21::num_edges][Prism21::nodes_per_edge] = - { - {0, 1, 6}, // Edge 0 - {1, 2, 7}, // Edge 1 - {0, 2, 8}, // Edge 2 - {0, 3, 9}, // Edge 3 - {1, 4, 10}, // Edge 4 - {2, 5, 11}, // Edge 5 - {3, 4, 12}, // Edge 6 - {4, 5, 13}, // Edge 7 - {3, 5, 14} // Edge 8 - }; - // ------------------------------------------------------------ // Prism21 class member functions diff --git a/src/geom/cell_prism6.C b/src/geom/cell_prism6.C index a520b137cc..6a972dc5c4 100644 --- a/src/geom/cell_prism6.C +++ b/src/geom/cell_prism6.C @@ -110,15 +110,6 @@ const int Prism6::num_nodes; const int Prism6::nodes_per_side; const int Prism6::nodes_per_edge; -const unsigned int Prism6::side_nodes_map[Prism6::num_sides][Prism6::nodes_per_side] = - { - {0, 2, 1, 99}, // Side 0 - {0, 1, 4, 3}, // Side 1 - {1, 2, 5, 4}, // Side 2 - {2, 0, 3, 5}, // Side 3 - {3, 4, 5, 99} // Side 4 - }; - const unsigned int Prism6::side_elems_map[Prism6::num_sides][Prism6::nodes_per_side] = { {0, 1, 2, 3}, // Side 0 @@ -128,19 +119,6 @@ const unsigned int Prism6::side_elems_map[Prism6::num_sides][Prism6::nodes_per_s {4, 5, 6, 7} // Side 4 }; -const unsigned int Prism6::edge_nodes_map[Prism6::num_edges][Prism6::nodes_per_edge] = - { - {0, 1}, // Edge 0 - {1, 2}, // Edge 1 - {0, 2}, // Edge 2 - {0, 3}, // Edge 3 - {1, 4}, // Edge 4 - {2, 5}, // Edge 5 - {3, 4}, // Edge 6 - {4, 5}, // Edge 7 - {3, 5} // Edge 8 - }; - // ------------------------------------------------------------ // Prism6 class member functions diff --git a/src/geom/cell_pyramid13.C b/src/geom/cell_pyramid13.C index 2f1e878af0..068252f088 100644 --- a/src/geom/cell_pyramid13.C +++ b/src/geom/cell_pyramid13.C @@ -36,27 +36,6 @@ const int Pyramid13::num_nodes; const int Pyramid13::nodes_per_side; const int Pyramid13::nodes_per_edge; -const unsigned int Pyramid13::side_nodes_map[Pyramid13::num_sides][Pyramid13::nodes_per_side] = - { - {0, 1, 4, 5, 10, 9, 99, 99}, // Side 0 (front) - {1, 2, 4, 6, 11, 10, 99, 99}, // Side 1 (right) - {2, 3, 4, 7, 12, 11, 99, 99}, // Side 2 (back) - {3, 0, 4, 8, 9, 12, 99, 99}, // Side 3 (left) - {0, 3, 2, 1, 8, 7, 6, 5} // Side 4 (base) - }; - -const unsigned int Pyramid13::edge_nodes_map[Pyramid13::num_edges][Pyramid13::nodes_per_edge] = - { - {0, 1, 5}, // Edge 0 - {1, 2, 6}, // Edge 1 - {2, 3, 7}, // Edge 2 - {0, 3, 8}, // Edge 3 - {0, 4, 9}, // Edge 4 - {1, 4, 10}, // Edge 5 - {2, 4, 11}, // Edge 6 - {3, 4, 12} // Edge 7 - }; - // ------------------------------------------------------------ // Pyramid13 class member functions diff --git a/src/geom/cell_pyramid14.C b/src/geom/cell_pyramid14.C index 73b056e0e2..769ca89862 100644 --- a/src/geom/cell_pyramid14.C +++ b/src/geom/cell_pyramid14.C @@ -36,27 +36,6 @@ const int Pyramid14::num_nodes; const int Pyramid14::nodes_per_side; const int Pyramid14::nodes_per_edge; -const unsigned int Pyramid14::side_nodes_map[Pyramid14::num_sides][Pyramid14::nodes_per_side] = - { - {0, 1, 4, 5, 10, 9, 99, 99, 99}, // Side 0 (front) - {1, 2, 4, 6, 11, 10, 99, 99, 99}, // Side 1 (right) - {2, 3, 4, 7, 12, 11, 99, 99, 99}, // Side 2 (back) - {3, 0, 4, 8, 9, 12, 99, 99, 99}, // Side 3 (left) - {0, 3, 2, 1, 8, 7, 6, 5, 13} // Side 4 (base) - }; - -const unsigned int Pyramid14::edge_nodes_map[Pyramid14::num_edges][Pyramid14::nodes_per_edge] = - { - {0, 1, 5}, // Edge 0 - {1, 2, 6}, // Edge 1 - {2, 3, 7}, // Edge 2 - {0, 3, 8}, // Edge 3 - {0, 4, 9}, // Edge 4 - {1, 4, 10}, // Edge 5 - {2, 4, 11}, // Edge 6 - {3, 4, 12} // Edge 7 - }; - // ------------------------------------------------------------ // Pyramid14 class member functions diff --git a/src/geom/cell_pyramid18.C b/src/geom/cell_pyramid18.C index 12f7ad69f5..5a4d89af1e 100644 --- a/src/geom/cell_pyramid18.C +++ b/src/geom/cell_pyramid18.C @@ -36,27 +36,6 @@ const int Pyramid18::num_nodes; const int Pyramid18::nodes_per_side; const int Pyramid18::nodes_per_edge; -const unsigned int Pyramid18::side_nodes_map[Pyramid18::num_sides][Pyramid18::nodes_per_side] = - { - {0, 1, 4, 5, 10, 9, 14, 99, 99}, // Side 0 (front) - {1, 2, 4, 6, 11, 10, 15, 99, 99}, // Side 1 (right) - {2, 3, 4, 7, 12, 11, 16, 99, 99}, // Side 2 (back) - {3, 0, 4, 8, 9, 12, 17, 99, 99}, // Side 3 (left) - {0, 3, 2, 1, 8, 7, 6, 5, 13} // Side 4 (base) - }; - -const unsigned int Pyramid18::edge_nodes_map[Pyramid18::num_edges][Pyramid18::nodes_per_edge] = - { - {0, 1, 5}, // Edge 0 - {1, 2, 6}, // Edge 1 - {2, 3, 7}, // Edge 2 - {0, 3, 8}, // Edge 3 - {0, 4, 9}, // Edge 4 - {1, 4, 10}, // Edge 5 - {2, 4, 11}, // Edge 6 - {3, 4, 12} // Edge 7 - }; - // ------------------------------------------------------------ // Pyramid18 class member functions diff --git a/src/geom/cell_pyramid5.C b/src/geom/cell_pyramid5.C index eb539aa591..c8c997ba95 100644 --- a/src/geom/cell_pyramid5.C +++ b/src/geom/cell_pyramid5.C @@ -37,27 +37,6 @@ const int Pyramid5::num_nodes; const int Pyramid5::nodes_per_side; const int Pyramid5::nodes_per_edge; -const unsigned int Pyramid5::side_nodes_map[Pyramid5::num_sides][Pyramid5::nodes_per_side] = - { - {0, 1, 4, 99}, // Side 0 - {1, 2, 4, 99}, // Side 1 - {2, 3, 4, 99}, // Side 2 - {3, 0, 4, 99}, // Side 3 - {0, 3, 2, 1} // Side 4 - }; - -const unsigned int Pyramid5::edge_nodes_map[Pyramid5::num_edges][Pyramid5::nodes_per_edge] = - { - {0, 1}, // Edge 0 - {1, 2}, // Edge 1 - {2, 3}, // Edge 2 - {0, 3}, // Edge 3 - {0, 4}, // Edge 4 - {1, 4}, // Edge 5 - {2, 4}, // Edge 6 - {3, 4} // Edge 7 - }; - // ------------------------------------------------------------ // Pyramid5 class member functions diff --git a/src/geom/cell_tet10.C b/src/geom/cell_tet10.C index 511eba251c..9896a24e41 100644 --- a/src/geom/cell_tet10.C +++ b/src/geom/cell_tet10.C @@ -34,24 +34,6 @@ const int Tet10::num_nodes; const int Tet10::nodes_per_side; const int Tet10::nodes_per_edge; -const unsigned int Tet10::side_nodes_map[Tet10::num_sides][Tet10::nodes_per_side] = - { - {0, 2, 1, 6, 5, 4}, // Side 0 - {0, 1, 3, 4, 8, 7}, // Side 1 - {1, 2, 3, 5, 9, 8}, // Side 2 - {2, 0, 3, 6, 7, 9} // Side 3 - }; - -const unsigned int Tet10::edge_nodes_map[Tet10::num_edges][Tet10::nodes_per_edge] = - { - {0, 1, 4}, // Edge 0 - {1, 2, 5}, // Edge 1 - {0, 2, 6}, // Edge 2 - {0, 3, 7}, // Edge 3 - {1, 3, 8}, // Edge 4 - {2, 3, 9} // Edge 5 - }; - // ------------------------------------------------------------ // Tet10 class member functions diff --git a/src/geom/cell_tet14.C b/src/geom/cell_tet14.C index b214ee1c36..91d3a9c653 100644 --- a/src/geom/cell_tet14.C +++ b/src/geom/cell_tet14.C @@ -42,24 +42,6 @@ const int Tet14::num_nodes; const int Tet14::nodes_per_side; const int Tet14::nodes_per_edge; -const unsigned int Tet14::side_nodes_map[Tet14::num_sides][Tet14::nodes_per_side] = - { - {0, 2, 1, 6, 5, 4, 10}, // Side 0 - {0, 1, 3, 4, 8, 7, 11}, // Side 1 - {1, 2, 3, 5, 9, 8, 12}, // Side 2 - {2, 0, 3, 6, 7, 9, 13} // Side 3 - }; - -const unsigned int Tet14::edge_nodes_map[Tet14::num_edges][Tet14::nodes_per_edge] = - { - {0, 1, 4}, // Edge 0 - {1, 2, 5}, // Edge 1 - {0, 2, 6}, // Edge 2 - {0, 3, 7}, // Edge 3 - {1, 3, 8}, // Edge 4 - {2, 3, 9} // Edge 5 - }; - // ------------------------------------------------------------ // Tet14 class member functions diff --git a/src/geom/cell_tet4.C b/src/geom/cell_tet4.C index fafde46173..63b8216efa 100644 --- a/src/geom/cell_tet4.C +++ b/src/geom/cell_tet4.C @@ -35,24 +35,6 @@ const int Tet4::num_nodes; const int Tet4::nodes_per_side; const int Tet4::nodes_per_edge; -const unsigned int Tet4::side_nodes_map[Tet4::num_sides][Tet4::nodes_per_side] = - { - {0, 2, 1}, // Side 0 - {0, 1, 3}, // Side 1 - {1, 2, 3}, // Side 2 - {2, 0, 3} // Side 3 - }; - -const unsigned int Tet4::edge_nodes_map[Tet4::num_edges][Tet4::nodes_per_edge] = - { - {0, 1}, // Edge 0 - {1, 2}, // Edge 1 - {0, 2}, // Edge 2 - {0, 3}, // Edge 3 - {1, 3}, // Edge 4 - {2, 3} // Edge 5 - }; - // ------------------------------------------------------------ // Tet4 class member functions diff --git a/src/geom/elem.C b/src/geom/elem.C index 9440ddd930..85a4896549 100644 --- a/src/geom/elem.C +++ b/src/geom/elem.C @@ -161,67 +161,6 @@ const unsigned int Elem::type_to_dim_map [] = const unsigned int Elem::max_n_nodes; -const unsigned int Elem::type_to_n_nodes_map [] = - { - 2, // EDGE2 - 3, // EDGE3 - 4, // EDGE4 - - 3, // TRI3 - 6, // TRI6 - - 4, // QUAD4 - 8, // QUAD8 - 9, // QUAD9 - - 4, // TET4 - 10, // TET10 - - 8, // HEX8 - 20, // HEX20 - 27, // HEX27 - - 6, // PRISM6 - 15, // PRISM15 - 18, // PRISM18 - - 5, // PYRAMID5 - 13, // PYRAMID13 - 14, // PYRAMID14 - - 2, // INFEDGE2 - - 4, // INFQUAD4 - 6, // INFQUAD6 - - 8, // INFHEX8 - 16, // INFHEX16 - 18, // INFHEX18 - - 6, // INFPRISM6 - 12, // INFPRISM12 - - 1, // NODEELEM - - 0, // REMOTEELEM - - 3, // TRI3SUBDIVISION - 3, // TRISHELL3 - 4, // QUADSHELL4 - 8, // QUADSHELL8 - - 7, // TRI7 - 14, // TET14 - 20, // PRISM20 - 21, // PRISM21 - 18, // PYRAMID18 - - 9, // QUADSHELL9 - - invalid_uint, // C0POLYGON - invalid_uint, // C0POLYHEDRON - }; - const unsigned int Elem::type_to_n_sides_map [] = { 2, // EDGE2 @@ -3164,6 +3103,307 @@ ElemType Elem::first_order_equivalent_type (const ElemType et) } +// Most of our elements have the same topology on every side, but the +// prisms and pyramids have both triangular and quadrilateral faces, and +// a 2D or 3D infinite element's side 0 is its finite base while the rest +// are infinite. The polygons and polyhedra have one side type but no +// fixed side count, so they can answer here too while their side index +// goes unchecked; ask an actual Elem when you have one. +ElemType Elem::side_type (const ElemType t, + const unsigned int s) +{ + libmesh_assert_less (s, type_to_n_sides_map[t]); + + switch (t) + { + case EDGE2: + case EDGE3: + case EDGE4: + return NODEELEM; + case TRI3: + case TRISHELL3: + case QUAD4: + case QUADSHELL4: + // A polygon's sides are all edges and a polyhedron's are all + // polygons, however many a given element turns out to have + case C0POLYGON: + return EDGE2; + case TRI6: + case TRI7: + case QUAD8: + case QUADSHELL8: + case QUAD9: + case QUADSHELL9: + return EDGE3; + case TET4: + return TRI3; + case TET10: + return TRI6; + case TET14: + return TRI7; + case HEX8: + return QUAD4; + case HEX20: + return QUAD8; + case HEX27: + return QUAD9; + // Sides 0 and 4 of a prism are its triangles + case PRISM6: + return (s == 0 || s == 4) ? TRI3 : QUAD4; + case PRISM15: + return (s == 0 || s == 4) ? TRI6 : QUAD8; + case PRISM18: + return (s == 0 || s == 4) ? TRI6 : QUAD9; + case PRISM20: + case PRISM21: + return (s == 0 || s == 4) ? TRI7 : QUAD9; + // Side 4 of a pyramid is its quadrilateral base + case PYRAMID5: + return (s < 4) ? TRI3 : QUAD4; + case PYRAMID13: + return (s < 4) ? TRI6 : QUAD8; + case PYRAMID14: + return (s < 4) ? TRI6 : QUAD9; + case PYRAMID18: + return (s < 4) ? TRI7 : QUAD9; + case C0POLYHEDRON: + return C0POLYGON; +#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS + // An InfEdge2's sides are its end nodes, like any other 1D element's + case INFEDGE2: + return NODEELEM; + // A 2D or 3D infinite element's side 0 is the finite base it was + // built from; its remaining sides run out to infinity with it + case INFQUAD4: + return (s == 0) ? EDGE2 : INFEDGE2; + case INFQUAD6: + return (s == 0) ? EDGE3 : INFEDGE2; + case INFHEX8: + return (s == 0) ? QUAD4 : INFQUAD4; + case INFHEX16: + return (s == 0) ? QUAD8 : INFQUAD6; + case INFHEX18: + return (s == 0) ? QUAD9 : INFQUAD6; + case INFPRISM6: + return (s == 0) ? TRI3 : INFQUAD4; + case INFPRISM12: + return (s == 0) ? TRI6 : INFQUAD6; +#endif + default: + libmesh_error_msg("No side type for element type " << Utility::enum_to_string(t)); + } + + return INVALID_ELEM; +} + + + +// This reads the same constexpr side_nodes_map tables the element +// classes use for their virtual local_side_node(), so that code which +// knows an element type but has no Elem to call a virtual function on -- +// building a side of an element that doesn't exist yet, or a device +// kernel that has only the type -- can still get at the reference +// element topology. +// +// The polygons and polyhedra have no such table to read, and the +// infinite elements' maps have no caller here yet, so both error out: +// asking for a topology we can't answer is a programming error rather +// than something to signal with a return value. +// A finite element's edges all have the same type, unlike its sides, +// which differ on the prisms and pyramids, so this takes no edge index. +// The infinite elements are the exception -- the edges of their finite +// base are not the ones running out to infinity -- and like the rest of +// the static topology lookups they are left to their virtual overrides. +ElemType Elem::edge_type (const ElemType t) +{ + switch (t) + { + // 1D elements have no edges + case EDGE2: + case EDGE3: + case EDGE4: + return INVALID_ELEM; + // A 2D element's edges are its sides + case TRI3: + case TRISHELL3: + case TRI6: + case TRI7: + case QUAD4: + case QUADSHELL4: + case QUAD8: + case QUADSHELL8: + case QUAD9: + case QUADSHELL9: + case C0POLYGON: + return side_type(t, 0); + // A first-order 3D element's edges hold their two vertices + case TET4: + case HEX8: + case PRISM6: + case PYRAMID5: + case C0POLYHEDRON: + return EDGE2; + // and a second-order element's add the midpoint + case TET10: + case TET14: + case HEX20: + case HEX27: + case PRISM15: + case PRISM18: + case PRISM20: + case PRISM21: + case PYRAMID13: + case PYRAMID14: + case PYRAMID18: + return EDGE3; + default: + libmesh_error_msg("No edge type for element type " << Utility::enum_to_string(t)); + } + + return INVALID_ELEM; +} + + + +unsigned int Elem::local_side_node (const ElemType t, + const unsigned int side, + const unsigned int side_node) +{ + // The nodes on a side are the nodes of the side's own element type, so + // we compose the two lookups here rather than tabulating side node + // counts a second time. That also keeps us inside the meaningful part + // of a row: the prisms' and pyramids' triangular sides carry fewer + // nodes than their quadrilateral ones, and their rows are padded out + // to the longer length + libmesh_assert_less (side, type_to_n_sides_map[t]); + libmesh_assert_less (side_node, type_to_n_nodes_map[side_type(t, side)]); + + switch (t) + { + // A 1D element's sides are its end nodes + case EDGE2: + case EDGE3: + case EDGE4: + return side; + // The shell elements are numbered like the elements they shadow + case TRI3: + case TRISHELL3: + return Tri3::side_nodes_map[side][side_node]; + case TRI6: + return Tri6::side_nodes_map[side][side_node]; + case TRI7: + return Tri7::side_nodes_map[side][side_node]; + case QUAD4: + case QUADSHELL4: + return Quad4::side_nodes_map[side][side_node]; + case QUAD8: + case QUADSHELL8: + return Quad8::side_nodes_map[side][side_node]; + case QUAD9: + case QUADSHELL9: + return Quad9::side_nodes_map[side][side_node]; + case TET4: + return Tet4::side_nodes_map[side][side_node]; + case TET10: + return Tet10::side_nodes_map[side][side_node]; + case TET14: + return Tet14::side_nodes_map[side][side_node]; + case HEX8: + return Hex8::side_nodes_map[side][side_node]; + case HEX20: + return Hex20::side_nodes_map[side][side_node]; + case HEX27: + return Hex27::side_nodes_map[side][side_node]; + case PRISM6: + return Prism6::side_nodes_map[side][side_node]; + case PRISM15: + return Prism15::side_nodes_map[side][side_node]; + case PRISM18: + return Prism18::side_nodes_map[side][side_node]; + case PRISM20: + return Prism20::side_nodes_map[side][side_node]; + case PRISM21: + return Prism21::side_nodes_map[side][side_node]; + case PYRAMID5: + return Pyramid5::side_nodes_map[side][side_node]; + case PYRAMID13: + return Pyramid13::side_nodes_map[side][side_node]; + case PYRAMID14: + return Pyramid14::side_nodes_map[side][side_node]; + case PYRAMID18: + return Pyramid18::side_nodes_map[side][side_node]; + default: + libmesh_error_msg("No static side node map for element type " << Utility::enum_to_string(t)); + } + + return invalid_uint; +} + + + +// A 2D element's edges are its sides -- Face::local_edge_node() defines +// the two to be the same thing -- so we answer those from the side map +// and keep one copy of the 2D numbering. 1D elements have no edges to +// ask about, and like local_side_node() above we error rather than +// return a flag for a type we have no table for. +unsigned int Elem::local_edge_node (const ElemType t, + const unsigned int edge, + const unsigned int edge_node) +{ + libmesh_assert_less (edge, type_to_n_edges_map[t]); + libmesh_assert_less (edge_node, type_to_n_nodes_map[edge_type(t)]); + + switch (t) + { + case TRI3: + case TRISHELL3: + case TRI6: + case TRI7: + case QUAD4: + case QUADSHELL4: + case QUAD8: + case QUADSHELL8: + case QUAD9: + case QUADSHELL9: + return local_side_node(t, edge, edge_node); + case TET4: + return Tet4::edge_nodes_map[edge][edge_node]; + case TET10: + return Tet10::edge_nodes_map[edge][edge_node]; + case TET14: + return Tet14::edge_nodes_map[edge][edge_node]; + case HEX8: + return Hex8::edge_nodes_map[edge][edge_node]; + case HEX20: + return Hex20::edge_nodes_map[edge][edge_node]; + case HEX27: + return Hex27::edge_nodes_map[edge][edge_node]; + case PRISM6: + return Prism6::edge_nodes_map[edge][edge_node]; + case PRISM15: + return Prism15::edge_nodes_map[edge][edge_node]; + case PRISM18: + return Prism18::edge_nodes_map[edge][edge_node]; + case PRISM20: + return Prism20::edge_nodes_map[edge][edge_node]; + case PRISM21: + return Prism21::edge_nodes_map[edge][edge_node]; + case PYRAMID5: + return Pyramid5::edge_nodes_map[edge][edge_node]; + case PYRAMID13: + return Pyramid13::edge_nodes_map[edge][edge_node]; + case PYRAMID14: + return Pyramid14::edge_nodes_map[edge][edge_node]; + case PYRAMID18: + return Pyramid18::edge_nodes_map[edge][edge_node]; + default: + libmesh_error_msg("No static edge node map for element type " << Utility::enum_to_string(t)); + } + + return invalid_uint; +} + + ElemType Elem::second_order_equivalent_type (const ElemType et, const bool full_ordered) diff --git a/src/geom/face_quad4.C b/src/geom/face_quad4.C index 4a1d2356f8..5cf7955b33 100644 --- a/src/geom/face_quad4.C +++ b/src/geom/face_quad4.C @@ -32,14 +32,6 @@ namespace libMesh const int Quad4::num_nodes; const int Quad4::nodes_per_side; -const unsigned int Quad4::side_nodes_map[Quad4::num_sides][Quad4::nodes_per_side] = - { - {0, 1}, // Side 0 - {1, 2}, // Side 1 - {2, 3}, // Side 2 - {3, 0} // Side 3 - }; - #ifdef LIBMESH_ENABLE_AMR const Real Quad4::_embedding_matrix[Quad4::num_children][Quad4::num_nodes][Quad4::num_nodes] = diff --git a/src/geom/face_quad8.C b/src/geom/face_quad8.C index 720f77dd4d..7a7cbfd2a8 100644 --- a/src/geom/face_quad8.C +++ b/src/geom/face_quad8.C @@ -32,14 +32,6 @@ namespace libMesh const int Quad8::num_nodes; const int Quad8::nodes_per_side; -const unsigned int Quad8::side_nodes_map[Quad8::num_sides][Quad8::nodes_per_side] = - { - {0, 1, 4}, // Side 0 - {1, 2, 5}, // Side 1 - {2, 3, 6}, // Side 2 - {3, 0, 7} // Side 3 - }; - #ifdef LIBMESH_ENABLE_AMR diff --git a/src/geom/face_quad9.C b/src/geom/face_quad9.C index 7182b023a4..ddc153d832 100644 --- a/src/geom/face_quad9.C +++ b/src/geom/face_quad9.C @@ -32,14 +32,6 @@ namespace libMesh const int Quad9::num_nodes; const int Quad9::nodes_per_side; -const unsigned int Quad9::side_nodes_map[Quad9::num_sides][Quad9::nodes_per_side] = - { - {0, 1, 4}, // Side 0 - {1, 2, 5}, // Side 1 - {2, 3, 6}, // Side 2 - {3, 0, 7} // Side 3 - }; - #ifdef LIBMESH_ENABLE_AMR diff --git a/src/geom/face_tri3.C b/src/geom/face_tri3.C index 39909e6214..fc7faa3440 100644 --- a/src/geom/face_tri3.C +++ b/src/geom/face_tri3.C @@ -31,13 +31,6 @@ namespace libMesh const int Tri3::num_nodes; const int Tri3::nodes_per_side; -const unsigned int Tri3::side_nodes_map[Tri3::num_sides][Tri3::nodes_per_side] = - { - {0, 1}, // Side 0 - {1, 2}, // Side 1 - {2, 0} // Side 2 - }; - #ifdef LIBMESH_ENABLE_AMR const Real Tri3::_embedding_matrix[Tri3::num_children][Tri3::num_nodes][Tri3::num_nodes] = diff --git a/src/geom/face_tri6.C b/src/geom/face_tri6.C index fd035d7e90..9985859137 100644 --- a/src/geom/face_tri6.C +++ b/src/geom/face_tri6.C @@ -32,13 +32,6 @@ namespace libMesh const int Tri6::num_nodes; const int Tri6::nodes_per_side; -const unsigned int Tri6::side_nodes_map[Tri6::num_sides][Tri6::nodes_per_side] = - { - {0, 1, 3}, // Side 0 - {1, 2, 4}, // Side 1 - {2, 0, 5} // Side 2 - }; - #ifdef LIBMESH_ENABLE_AMR diff --git a/src/geom/face_tri7.C b/src/geom/face_tri7.C index e30c72ced5..c150153787 100644 --- a/src/geom/face_tri7.C +++ b/src/geom/face_tri7.C @@ -38,13 +38,6 @@ namespace libMesh const int Tri7::num_nodes; const int Tri7::nodes_per_side; -const unsigned int Tri7::side_nodes_map[Tri7::num_sides][Tri7::nodes_per_side] = - { - {0, 1, 3}, // Side 0 - {1, 2, 4}, // Side 1 - {2, 0, 5} // Side 2 - }; - #ifdef LIBMESH_ENABLE_AMR diff --git a/tests/geom/elem_test.C b/tests/geom/elem_test.C index 943e4f8554..fc3ceeef70 100644 --- a/tests/geom/elem_test.C +++ b/tests/geom/elem_test.C @@ -75,6 +75,138 @@ public: } } + // The Elem::side_type/local_side_node/local_edge_node overloads that + // take an ElemType have to agree with the virtual versions they shadow; + // an element of each type is the only honest way to check that. + void test_static_topology() + { + LOG_UNIT_TEST; + + for (const auto & elem : this->_mesh->active_local_element_ptr_range()) + { + const ElemType type = elem->type(); + + for (const auto s : elem->side_index_range()) + { + const ElemType side_type = Elem::side_type(type, s); + CPPUNIT_ASSERT_EQUAL(elem->side_type(s), side_type); + + // A polytope's side type follows from its element type even + // though its side count does not, but it has no static node + // map to check; neither are the infinite elements' maps read + // by these lookups, so both stop at the side type. + if (elem->runtime_topology() || elem->infinite()) + continue; + + const auto nodes = elem->nodes_on_side(s); + CPPUNIT_ASSERT_EQUAL(std::size_t(Elem::type_to_n_nodes_map[side_type]), + nodes.size()); + for (auto n : index_range(nodes)) + { + CPPUNIT_ASSERT_EQUAL(elem->local_side_node(s, n), + Elem::local_side_node(type, s, n)); + CPPUNIT_ASSERT_EQUAL(nodes[n], + Elem::local_side_node(type, s, n)); + } + } + + // 1D elements have no edges. A 2D element's edges are its sides, + // but nodes_on_edge() and local_edge_node() are not the APIs the + // loop above exercised, so those are checked here too. + if (elem->infinite() || elem->dim() < 2) + continue; + + for (const auto e : elem->edge_index_range()) + { + // There is no virtual edge_type() to check against, but the + // element classes state the same fact when they build an + // edge, so compare with that instead + CPPUNIT_ASSERT_EQUAL(elem->build_edge_ptr(e)->type(), + Elem::edge_type(type)); + + if (elem->runtime_topology()) + continue; + + const auto nodes = elem->nodes_on_edge(e); + CPPUNIT_ASSERT_EQUAL(std::size_t(Elem::type_to_n_nodes_map[Elem::edge_type(type)]), + nodes.size()); + for (auto n : index_range(nodes)) + { + CPPUNIT_ASSERT_EQUAL(elem->local_edge_node(e, n), + Elem::local_edge_node(type, e, n)); + CPPUNIT_ASSERT_EQUAL(nodes[n], + Elem::local_edge_node(type, e, n)); + } + } + } + } + + void test_higher_order_node_placement() + { + LOG_UNIT_TEST; + + // The side/edge map consistency is a static_assert in each element + // header now; what's left to check at runtime is where the + // higher-order nodes sit. Every non-vertex reference node is at + // the centroid of the vertices of its edge, its face, or the whole + // element, except for EDGE4's nodes, which trisect it. + for (const auto & elem : this->_mesh->active_local_element_ptr_range()) + { + if (elem->infinite() || elem->runtime_topology()) + continue; + + const ElemType type = elem->type(); + + for (const auto i : elem->node_index_range()) + { + if (elem->is_vertex(i)) + continue; + + // EDGE4's interior nodes trisect it instead + if (type == EDGE4) + { + LIBMESH_ASSERT_REALVEC_EQUAL(Point(i == 2 ? Real(-1)/3 : Real(1)/3), + elem->master_point(i), + TOLERANCE*TOLERANCE); + continue; + } + + // Find the smallest subentity the node belongs to: an edge + // if it sits on one, else a face, else the element itself. + // Its vertices are what the node should be the centroid of. + std::vector subentity_vertices; + if (elem->dim() > 1) + for (const auto e : elem->edge_index_range()) + { + const auto nodes = elem->nodes_on_edge(e); + if (std::find(nodes.begin(), nodes.end(), i) != nodes.end()) + subentity_vertices = {nodes[0], nodes[1]}; + } + if (subentity_vertices.empty()) + for (const auto s : elem->side_index_range()) + { + const auto nodes = elem->nodes_on_side(s); + if (std::find(nodes.begin(), nodes.end(), i) != nodes.end()) + for (auto n : nodes) + if (elem->is_vertex(n)) + subentity_vertices.push_back(n); + } + if (subentity_vertices.empty()) + for (const auto n : elem->node_index_range()) + if (elem->is_vertex(n)) + subentity_vertices.push_back(n); + + Point centroid; + for (auto v : subentity_vertices) + centroid += elem->master_point(v); + centroid /= Real(subentity_vertices.size()); + + LIBMESH_ASSERT_REALVEC_EQUAL(centroid, elem->master_point(i), + TOLERANCE*TOLERANCE); + } + } + } + void test_quality() { LOG_UNIT_TEST; @@ -970,6 +1102,8 @@ public: #define ELEMTEST \ CPPUNIT_TEST( test_bounding_box ); \ CPPUNIT_TEST( test_ref_elem ); \ + CPPUNIT_TEST( test_static_topology ); \ + CPPUNIT_TEST( test_higher_order_node_placement ); \ CPPUNIT_TEST( test_quality ); \ CPPUNIT_TEST( test_node_edge_map_consistency ); \ CPPUNIT_TEST( test_maps ); \