From 1e89e7d6a9ed832b27cc910c9a07204bea8bcbed Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 15 Jan 2021 09:07:33 -0700 Subject: [PATCH] Move several free functions into class ReferenceCell::Type. --- include/deal.II/grid/connectivity.h | 5 +- include/deal.II/grid/reference_cell.h | 196 ++++++++++++-------- source/base/qprojector.cc | 12 +- source/fe/fe_nothing.cc | 3 +- source/fe/mapping_fe.cc | 17 +- source/fe/mapping_fe_field.cc | 8 +- source/grid/manifold.cc | 9 +- source/grid/reference_cell.cc | 10 +- tests/simplex/unit_tangential_vectors_01.cc | 4 +- 9 files changed, 147 insertions(+), 117 deletions(-) diff --git a/include/deal.II/grid/connectivity.h b/include/deal.II/grid/connectivity.h index 07b1fc376a..d049fa5f09 100644 --- a/include/deal.II/grid/connectivity.h +++ b/include/deal.II/grid/connectivity.h @@ -1164,9 +1164,8 @@ namespace internal { col_d[offset_i] = counter; orientations[offset_i] = - ReferenceCell::compute_orientation(ad_entity_types[offset_i], - ad_entity_vertices[offset_i], - ref_indices); + ad_entity_types[offset_i].compute_orientation( + ad_entity_vertices[offset_i], ref_indices); } } ptr_0.push_back(col_0.size()); diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index b7386fed1b..636e671abb 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -20,6 +20,7 @@ #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -73,6 +74,66 @@ namespace ReferenceCell */ Type(const CellKinds kind); + /** + * Return the dimension of the reference cell represented by the current + * object. + */ + unsigned int + get_dimension() const; + + /** + * Compute the value of the $i$-th linear shape function at location $\xi$ + * for the current reference-cell type. + */ + template + double + d_linear_shape_function(const Point &xi, const unsigned int i) const; + + /** + * Compute the gradient of the $i$-th linear shape function at location + * $\xi$ for the current reference-cell type. + */ + template + Tensor<1, dim> + d_linear_shape_function_gradient(const Point & xi, + const unsigned int i) const; + + /* + * Return $i$-th unit tangential vector of a face of the reference cell. + * The vectors are arranged such that the + * cross product between the two vectors returns the unit normal vector. + * + * @pre $i$ must be between zero and `dim-1`. + */ + template + Tensor<1, dim> + unit_tangential_vectors(const unsigned int face_no, + const unsigned int i) const; + + /** + * Determine the orientation of the current entity described by its + * vertices @p var_1 relative to an entity described by @p var_0. + */ + template + unsigned char + compute_orientation(const std::array &vertices_0, + const std::array &vertices_1) const; + + /** + * Inverse function of compute_orientation(). + */ + template + std::array + permute_according_orientation(const std::array &vertices, + const unsigned int orientation) const; + + /** + * Return a text representation of the reference cell represented by the + * current object. + */ + std::string + to_string() const; + /** * Conversion operator to an integer. */ @@ -259,13 +320,10 @@ namespace ReferenceCell - /** - * Return the dimension of the given reference-cell type @p type. - */ inline unsigned int - get_dimension(const Type &type) + Type::get_dimension() const { - switch (type) + switch (kind) { case Type::Vertex: return 0; @@ -284,13 +342,15 @@ namespace ReferenceCell } } + + /** * Convert the given reference cell type to a string. */ inline std::string - to_string(const Type &type) + Type::to_string() const { - switch (type) + switch (kind) { case Type::Vertex: return "Vertex"; @@ -393,21 +453,17 @@ namespace ReferenceCell } - /** - * Compute the value of the $i$-th linear shape function at location $\xi$ for - * a given reference-cell type. - */ + template inline double - d_linear_shape_function(const Type & reference_cell, - const Point & xi, - const unsigned int i) + Type::d_linear_shape_function(const Point & xi, + const unsigned int i) const { - if (reference_cell == get_hypercube(dim)) + AssertDimension(dim, get_dimension()); + if (*this == get_hypercube(dim)) return GeometryInfo::d_linear_shape_function(xi, i); - if (reference_cell == - Type::Tri) // see also Simplex::ScalarPolynomial::compute_value + if (*this == Type::Tri) // see also Simplex::ScalarPolynomial::compute_value { switch (i) { @@ -420,8 +476,7 @@ namespace ReferenceCell } } - if (reference_cell == - Type::Tet) // see also Simplex::ScalarPolynomial::compute_value + if (*this == Type::Tet) // see also Simplex::ScalarPolynomial::compute_value { switch (i) { @@ -437,19 +492,18 @@ namespace ReferenceCell } } - if (reference_cell == + if (*this == Type::Wedge) // see also Simplex::ScalarWedgePolynomial::compute_value { - return d_linear_shape_function<2>(Type::Tri, - Point<2>(xi[std::min(0, dim - 1)], - xi[std::min(1, dim - 1)]), - i % 3) * - d_linear_shape_function<1>(Type::Line, - Point<1>(xi[std::min(2, dim - 1)]), - i / 3); + return Type(Type::Tri).d_linear_shape_function<2>( + Point<2>(xi[std::min(0, dim - 1)], xi[std::min(1, dim - 1)]), + i % 3) * + Type(Type::Line) + .d_linear_shape_function<1>(Point<1>(xi[std::min(2, dim - 1)]), + i / 3); } - if (reference_cell == + if (*this == Type::Pyramid) // see also // Simplex::ScalarPyramidPolynomial::compute_value { @@ -486,21 +540,18 @@ namespace ReferenceCell return 0.0; } - /** - * Compute the gradient of the $i$-th linear shape function at location $\xi$ - * for a given reference-cell type. - */ + + template inline Tensor<1, dim> - d_linear_shape_function_gradient(const Type & reference_cell, - const Point & xi, - const unsigned int i) + Type::d_linear_shape_function_gradient(const Point & xi, + const unsigned int i) const { - if (reference_cell == get_hypercube(dim)) + AssertDimension(dim, get_dimension()); + if (*this == get_hypercube(dim)) return GeometryInfo::d_linear_shape_function_gradient(xi, i); - if (reference_cell == - Type::Tri) // see also Simplex::ScalarPolynomial::compute_grad + if (*this == Type::Tri) // see also Simplex::ScalarPolynomial::compute_grad { switch (i) { @@ -518,26 +569,21 @@ namespace ReferenceCell return Point(+0.0, +0.0, +0.0); } - /** - * Return i-th unit tangential vector of a face of the reference cell. - * The vectors are arranged such that the - * cross product between the two vectors returns the unit normal vector. - */ + template inline Tensor<1, dim> - unit_tangential_vectors(const Type & reference_cell, - const unsigned int face_no, - const unsigned int i) + Type::unit_tangential_vectors(const unsigned int face_no, + const unsigned int i) const { - AssertDimension(dim, get_dimension(reference_cell)); + AssertDimension(dim, get_dimension()); AssertIndexRange(i, dim - 1); - if (reference_cell == get_hypercube(dim)) + if (*this == get_hypercube(dim)) { AssertIndexRange(face_no, GeometryInfo::faces_per_cell); return GeometryInfo::unit_tangential_vectors[face_no][i]; } - else if (reference_cell == Type::Tri) + else if (*this == Type::Tri) { AssertIndexRange(face_no, 3); static const std::array, 3> table = { @@ -547,7 +593,7 @@ namespace ReferenceCell return table[face_no]; } - else if (reference_cell == Type::Tet) + else if (*this == Type::Tet) { AssertIndexRange(face_no, 4); static const std::array, 2>, 4> table = { @@ -563,7 +609,7 @@ namespace ReferenceCell return table[face_no][i]; } - else if (reference_cell == Type::Wedge) + else if (*this == Type::Wedge) { AssertIndexRange(face_no, 5); static const std::array, 2>, 5> table = { @@ -576,7 +622,7 @@ namespace ReferenceCell return table[face_no][i]; } - else if (reference_cell == Type::Pyramid) + else if (*this == Type::Pyramid) { AssertIndexRange(face_no, 5); static const std::array, 2>, 5> table = { @@ -605,7 +651,7 @@ namespace ReferenceCell inline Tensor<1, dim> unit_normal_vectors(const Type &reference_cell, const unsigned int face_no) { - AssertDimension(dim, get_dimension(reference_cell)); + AssertDimension(dim, reference_cell.get_dimension()); if (reference_cell == get_hypercube(dim)) { @@ -615,7 +661,7 @@ namespace ReferenceCell else if (dim == 2) { const auto tangential = - unit_tangential_vectors(reference_cell, face_no, 0); + reference_cell.unit_tangential_vectors(face_no, 0); Tensor<1, dim> result; @@ -627,8 +673,8 @@ namespace ReferenceCell else if (dim == 3) { return cross_product_3d( - unit_tangential_vectors(reference_cell, face_no, 0), - unit_tangential_vectors(reference_cell, face_no, 1)); + reference_cell.unit_tangential_vectors(face_no, 0), + reference_cell.unit_tangential_vectors(face_no, 1)); } Assert(false, ExcNotImplemented()); @@ -1945,19 +1991,14 @@ namespace ReferenceCell - /** - * Determine the orientation of an entity of @p type described by its - * vertices @p var_1 relative to an entity described by @p var_0. - */ template inline unsigned char - compute_orientation(const ReferenceCell::Type entity_type, - const std::array & vertices_0, - const std::array & vertices_1) + Type::compute_orientation(const std::array &vertices_0, + const std::array &vertices_1) const { AssertIndexRange( - ReferenceCell::internal::Info::get_cell(entity_type).n_vertices(), N + 1); - if (entity_type == ReferenceCell::Type::Line) + ReferenceCell::internal::Info::get_cell(*this).n_vertices(), N + 1); + if (*this == ReferenceCell::Type::Line) { const std::array i{{vertices_0[0], vertices_0[1]}}; const std::array j{{vertices_1[0], vertices_1[1]}}; @@ -1970,7 +2011,7 @@ namespace ReferenceCell if (i == std::array{{j[1], j[0]}}) return 0; } - else if (entity_type == ReferenceCell::Type::Tri) + else if (*this == ReferenceCell::Type::Tri) { const std::array i{{vertices_0[0], vertices_0[1], vertices_0[2]}}; const std::array j{{vertices_1[0], vertices_1[1], vertices_1[2]}}; @@ -1999,7 +2040,7 @@ namespace ReferenceCell if (i == std::array{{j[1], j[0], j[2]}}) return 4; } - else if (entity_type == ReferenceCell::Type::Quad) + else if (*this == ReferenceCell::Type::Quad) { const std::array i{ {vertices_0[0], vertices_0[1], vertices_0[2], vertices_0[3]}}; @@ -2039,25 +2080,22 @@ namespace ReferenceCell return 6; } - Assert( - false, - (internal::NoPermutation(entity_type, vertices_0, vertices_1))); + Assert(false, + (internal::NoPermutation(*this, vertices_0, vertices_1))); return -1; } - /** - * Inverse function of compute_orientation(). - */ + + template inline std::array - permute_according_orientation(const ReferenceCell::Type entity_type, - const std::array & vertices, - const unsigned int orientation) + Type::permute_according_orientation(const std::array &vertices, + const unsigned int orientation) const { std::array temp; - if (entity_type == ReferenceCell::Type::Line) + if (*this == ReferenceCell::Type::Line) { switch (orientation) { @@ -2071,7 +2109,7 @@ namespace ReferenceCell Assert(false, ExcNotImplemented()); } } - else if (entity_type == ReferenceCell::Type::Tri) + else if (*this == ReferenceCell::Type::Tri) { switch (orientation) { @@ -2097,7 +2135,7 @@ namespace ReferenceCell Assert(false, ExcNotImplemented()); } } - else if (entity_type == ReferenceCell::Type::Quad) + else if (*this == ReferenceCell::Type::Quad) { switch (orientation) { diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc index cdcef09bcd..85a997f514 100644 --- a/source/base/qprojector.cc +++ b/source/base/qprojector.cc @@ -647,10 +647,8 @@ QProjector<3>::project_to_all_faces( [](const auto &face, const auto &orientation) -> std::vector> { std::array, 3> vertices; std::copy_n(face.first.begin(), face.first.size(), vertices.begin()); - const auto temp = - ReferenceCell::permute_according_orientation(ReferenceCell::Type::Tri, - vertices, - orientation); + const auto temp = ReferenceCell::Type(ReferenceCell::Type::Tri) + .permute_according_orientation(vertices, orientation); return std::vector>(temp.begin(), temp.begin() + face.first.size()); }; @@ -659,10 +657,8 @@ QProjector<3>::project_to_all_faces( [](const auto &face, const auto &orientation) -> std::vector> { std::array, 4> vertices; std::copy_n(face.first.begin(), face.first.size(), vertices.begin()); - const auto temp = - ReferenceCell::permute_according_orientation(ReferenceCell::Type::Quad, - vertices, - orientation); + const auto temp = ReferenceCell::Type(ReferenceCell::Type::Quad) + .permute_according_orientation(vertices, orientation); return std::vector>(temp.begin(), temp.begin() + face.first.size()); }; diff --git a/source/fe/fe_nothing.cc b/source/fe/fe_nothing.cc index 44009a9e48..99b7e325be 100644 --- a/source/fe/fe_nothing.cc +++ b/source/fe/fe_nothing.cc @@ -76,8 +76,7 @@ FE_Nothing::get_name() const std::vector name_components; if (this->reference_cell_type() != ReferenceCell::get_hypercube(dim)) - name_components.push_back( - ReferenceCell::to_string(this->reference_cell_type())); + name_components.push_back(this->reference_cell_type().to_string()); if (this->n_components() > 1) name_components.push_back(std::to_string(this->n_components())); if (dominate) diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index be1982c505..3b8415b4d8 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -158,17 +158,17 @@ MappingFE::InternalData::initialize_face( for (unsigned int i = 0; i < n_faces; ++i) { unit_tangentials[i].resize(n_original_q_points); - std::fill(unit_tangentials[i].begin(), - unit_tangentials[i].end(), - ReferenceCell::unit_tangential_vectors( - reference_cell_type, i, 0)); + std::fill( + unit_tangentials[i].begin(), + unit_tangentials[i].end(), + reference_cell_type.template unit_tangential_vectors(i, 0)); if (dim > 2) { unit_tangentials[n_faces + i].resize(n_original_q_points); std::fill(unit_tangentials[n_faces + i].begin(), unit_tangentials[n_faces + i].end(), - ReferenceCell::unit_tangential_vectors( - reference_cell_type, i, 1)); + reference_cell_type + .template unit_tangential_vectors(i, 1)); } } } @@ -871,9 +871,8 @@ MappingFE::MappingFE(const FiniteElement &fe) for (unsigned int point = 0; point < n_points; ++point) for (unsigned int i = 0; i < n_shape_functions; ++i) mapping_support_point_weights(point, i) = - ReferenceCell::d_linear_shape_function(reference_cell_type, - mapping_support_points[point], - i); + reference_cell_type.d_linear_shape_function( + mapping_support_points[point], i); } diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc index d440769a9a..bb21ab4fa8 100644 --- a/source/fe/mapping_fe_field.cc +++ b/source/fe/mapping_fe_field.cc @@ -580,16 +580,16 @@ MappingFEField::compute_face_data( data.unit_tangentials[i].resize(n_original_q_points); std::fill(data.unit_tangentials[i].begin(), data.unit_tangentials[i].end(), - ReferenceCell::unit_tangential_vectors( - reference_cell_type, i, 0)); + reference_cell_type + .template unit_tangential_vectors(i, 0)); if (dim > 2) { data.unit_tangentials[n_faces + i].resize( n_original_q_points); std::fill(data.unit_tangentials[n_faces + i].begin(), data.unit_tangentials[n_faces + i].end(), - ReferenceCell::unit_tangential_vectors( - reference_cell_type, i, 1)); + reference_cell_type + .template unit_tangential_vectors(i, 1)); } } } diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index 2fe79620a4..4dba26b001 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -913,16 +913,17 @@ FlatManifold::normal_vector( { Point F; for (const unsigned int v : face->vertex_indices()) - F += face->vertex(v) * ReferenceCell::d_linear_shape_function( - face_reference_cell_type, xi, v); + F += face->vertex(v) * + face_reference_cell_type.d_linear_shape_function(xi, v); for (unsigned int i = 0; i < facedim; ++i) { grad_F[i] = 0; for (const unsigned int v : face->vertex_indices()) grad_F[i] += - face->vertex(v) * ReferenceCell::d_linear_shape_function_gradient( - face_reference_cell_type, xi, v)[i]; + face->vertex(v) * + face_reference_cell_type.d_linear_shape_function_gradient(xi, + v)[i]; } Tensor<1, facedim> J; diff --git a/source/grid/reference_cell.cc b/source/grid/reference_cell.cc index 8ae0afa4f4..d7b2e33f1a 100644 --- a/source/grid/reference_cell.cc +++ b/source/grid/reference_cell.cc @@ -37,7 +37,7 @@ namespace ReferenceCell make_triangulation(const Type & reference_cell, Triangulation &tria) { - AssertDimension(dim, get_dimension(reference_cell)); + AssertDimension(dim, reference_cell.get_dimension()); if (reference_cell == get_hypercube(dim)) { @@ -110,7 +110,7 @@ namespace ReferenceCell std::unique_ptr> get_default_mapping(const Type &reference_cell, const unsigned int degree) { - AssertDimension(dim, get_dimension(reference_cell)); + AssertDimension(dim, reference_cell.get_dimension()); if (reference_cell == get_hypercube(dim)) return std::make_unique>(degree); @@ -136,7 +136,7 @@ namespace ReferenceCell const Mapping & get_default_linear_mapping(const Type &reference_cell) { - AssertDimension(dim, get_dimension(reference_cell)); + AssertDimension(dim, reference_cell.get_dimension()); if (reference_cell == get_hypercube(dim)) { @@ -195,7 +195,7 @@ namespace ReferenceCell get_gauss_type_quadrature(const Type & reference_cell, const unsigned n_points_1D) { - AssertDimension(dim, get_dimension(reference_cell)); + AssertDimension(dim, reference_cell.get_dimension()); if (reference_cell == get_hypercube(dim)) return QGauss(n_points_1D); @@ -215,7 +215,7 @@ namespace ReferenceCell Quadrature & get_nodal_type_quadrature(const Type &reference_cell) { - AssertDimension(dim, get_dimension(reference_cell)); + AssertDimension(dim, reference_cell.get_dimension()); const auto create_quadrature = [](const Type &reference_cell) { Triangulation tria; diff --git a/tests/simplex/unit_tangential_vectors_01.cc b/tests/simplex/unit_tangential_vectors_01.cc index fae80a25ba..99f5a36ba4 100644 --- a/tests/simplex/unit_tangential_vectors_01.cc +++ b/tests/simplex/unit_tangential_vectors_01.cc @@ -33,9 +33,7 @@ test(const ReferenceCell::Type &reference_cell) for (const auto face_no : ReferenceCell::internal::Info::get_cell(reference_cell).face_indices()) { - deallog << ReferenceCell::unit_normal_vectors(reference_cell, - face_no) - << std::endl; + deallog << reference_cell.unit_normal_vectors(face_no) << std::endl; for (unsigned int i = 0; i < dim - 1; ++i) deallog << ReferenceCell::unit_tangential_vectors(reference_cell, -- 2.39.5