From 2edb9bf1b284193cc7512ce2223a51f9f0e03f08 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 8 Feb 2021 14:06:24 -0700 Subject: [PATCH] Move some internal::ReferenceCell::Base functions into ReferenceCell itself. --- include/deal.II/fe/fe_tools.templates.h | 38 ++-- include/deal.II/grid/reference_cell.h | 176 +++++++++++++++++- .../deal.II/grid/tria_accessor.templates.h | 6 +- .../matrix_free/shape_info.templates.h | 11 +- source/base/data_out_base.cc | 11 +- source/dofs/dof_tools.cc | 17 +- source/fe/fe.cc | 25 ++- source/fe/fe_data.cc | 45 ++--- source/fe/fe_values.cc | 3 +- source/fe/mapping_fe.cc | 8 +- source/fe/mapping_fe_field.cc | 3 +- 11 files changed, 232 insertions(+), 111 deletions(-) diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h index 2b90de16c8..7952648c1b 100644 --- a/include/deal.II/fe/fe_tools.templates.h +++ b/include/deal.II/fe/fe_tools.templates.h @@ -282,11 +282,9 @@ namespace FETools // 3. Quads for (unsigned int quad_number = 0; quad_number < - (dim == 2 ? 1 : - (dim == 3 ? dealii::internal::ReferenceCell::get_cell( - fes.front()->reference_cell()) - .n_faces() : - 0)); + (dim == 2 ? + 1 : + (dim == 3 ? fes.front()->reference_cell().n_faces() : 0)); ++quad_number) { for (unsigned int base = 0; base < fes.size(); ++base) @@ -499,11 +497,9 @@ namespace FETools // 3. Quads for (unsigned int quad_number = 0; quad_number < - (dim == 2 ? 1 : - (dim == 3 ? dealii::internal::ReferenceCell::get_cell( - fes.front()->reference_cell()) - .n_faces() : - 0)); + (dim == 2 ? + 1 : + (dim == 3 ? fes.front()->reference_cell().n_faces() : 0)); ++quad_number) { unsigned int comp_start = 0; @@ -761,11 +757,7 @@ namespace FETools // 3. Quads for (unsigned int quad_number = 0; quad_number < - (dim == 2 ? 1 : - (dim == 3 ? dealii::internal::ReferenceCell::get_cell( - fe.reference_cell()) - .n_faces() : - 0)); + (dim == 2 ? 1 : (dim == 3 ? fe.reference_cell().n_faces() : 0)); ++quad_number) { unsigned int comp_start = 0; @@ -872,9 +864,7 @@ namespace FETools unsigned int total_index = 0; for (unsigned int vertex_number = 0; vertex_number < - dealii::internal::ReferenceCell::get_face(fe.reference_cell(), - face_no) - .n_vertices(); + fe.reference_cell().face_reference_cell(face_no).n_vertices(); ++vertex_number) { unsigned int comp_start = 0; @@ -932,9 +922,7 @@ namespace FETools // 2. Lines for (unsigned int line_number = 0; line_number < - dealii::internal::ReferenceCell::get_face(fe.reference_cell(), - face_no) - .n_lines(); + fe.reference_cell().face_reference_cell(face_no).n_lines(); ++line_number) { unsigned int comp_start = 0; @@ -1995,9 +1983,7 @@ namespace FETools { unsigned int face_dof = 0; for (unsigned int i = 0; - i < dealii::internal::ReferenceCell::get_face(fe.reference_cell(), - face_no) - .n_vertices(); + i < fe.reference_cell().face_reference_cell(face_no).n_vertices(); ++i) { const unsigned int offset_c = @@ -2015,9 +2001,7 @@ namespace FETools } for (unsigned int i = 1; - i <= dealii::internal::ReferenceCell::get_face(fe.reference_cell(), - face_no) - .n_lines(); + i <= fe.reference_cell().face_reference_cell(face_no).n_lines(); ++i) { const unsigned int offset_c = diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index 4c8adc8536..83718e97bc 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -203,6 +203,54 @@ public: const Quadrature & get_nodal_type_quadrature() const; + /** + * @} + */ + + /** + * @name Querying the number of building blocks of a reference cell + * @{ + */ + + /** + * Return the number of vertices that make up the reference + * cell in question. A vertex is a "corner" (a zero-dimensional + * object) of the reference cell. + */ + unsigned int + n_vertices() const; + + /** + * Return the number of lines that make up the reference + * cell in question. A line is an "edge" (a one-dimensional + * object) of the reference cell. + */ + unsigned int + n_lines() const; + + /** + * Return the number of faces that make up the reference + * cell in question. A face is a `(dim-1)`-dimensional + * object bounding the reference cell. + */ + unsigned int + n_faces() const; + + /** + * Return the reference-cell type of face @p face_no of the current + * object. For example, if the current object is + * ReferenceCells::Tetrahedron, then `face_no` must be between + * in the interval $[0,4)$ and the function will always return + * ReferenceCells::Triangle. If the current object is + * ReferenceCells::Hexahedron, then `face_no` must be between + * in the interval $[0,6)$ and the function will always return + * ReferenceCells::Quadrilateral. For wedges and pyramids, the + * returned object may be either ReferenceCells::Triangle or + * ReferenceCells::Quadrilateral, depending on the given index. + */ + ReferenceCell + face_reference_cell(const unsigned int face_no) const; + /** * @} */ @@ -550,6 +598,122 @@ ReferenceCell::get_dimension() const +inline unsigned int +ReferenceCell::n_vertices() const +{ + if (*this == ReferenceCells::Vertex) + return 1; + else if (*this == ReferenceCells::Line) + return 2; + else if (*this == ReferenceCells::Triangle) + return 3; + else if (*this == ReferenceCells::Quadrilateral) + return 4; + else if (*this == ReferenceCells::Tetrahedron) + return 4; + else if (*this == ReferenceCells::Pyramid) + return 5; + else if (*this == ReferenceCells::Wedge) + return 6; + else if (*this == ReferenceCells::Hexahedron) + return 8; + + Assert(false, ExcNotImplemented()); + return numbers::invalid_unsigned_int; +} + + + +inline unsigned int +ReferenceCell::n_lines() const +{ + if (*this == ReferenceCells::Vertex) + return 0; + else if (*this == ReferenceCells::Line) + return 1; + else if (*this == ReferenceCells::Triangle) + return 3; + else if (*this == ReferenceCells::Quadrilateral) + return 4; + else if (*this == ReferenceCells::Tetrahedron) + return 6; + else if (*this == ReferenceCells::Pyramid) + return 7; + else if (*this == ReferenceCells::Wedge) + return 9; + else if (*this == ReferenceCells::Hexahedron) + return 12; + + Assert(false, ExcNotImplemented()); + return numbers::invalid_unsigned_int; +} + + + +inline unsigned int +ReferenceCell::n_faces() const +{ + if (*this == ReferenceCells::Vertex) + return 0; + else if (*this == ReferenceCells::Line) + return 2; + else if (*this == ReferenceCells::Triangle) + return 3; + else if (*this == ReferenceCells::Quadrilateral) + return 4; + else if (*this == ReferenceCells::Tetrahedron) + return 4; + else if (*this == ReferenceCells::Pyramid) + return 5; + else if (*this == ReferenceCells::Wedge) + return 5; + else if (*this == ReferenceCells::Hexahedron) + return 6; + + Assert(false, ExcNotImplemented()); + return numbers::invalid_unsigned_int; +} + + + +inline ReferenceCell +ReferenceCell::face_reference_cell(const unsigned int face_no) const +{ + AssertIndexRange(face_no, n_faces()); + + if (*this == ReferenceCells::Vertex) + return ReferenceCells::Invalid; + else if (*this == ReferenceCells::Line) + return ReferenceCells::Vertex; + else if (*this == ReferenceCells::Triangle) + return ReferenceCells::Line; + else if (*this == ReferenceCells::Quadrilateral) + return ReferenceCells::Line; + else if (*this == ReferenceCells::Tetrahedron) + return ReferenceCells::Triangle; + else if (*this == ReferenceCells::Pyramid) + { + if (face_no == 0) + return ReferenceCells::Quadrilateral; + else + return ReferenceCells::Triangle; + } + else if (*this == ReferenceCells::Wedge) + { + if (face_no > 1) + return ReferenceCells::Quadrilateral; + else + return ReferenceCells::Triangle; + } + else if (*this == ReferenceCells::Hexahedron) + return ReferenceCells::Quadrilateral; + + Assert(false, ExcNotImplemented()); + return ReferenceCells::Invalid; +} + + + namespace ReferenceCells { template @@ -887,6 +1051,7 @@ namespace internal */ virtual ~Base() = default; + private: /** * Number of vertices. */ @@ -918,6 +1083,7 @@ namespace internal return 0; } + public: /** * Return an object that can be thought of as an array containing all * indices from zero to n_vertices(). @@ -1027,6 +1193,7 @@ namespace internal return true; } + private: /** * Return reference-cell type of face @p face_no. */ @@ -1039,6 +1206,7 @@ namespace internal return ReferenceCells::Invalid; } + public: /** * Map face line number to cell line number. */ @@ -2011,7 +2179,7 @@ namespace internal inline const internal::ReferenceCell::Base & get_face(const dealii::ReferenceCell &type, const unsigned int face_no) { - return get_cell(get_cell(type).face_reference_cell(face_no)); + return get_cell(type.face_reference_cell(face_no)); } } // namespace ReferenceCell @@ -2049,8 +2217,7 @@ namespace internal { out << "["; - const unsigned int n_vertices = - internal::ReferenceCell::get_cell(entity_type).n_vertices(); + const unsigned int n_vertices = entity_type.n_vertices(); for (unsigned int i = 0; i < n_vertices; ++i) { @@ -2095,8 +2262,7 @@ inline unsigned char ReferenceCell::compute_orientation(const std::array &vertices_0, const std::array &vertices_1) const { - AssertIndexRange(internal::ReferenceCell::get_cell(*this).n_vertices(), - N + 1); + AssertIndexRange(n_vertices(), N + 1); if (*this == ReferenceCells::Line) { const std::array i{{vertices_0[0], vertices_0[1]}}; diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h index b99b032003..a59999e0b0 100644 --- a/include/deal.II/grid/tria_accessor.templates.h +++ b/include/deal.II/grid/tria_accessor.templates.h @@ -2195,7 +2195,7 @@ template unsigned int TriaAccessor::n_vertices() const { - return this->reference_cell_info().n_vertices(); + return this->reference_cell().n_vertices(); } @@ -2204,7 +2204,7 @@ template unsigned int TriaAccessor::n_lines() const { - return this->reference_cell_info().n_lines(); + return this->reference_cell().n_lines(); } @@ -2215,7 +2215,7 @@ TriaAccessor::n_faces() const { AssertDimension(structdim, dim); - return this->reference_cell_info().n_faces(); + return this->reference_cell().n_faces(); } diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 4e77bb7e33..53618be269 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -344,10 +344,10 @@ namespace internal for (unsigned int face_no = 0; face_no < quad_face.size(); ++face_no) - n_max_vertices = std::max(n_max_vertices, - internal::ReferenceCell::get_face( - reference_cell_type, face_no) - .n_vertices()); + n_max_vertices = + std::max(n_max_vertices, + reference_cell_type.face_reference_cell(face_no) + .n_vertices()); const auto projected_quad_face = QProjector::project_to_all_faces(reference_cell_type, @@ -369,8 +369,7 @@ namespace internal { const unsigned int n_face_orientations = dim == 2 ? 2 : - (2 * internal::ReferenceCell::get_face( - reference_cell_type, f) + (2 * reference_cell_type.face_reference_cell(f) .n_vertices()); const unsigned int n_q_points_face = quad_face[f].size(); diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc index 746e2d9653..e9cdbe4b2e 100644 --- a/source/base/data_out_base.cc +++ b/source/base/data_out_base.cc @@ -927,9 +927,7 @@ namespace else { Assert(patch.n_subdivisions == 1, ExcNotImplemented()); - const auto &info = - internal::ReferenceCell::get_cell(patch.reference_cell); - n_nodes += info.n_vertices(); + n_nodes += patch.reference_cell.n_vertices(); n_cells += 1; } } @@ -7739,9 +7737,6 @@ DataOutBase::write_hdf5_parallel( // patches Assert(patches.size() > 0, ExcNoPatches()); - const auto &cell_info = - internal::ReferenceCell::get_cell(patches[0].reference_cell); - hid_t h5_mesh_file_id = -1, h5_solution_file_id, file_plist_id, plist_id; hid_t node_dataspace, node_dataset, node_file_dataspace, node_memory_dataspace; @@ -7835,7 +7830,7 @@ DataOutBase::write_hdf5_parallel( AssertThrow(node_dataspace >= 0, ExcIO()); cell_ds_dim[0] = global_node_cell_count[1]; - cell_ds_dim[1] = cell_info.n_vertices(); + cell_ds_dim[1] = patches[0].reference_cell.n_vertices(); cell_dataspace = H5Screate_simple(2, cell_ds_dim, nullptr); AssertThrow(cell_dataspace >= 0, ExcIO()); @@ -7896,7 +7891,7 @@ DataOutBase::write_hdf5_parallel( // And repeat for cells count[0] = local_node_cell_count[1]; - count[1] = cell_info.n_vertices(); + count[1] = patches[0].reference_cell.n_vertices(); offset[0] = global_node_cell_offsets[1]; offset[1] = 0; cell_memory_dataspace = H5Screate_simple(2, count, nullptr); diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc index d09735a6a0..30fc0849fe 100644 --- a/source/dofs/dof_tools.cc +++ b/source/dofs/dof_tools.cc @@ -689,16 +689,13 @@ namespace DoFTools const auto reference_cell = cell->reference_cell(); - const auto &cell_rc = - dealii::internal::ReferenceCell::get_cell(reference_cell); - const auto &face_rc = - dealii::internal::ReferenceCell::get_face(reference_cell, - face); - - const unsigned int n_vertices_per_cell = cell_rc.n_vertices(); - const unsigned int n_lines_per_cell = cell_rc.n_lines(); - const unsigned int n_vertices_per_face = face_rc.n_vertices(); - const unsigned int n_lines_per_face = face_rc.n_lines(); + const unsigned int n_vertices_per_cell = + reference_cell.n_vertices(); + const unsigned int n_lines_per_cell = reference_cell.n_lines(); + const unsigned int n_vertices_per_face = + reference_cell.face_reference_cell(face).n_vertices(); + const unsigned int n_lines_per_face = + reference_cell.face_reference_cell(face).n_lines(); const unsigned int dofs_per_face = fe.n_dofs_per_face(face); face_dof_indices.resize(dofs_per_face); diff --git a/source/fe/fe.cc b/source/fe/fe.cc index 2bd34699a4..9d40d929e3 100644 --- a/source/fe/fe.cc +++ b/source/fe/fe.cc @@ -156,12 +156,12 @@ FiniteElement::FiniteElement( for (unsigned int f = 0; f < this->n_unique_quads(); ++f) { - adjust_quad_dof_index_for_face_orientation_table[f] = Table<2, int>( - this->n_dofs_per_quad(f), - internal::ReferenceCell::get_cell(this->reference_cell()) - .face_reference_cell(f) == ReferenceCells::Quadrilateral ? - 8 : - 6); + adjust_quad_dof_index_for_face_orientation_table[f] = + Table<2, int>(this->n_dofs_per_quad(f), + this->reference_cell().face_reference_cell(f) == + ReferenceCells::Quadrilateral ? + 8 : + 6); adjust_quad_dof_index_for_face_orientation_table[f].fill(0); } } @@ -575,7 +575,7 @@ FiniteElement::face_to_cell_index(const unsigned int face_index, internal::ReferenceCell::get_cell(this->reference_cell()); AssertIndexRange(face_index, this->n_dofs_per_face(face)); - AssertIndexRange(face, refence_cell.n_faces()); + AssertIndexRange(face, this->reference_cell().n_faces()); // TODO: we could presumably solve the 3d case below using the // adjust_quad_dof_index_for_face_orientation_table field. for the @@ -684,12 +684,11 @@ FiniteElement::adjust_quad_dof_index_for_face_orientation( AssertIndexRange(index, this->n_dofs_per_quad(face)); Assert(adjust_quad_dof_index_for_face_orientation_table [this->n_unique_quads() == 1 ? 0 : face] - .n_elements() == - (internal::ReferenceCell::get_cell(this->reference_cell()) - .face_reference_cell(face) == ReferenceCells::Quadrilateral ? - 8 : - 6) * - this->n_dofs_per_quad(face), + .n_elements() == (this->reference_cell().face_reference_cell( + face) == ReferenceCells::Quadrilateral ? + 8 : + 6) * + this->n_dofs_per_quad(face), ExcInternalError()); return index + adjust_quad_dof_index_for_face_orientation_table diff --git a/source/fe/fe_data.cc b/source/fe/fe_data.cc index 80c7efbb09..b7a980e6cb 100644 --- a/source/fe/fe_data.cc +++ b/source/fe/fe_data.cc @@ -54,65 +54,50 @@ namespace internal // first_line_index const unsigned int first_line_index = - (internal::ReferenceCell::get_cell(cell_type).n_vertices() * - dofs_per_vertex); + (cell_type.n_vertices() * dofs_per_vertex); result.object_index[1][0] = first_line_index; // first_quad_index const unsigned int first_quad_index = - (first_line_index + - internal::ReferenceCell::get_cell(cell_type).n_lines() * dofs_per_line); + (first_line_index + cell_type.n_lines() * dofs_per_line); result.object_index[2][0] = first_quad_index; // first_hex_index result.object_index[3][0] = (first_quad_index + - (dim == 2 ? - 1 : - (dim == 3 ? internal::ReferenceCell::get_cell(cell_type).n_faces() : - 0)) * - dofs_per_quad); + (dim == 2 ? 1 : (dim == 3 ? cell_type.n_faces() : 0)) * dofs_per_quad); // first_face_line_index result.first_object_index_on_face[1][0] = - (internal::ReferenceCell::get_face(cell_type, face_no).n_vertices() * - dofs_per_vertex); + (cell_type.face_reference_cell(face_no).n_vertices() * dofs_per_vertex); // first_face_quad_index result.first_object_index_on_face[2][0] = - ((dim == 3 ? - internal::ReferenceCell::get_face(cell_type, face_no).n_vertices() * - dofs_per_vertex : - internal::ReferenceCell::get_cell(cell_type).n_vertices() * - dofs_per_vertex) + - internal::ReferenceCell::get_face(cell_type, face_no).n_lines() * - dofs_per_line); + ((dim == 3 ? cell_type.face_reference_cell(face_no).n_vertices() * + dofs_per_vertex : + cell_type.n_vertices() * dofs_per_vertex) + + cell_type.face_reference_cell(face_no).n_lines() * dofs_per_line); // dofs_per_face result.dofs_per_object_inclusive[dim - 1][0] = - (internal::ReferenceCell::get_face(cell_type, face_no).n_vertices() * - dofs_per_vertex + - internal::ReferenceCell::get_face(cell_type, face_no).n_lines() * - dofs_per_line + + (cell_type.face_reference_cell(face_no).n_vertices() * dofs_per_vertex + + cell_type.face_reference_cell(face_no).n_lines() * dofs_per_line + (dim == 3 ? 1 : 0) * dofs_per_quad); // dofs_per_cell result.dofs_per_object_inclusive[dim][0] = - (internal::ReferenceCell::get_cell(cell_type).n_vertices() * - dofs_per_vertex + - internal::ReferenceCell::get_cell(cell_type).n_lines() * dofs_per_line + - (dim == 2 ? - 1 : - (dim == 3 ? internal::ReferenceCell::get_cell(cell_type).n_faces() : - 0)) * - dofs_per_quad + + (cell_type.n_vertices() * dofs_per_vertex + + cell_type.n_lines() * dofs_per_line + + (dim == 2 ? 1 : (dim == 3 ? cell_type.n_faces() : 0)) * dofs_per_quad + (dim == 3 ? 1 : 0) * dofs_per_hex); return result; } } // namespace internal + + template FiniteElementData::FiniteElementData( const std::vector &dofs_per_object, diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc index 70d3362853..6284724fe2 100644 --- a/source/fe/fe_values.cc +++ b/source/fe/fe_values.cc @@ -4637,8 +4637,7 @@ FEFaceValuesBase::FEFaceValuesBase( , quadrature(quadrature) { Assert(quadrature.size() == 1 || - quadrature.size() == - internal::ReferenceCell::get_cell(fe.reference_cell()).n_faces(), + quadrature.size() == fe.reference_cell().n_faces(), ExcInternalError()); } diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index d7fffde4d1..577aa3924d 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -154,8 +154,7 @@ MappingFE::InternalData::initialize_face( // Compute tangentials to the unit cell. const auto reference_cell = this->fe.reference_cell(); - const auto n_faces = - internal::ReferenceCell::get_cell(reference_cell).n_faces(); + const auto n_faces = reference_cell.n_faces(); for (unsigned int i = 0; i < n_faces; ++i) { @@ -862,9 +861,8 @@ MappingFE::MappingFE(const FiniteElement &fe) const auto reference_cell = fe.reference_cell(); - const unsigned int n_points = mapping_support_points.size(); - const unsigned int n_shape_functions = - internal::ReferenceCell::get_cell(reference_cell).n_vertices(); + const unsigned int n_points = mapping_support_points.size(); + const unsigned int n_shape_functions = reference_cell.n_vertices(); this->mapping_support_point_weights = Table<2, double>(n_points, n_shape_functions); diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc index 9859e547ec..52fd7f8f4d 100644 --- a/source/fe/mapping_fe_field.cc +++ b/source/fe/mapping_fe_field.cc @@ -591,8 +591,7 @@ MappingFEField::compute_face_data( // TODO: only a single reference cell type possible... const auto reference_cell = this->euler_dof_handler->get_fe().reference_cell(); - const auto n_faces = - internal::ReferenceCell::get_cell(reference_cell).n_faces(); + const auto n_faces = reference_cell.n_faces(); // Compute tangentials to the unit cell. for (unsigned int i = 0; i < n_faces; ++i) -- 2.39.5