From 793f712d61a18bd25df5deefe7e3f3f5bc94d067 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 31 Jul 2020 22:20:28 +0200 Subject: [PATCH] Generalize FiniteElementData s.t. faces can different nr. of dofs --- include/deal.II/fe/fe_base.h | 212 +++++++++++++++++++++++++++++++---- source/fe/fe_data.cc | 195 +++++++++++++++++++++++--------- 2 files changed, 332 insertions(+), 75 deletions(-) diff --git a/include/deal.II/fe/fe_base.h b/include/deal.II/fe/fe_base.h index c8db30f119..8c126dff2a 100644 --- a/include/deal.II/fe/fe_base.h +++ b/include/deal.II/fe/fe_base.h @@ -128,6 +128,67 @@ namespace FiniteElementDomination inline Domination operator&(const Domination d1, const Domination d2); } // namespace FiniteElementDomination +namespace internal +{ + /** + * Internal data structure for setting up FiniteElementData. It stores for + * each object the (inclusive/exclusive) number of degrees of freedoms, as + * well as, the index of its first degree of freedom within a cell and the + * index of the first d-dimensional object within each face. + * + * The information is saved as a vector of vectors. One can query the + * inclusive number of dofs of the i-th d-dimensional object via: + * dofs_per_object_inclusive[d][i]. + * + * As an example, the data is shown for a quadratic wedge. Which consists of + * 6 vertices, 9 lines, and 5 faces (two triangles and three quadrilaterals). + * + * vertices lines faces cell + * dpo_excl 1 1 1 1 1 1 | 1 1 1 1 1 1 1 1 1 | 0 0 1 1 1 | 0 + * dpo_incl 1 1 1 1 1 1 | 3 3 3 3 3 3 3 3 3 | 6 6 9 9 9 | 18 + * obj_index 0 1 2 3 4 5 | 6 7 8 9 10 11 12 13 14 | 15 15 15 16 17 | 18 + * + * Since the above table looks as follows for: + * + * - a triangle: + * + * dpo_excl 1 1 1 | 1 1 1 | 0 + * obj_index 0 1 2 | 3 4 5 | 6 + * + * - quadrilateral: + * + * dpo_excl 1 1 1 1 | 1 1 1 1 | 1 + * obj_index 0 1 2 3 | 4 5 6 7 | 8 + * + * The index of the first d-dimensional object within each face results as: + * + * vertices lines face + * first_obj_index_on_face 0 0 0 0 0 | 3 3 4 4 4 | 6 6 8 8 8 + * + */ + struct GenericDoFsPerObject + { + /** + * Exclusive number of degrees of freedom per object. + */ + std::vector> dofs_per_object_exclusive; + + /** + * Inclusive number of degrees of freedom per object. + */ + std::vector> dofs_per_object_inclusive; + + /** + * First index of an object. + */ + std::vector> object_index; + + /** + * First index of an object within a face. + */ + std::vector> first_object_index_on_face; + }; +} // namespace internal /** * A class that declares a number of scalar constant variables that describe @@ -225,6 +286,18 @@ private: */ const ReferenceCell::Type cell_type; + /** + * Number of unique quads. If all quads have the same type, the value is + * one; else it equals the number of quads. + */ + const unsigned int number_unique_quads; + + /** + * Number of unique faces. If all faces have the same type, the value is + * one; else it equals the number of faces. + */ + const unsigned int number_unique_faces; + public: /** * Number of degrees of freedom on a vertex. @@ -237,6 +310,14 @@ public: */ const unsigned int dofs_per_line; +private: + /** + * Number of degrees of freedom on quads. If all quads have the same + * number of degrees freedoms the values equal dofs_per_quad. + */ + const std::vector n_dofs_on_quad; + +public: /** * Number of degrees of freedom in a quadrilateral; not including the * degrees of freedom on the lines and vertices of the quadrilateral. @@ -261,6 +342,15 @@ public: */ const unsigned int first_line_index; +private: + /** + * First index of a quad. If all quads have the same number of degrees of + * freedom, only the first index of the first quad is stored since the + * indices of the others can be simply recomputed. + */ + const std::vector first_index_of_quads; + +public: /** * First index of dof on a quad. */ @@ -271,16 +361,38 @@ public: */ const unsigned int first_hex_index; +private: + /** + * Index of the first line of all faces. + */ + const std::vector first_line_index_of_faces; + +public: /** * First index of dof on a line for face data. */ const unsigned int first_face_line_index; +private: + /** + * Index of the first quad of all faces. + */ + const std::vector first_quad_index_of_faces; + +public: /** * First index of dof on a quad for face data. */ const unsigned int first_face_quad_index; +private: + /** + * Number of degrees of freedom on faces. If all faces have the same + * number of degrees freedoms the values equal dofs_per_quad. + */ + const std::vector n_dofs_on_face; + +public: /** * Number of degrees of freedom on a face. This is the accumulated number of * degrees of freedom on all the objects of dimension up to dim-1 @@ -385,12 +497,39 @@ public: const Conformity conformity = unknown, const BlockIndices &block_indices = BlockIndices()); + /** + * The same as above but instead of passing a vector containing the degrees + * of freedoms per object a struct of type GenericDoFsPerObject. This allows + * that 2D objects might have different number of degrees of freedoms, which + * is particular useful for cells with triangles and quadrilaterals as faces. + */ + FiniteElementData(const internal::GenericDoFsPerObject &data, + const ReferenceCell::Type cell_type, + const unsigned int n_components, + const unsigned int degree, + const Conformity conformity = unknown, + const BlockIndices &block_indices = BlockIndices()); + /** * Return type of reference cell. */ ReferenceCell::Type reference_cell_type() const; + /** + * Number of unique quads. If all quads have the same type, the value is + * one; else it equals the number of quads. + */ + unsigned int + n_unique_quads() const; + + /** + * Number of unique faces. If all faces have the same type, the value is + * one; else it equals the number of faces. + */ + unsigned int + n_unique_faces() const; + /** * Number of dofs per vertex. */ @@ -407,7 +546,7 @@ public: * Number of dofs per quad. Not including dofs on lower dimensional objects. */ unsigned int - n_dofs_per_quad() const; + n_dofs_per_quad(unsigned int face_no = 0) const; /** * Maximum number of dofs per quad. Not including dofs on lower dimensional @@ -427,7 +566,7 @@ public: * dimensional objects. */ unsigned int - n_dofs_per_face() const; + n_dofs_per_face(unsigned int face_no = 0, unsigned int child = 0) const; /** * Maximum number of dofs per face, accumulating degrees of freedom of all @@ -453,7 +592,7 @@ public: */ template unsigned int - n_dofs_per_object() const; + n_dofs_per_object(const unsigned int i = 0) const; /** * Number of components. See @@ -511,7 +650,7 @@ public: * Return first index of dof on a quad. */ unsigned int - get_first_quad_index() const; + get_first_quad_index(const unsigned int quad_no = 0) const; /** * Return first index of dof on a hexahedron. @@ -523,13 +662,13 @@ public: * Return first index of dof on a line for face data. */ unsigned int - get_first_face_line_index() const; + get_first_face_line_index(const unsigned int face_no = 0) const; /** * Return first index of dof on a quad for face data. */ unsigned int - get_first_face_quad_index() const; + get_first_face_quad_index(const unsigned int face_no = 0) const; }; @@ -595,6 +734,25 @@ FiniteElementData::reference_cell_type() const } + +template +inline unsigned int +FiniteElementData::n_unique_quads() const +{ + return number_unique_quads; +} + + + +template +inline unsigned int +FiniteElementData::n_unique_faces() const +{ + return number_unique_faces; +} + + + template inline unsigned int FiniteElementData::n_dofs_per_vertex() const @@ -615,9 +773,9 @@ FiniteElementData::n_dofs_per_line() const template inline unsigned int -FiniteElementData::n_dofs_per_quad() const +FiniteElementData::n_dofs_per_quad(unsigned int face_no) const { - return dofs_per_quad; + return n_dofs_on_quad[n_dofs_on_quad.size() == 1 ? 0 : face_no]; } @@ -642,9 +800,12 @@ FiniteElementData::n_dofs_per_hex() const template inline unsigned int -FiniteElementData::n_dofs_per_face() const +FiniteElementData::n_dofs_per_face(unsigned int face_no, + unsigned int child_no) const { - return dofs_per_face; + (void)child_no; + + return n_dofs_on_face[n_dofs_on_face.size() == 1 ? 0 : face_no]; } @@ -670,18 +831,18 @@ FiniteElementData::n_dofs_per_cell() const template template inline unsigned int -FiniteElementData::n_dofs_per_object() const +FiniteElementData::n_dofs_per_object(const unsigned int i) const { switch (structdim) { case 0: - return dofs_per_vertex; + return n_dofs_per_vertex(); case 1: - return dofs_per_line; + return n_dofs_per_line(); case 2: - return dofs_per_quad; + return n_dofs_per_quad((structdim == 2 && dim == 3) ? i : 0); case 3: - return dofs_per_hex; + return n_dofs_per_hex(); default: Assert(false, ExcInternalError()); } @@ -743,9 +904,12 @@ FiniteElementData::get_first_line_index() const template unsigned int -FiniteElementData::get_first_quad_index() const +FiniteElementData::get_first_quad_index(const unsigned int quad_no) const { - return first_quad_index; + if (first_index_of_quads.size() == 1) + return first_index_of_quads[0] + quad_no * n_dofs_per_quad(0); + else + return first_index_of_quads[quad_no]; } template @@ -757,16 +921,22 @@ FiniteElementData::get_first_hex_index() const template unsigned int -FiniteElementData::get_first_face_line_index() const +FiniteElementData::get_first_face_line_index( + const unsigned int face_no) const { - return first_face_line_index; + return first_line_index_of_faces[first_line_index_of_faces.size() == 1 ? + 0 : + face_no]; } template unsigned int -FiniteElementData::get_first_face_quad_index() const +FiniteElementData::get_first_face_quad_index( + const unsigned int face_no) const { - return first_face_quad_index; + return first_quad_index_of_faces[first_quad_index_of_faces.size() == 1 ? + 0 : + face_no]; } diff --git a/source/fe/fe_data.cc b/source/fe/fe_data.cc index 21e16f708a..497f398444 100644 --- a/source/fe/fe_data.cc +++ b/source/fe/fe_data.cc @@ -19,6 +19,107 @@ DEAL_II_NAMESPACE_OPEN +namespace +{ + internal::GenericDoFsPerObject + expand(const unsigned int dim, + const std::vector &dofs_per_object, + const ReferenceCell::Type cell_type) + { + internal::GenericDoFsPerObject result; + + const unsigned int face_no = 0; + + result.dofs_per_object_exclusive.resize(4, std::vector(1)); + result.dofs_per_object_inclusive.resize(4, std::vector(1)); + result.object_index.resize(4, std::vector(1)); + result.first_object_index_on_face.resize(3, std::vector(1)); + + // dofs_per_vertex + const unsigned int dofs_per_vertex = dofs_per_object[0]; + result.dofs_per_object_exclusive[0][0] = dofs_per_vertex; + + // dofs_per_line + const unsigned int dofs_per_line = dofs_per_object[1]; + result.dofs_per_object_exclusive[1][0] = dofs_per_line; + + // dofs_per_quad + const unsigned int dofs_per_quad = dim > 1 ? dofs_per_object[2] : 0; + result.dofs_per_object_exclusive[2][0] = dofs_per_quad; + + // dofs_per_hex + const unsigned int dofs_per_hex = dim > 2 ? dofs_per_object[3] : 0; + result.dofs_per_object_exclusive[3][0] = dofs_per_hex; + + + // first_line_index + const unsigned int first_line_index = + (ReferenceCell::internal::Info::get_cell(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 + + ReferenceCell::internal::Info::get_cell(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 ? + ReferenceCell::internal::Info::get_cell(cell_type).n_faces() : + 0)) * + dofs_per_quad); + + // first_face_line_index + result.first_object_index_on_face[1][0] = + (ReferenceCell::internal::Info::get_face(cell_type, face_no) + .n_vertices() * + dofs_per_vertex); + + // first_face_quad_index + result.first_object_index_on_face[2][0] = + ((dim == 3 ? + ReferenceCell::internal::Info::get_face(cell_type, face_no) + .n_vertices() * + dofs_per_vertex : + ReferenceCell::internal::Info::get_cell(cell_type).n_vertices() * + dofs_per_vertex) + + ReferenceCell::internal::Info::get_face(cell_type, face_no).n_lines() * + dofs_per_line); + + // dofs_per_face + result.dofs_per_object_inclusive[dim - 1][0] = + (ReferenceCell::internal::Info::get_face(cell_type, face_no) + .n_vertices() * + dofs_per_vertex + + ReferenceCell::internal::Info::get_face(cell_type, face_no).n_lines() * + dofs_per_line + + (dim == 3 ? 1 : 0) * dofs_per_quad); + + + // dofs_per_cell + result.dofs_per_object_inclusive[dim][0] = + (ReferenceCell::internal::Info::get_cell(cell_type).n_vertices() * + dofs_per_vertex + + ReferenceCell::internal::Info::get_cell(cell_type).n_lines() * + dofs_per_line + + (dim == 2 ? + 1 : + (dim == 3 ? + ReferenceCell::internal::Info::get_cell(cell_type).n_faces() : + 0)) * + dofs_per_quad + + (dim == 3 ? 1 : 0) * dofs_per_hex); + + return result; + } +} // namespace + template FiniteElementData::FiniteElementData( const std::vector &dofs_per_object, @@ -46,67 +147,53 @@ FiniteElementData::FiniteElementData( const unsigned int degree, const Conformity conformity, const BlockIndices & block_indices) + : FiniteElementData(expand(dim, dofs_per_object, cell_type), + cell_type, + n_components, + degree, + conformity, + block_indices) +{} + +template +FiniteElementData::FiniteElementData( + const internal::GenericDoFsPerObject &data, + const ReferenceCell::Type cell_type, + const unsigned int n_components, + const unsigned int degree, + const Conformity conformity, + const BlockIndices & block_indices) : cell_type(cell_type) - , dofs_per_vertex(dofs_per_object[0]) - , dofs_per_line(dofs_per_object[1]) - , dofs_per_quad(dim > 1 ? dofs_per_object[2] : 0) - , dofs_per_quad_max(dofs_per_quad) - , dofs_per_hex(dim > 2 ? dofs_per_object[3] : 0) - , first_line_index( - ReferenceCell::internal::Info::get_cell(cell_type).n_vertices() * - dofs_per_vertex) - , first_quad_index( - first_line_index + - ReferenceCell::internal::Info::get_cell(cell_type).n_lines() * - dofs_per_line) - , first_hex_index( - first_quad_index + - (dim == 2 ? - 1 : - (dim == 3 ? - ReferenceCell::internal::Info::get_cell(cell_type).n_faces() : - 0)) * - dofs_per_quad) - , first_face_line_index( - ReferenceCell::internal::Info::get_face(cell_type, 0).n_vertices() * - dofs_per_vertex) - , first_face_quad_index( - (dim == 3 ? - ReferenceCell::internal::Info::get_face(cell_type, 0).n_vertices() * - dofs_per_vertex : - ReferenceCell::internal::Info::get_cell(cell_type).n_vertices() * - dofs_per_vertex) + - ReferenceCell::internal::Info::get_face(cell_type, 0).n_lines() * - dofs_per_line) - , dofs_per_face( - ReferenceCell::internal::Info::get_face(cell_type, 0).n_vertices() * - dofs_per_vertex + - ReferenceCell::internal::Info::get_face(cell_type, 0).n_lines() * - dofs_per_line + - (dim == 3 ? 1 : 0) * dofs_per_quad) - , dofs_per_face_max(dofs_per_face) - , dofs_per_cell( - ReferenceCell::internal::Info::get_cell(cell_type).n_vertices() * - dofs_per_vertex + - ReferenceCell::internal::Info::get_cell(cell_type).n_lines() * - dofs_per_line + - (dim == 2 ? - 1 : - (dim == 3 ? - ReferenceCell::internal::Info::get_cell(cell_type).n_faces() : - 0)) * - dofs_per_quad + - (dim == 3 ? 1 : 0) * dofs_per_hex) + , number_unique_quads(data.dofs_per_object_inclusive[2].size()) + , number_unique_faces(data.dofs_per_object_inclusive[dim - 1].size()) + , dofs_per_vertex(data.dofs_per_object_exclusive[0][0]) + , dofs_per_line(data.dofs_per_object_exclusive[1][0]) + , n_dofs_on_quad(dim > 1 ? data.dofs_per_object_exclusive[2] : + std::vector{0}) + , dofs_per_quad(n_dofs_on_quad[0]) + , dofs_per_quad_max( + *max_element(n_dofs_on_quad.begin(), n_dofs_on_quad.end())) + , dofs_per_hex(dim > 2 ? data.dofs_per_object_exclusive[3][0] : 0) + , first_line_index(data.object_index[1][0]) + , first_index_of_quads(data.object_index[2]) + , first_quad_index(first_index_of_quads[0]) + , first_hex_index(data.object_index[3][0]) + , first_line_index_of_faces(data.first_object_index_on_face[1]) + , first_face_line_index(first_line_index_of_faces[0]) + , first_quad_index_of_faces(data.first_object_index_on_face[2]) + , first_face_quad_index(first_quad_index_of_faces[0]) + , n_dofs_on_face(data.dofs_per_object_inclusive[dim - 1]) + , dofs_per_face(n_dofs_on_face[0]) + , dofs_per_face_max( + *max_element(n_dofs_on_face.begin(), n_dofs_on_face.end())) + , dofs_per_cell(data.dofs_per_object_inclusive[dim][0]) , components(n_components) , degree(degree) , conforming_space(conformity) , block_indices_data(block_indices.size() == 0 ? BlockIndices(1, dofs_per_cell) : block_indices) -{ - Assert(dofs_per_object.size() == dim + 1, - ExcDimensionMismatch(dofs_per_object.size() - 1, dim)); -} +{} -- 2.39.5