From 20a624fdef8862743006780886e0cc06d51300f5 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 27 May 2022 12:21:03 +0200 Subject: [PATCH] Generalize FETools::Compositing::multiply_dof_numbers() --- include/deal.II/fe/fe_base.h | 32 +++++--- include/deal.II/fe/fe_tools.templates.h | 87 +++++++++++++++------- source/fe/fe_data.cc | 26 +++++-- tests/fe/generic_dofs_per_object_01.output | 48 ++++++------ 4 files changed, 125 insertions(+), 68 deletions(-) diff --git a/include/deal.II/fe/fe_base.h b/include/deal.II/fe/fe_base.h index 84eb268511..519f2045ef 100644 --- a/include/deal.II/fe/fe_base.h +++ b/include/deal.II/fe/fe_base.h @@ -980,9 +980,9 @@ internal::GenericDoFsPerObject::generate(const FiniteElementData &fe) internal::GenericDoFsPerObject result; - result.dofs_per_object_exclusive.resize(dim + 1); - result.dofs_per_object_inclusive.resize(dim + 1); - result.object_index.resize(dim + 1); + result.dofs_per_object_exclusive.resize(4); + result.dofs_per_object_inclusive.resize(4); + result.object_index.resize(4); unsigned int counter = 0; @@ -1023,24 +1023,32 @@ internal::GenericDoFsPerObject::generate(const FiniteElementData &fe) } { - result.dofs_per_object_exclusive[dim].emplace_back( - fe.template n_dofs_per_object()); + const auto c = fe.template n_dofs_per_object(); + + result.dofs_per_object_exclusive[dim].emplace_back(c); result.dofs_per_object_inclusive[dim].emplace_back(fe.n_dofs_per_cell()); result.object_index[dim].emplace_back(counter); + + counter += c; } - result.first_object_index_on_face.resize(dim); + for (unsigned int d = dim + 1; d <= 3; ++d) + { + result.dofs_per_object_exclusive[d].emplace_back(0); + result.dofs_per_object_inclusive[d].emplace_back(0); + result.object_index[d].emplace_back(counter); + } + + result.first_object_index_on_face.resize(3); for (unsigned int face_no : reference_cell.face_indices()) { result.first_object_index_on_face[0].emplace_back(0); - if (dim >= 2) - result.first_object_index_on_face[1].emplace_back( - fe.get_first_face_line_index(face_no)); + result.first_object_index_on_face[1].emplace_back( + fe.get_first_face_line_index(face_no)); - if (dim == 3) - result.first_object_index_on_face[2].emplace_back( - fe.get_first_face_quad_index(face_no)); + result.first_object_index_on_face[2].emplace_back( + fe.get_first_face_quad_index(face_no)); } return result; diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h index c8f7678dc6..8da91c7938 100644 --- a/include/deal.II/fe/fe_tools.templates.h +++ b/include/deal.II/fe/fe_tools.templates.h @@ -85,11 +85,6 @@ namespace FETools { AssertDimension(fes.size(), multiplicities.size()); - unsigned int multiplied_dofs_per_vertex = 0; - unsigned int multiplied_dofs_per_line = 0; - unsigned int multiplied_dofs_per_quad = 0; - unsigned int multiplied_dofs_per_hex = 0; - unsigned int multiplied_n_components = 0; unsigned int degree = 0; // degree is the maximal degree of the components @@ -103,23 +98,69 @@ namespace FETools break; } + dealii::internal::GenericDoFsPerObject dpo; + + std::vector dpos_in(fes.size()); + for (unsigned int i = 0; i < fes.size(); ++i) if (multiplicities[i] > 0) - { - // TODO: the implementation makes the assumption that all faces have - // the same number of dofs -> don't construct DPO but - // PrecomputedFiniteElementData - AssertDimension(fes[i]->n_unique_quads(), 1); - - multiplied_dofs_per_vertex += - fes[i]->n_dofs_per_vertex() * multiplicities[i]; - multiplied_dofs_per_line += - fes[i]->n_dofs_per_line() * multiplicities[i]; - multiplied_dofs_per_quad += - fes[i]->n_dofs_per_quad(0) * multiplicities[i]; - multiplied_dofs_per_hex += - fes[i]->n_dofs_per_hex() * multiplicities[i]; + dpos_in[i] = + dealii::internal::GenericDoFsPerObject::generate(*fes[i]); + + // helper function to fill a vector of GenericDoFsPerObject + // according to multiplicities + const auto fill_dpo_vector = + [&](const std::function> &( + dealii::internal::GenericDoFsPerObject &)> &get_vector) { + auto &vector_dst = get_vector(dpo); + + // allocate memory + for (unsigned int i = 0; i < fes.size(); ++i) + if (multiplicities[i] > 0) + { + const auto &vector_src = get_vector(dpos_in[i]); + + vector_dst.resize(vector_src.size()); + + for (unsigned int j = 0; j < vector_src.size(); ++j) + vector_dst[j].assign(vector_src[j].size(), 0); + + break; + } + + // fill vector according to multiplicities + for (unsigned int i = 0; i < fes.size(); ++i) + if (multiplicities[i] > 0) + { + const auto &vector_src = get_vector(dpos_in[i]); + for (unsigned int j = 0; j < vector_src.size(); ++j) + for (unsigned int k = 0; k < vector_src[j].size(); ++k) + vector_dst[j][k] += vector_src[j][k] * multiplicities[i]; + } + }; + + // go through each field of GenericDoFsPerObject + fill_dpo_vector( + [](auto &dpo) -> std::vector> & { + return dpo.dofs_per_object_exclusive; + }); + fill_dpo_vector( + [](auto &dpo) -> std::vector> & { + return dpo.dofs_per_object_inclusive; + }); + fill_dpo_vector( + [](auto &dpo) -> std::vector> & { + return dpo.object_index; + }); + fill_dpo_vector( + [](auto &dpo) -> std::vector> & { + return dpo.first_object_index_on_face; + }); + + for (unsigned int i = 0; i < fes.size(); ++i) + if (multiplicities[i] > 0) + { multiplied_n_components += fes[i]->n_components() * multiplicities[i]; @@ -150,14 +191,6 @@ namespace FETools total_conformity & fes[index]->conforming_space); } - std::vector dpo; - dpo.push_back(multiplied_dofs_per_vertex); - dpo.push_back(multiplied_dofs_per_line); - if (dim > 1) - dpo.push_back(multiplied_dofs_per_quad); - if (dim > 2) - dpo.push_back(multiplied_dofs_per_hex); - BlockIndices block_indices(0, 0); for (unsigned int base = 0; base < fes.size(); ++base) diff --git a/source/fe/fe_data.cc b/source/fe/fe_data.cc index fe8d7437ed..12e5a56f75 100644 --- a/source/fe/fe_data.cc +++ b/source/fe/fe_data.cc @@ -94,6 +94,19 @@ namespace internal return result; } + + unsigned int + number_unique_entries(const std::vector &vector) + { + if (std::all_of(vector.begin(), vector.end(), [&](const auto &e) { + return e == vector.front(); + })) + { + return 1; + } + else + return vector.size(); + } } // namespace internal @@ -133,6 +146,8 @@ FiniteElementData::FiniteElementData( block_indices) {} + + template FiniteElementData::FiniteElementData( const internal::GenericDoFsPerObject &data, @@ -142,16 +157,17 @@ FiniteElementData::FiniteElementData( const Conformity conformity, const BlockIndices & block_indices) : reference_cell_kind(reference_cell) - , number_unique_quads(data.dofs_per_object_inclusive[2].size()) - , number_unique_faces(data.dofs_per_object_inclusive[dim - 1].size()) + , number_unique_quads( + internal::number_unique_entries(data.dofs_per_object_inclusive[2])) + , number_unique_faces( + internal::number_unique_entries(data.dofs_per_object_inclusive[dim - 1])) , 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}) + , n_dofs_on_quad(data.dofs_per_object_exclusive[2]) , 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) + , dofs_per_hex(data.dofs_per_object_exclusive[3][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]) diff --git a/tests/fe/generic_dofs_per_object_01.output b/tests/fe/generic_dofs_per_object_01.output index 9ad3e7cf86..22e20dbfac 100644 --- a/tests/fe/generic_dofs_per_object_01.output +++ b/tests/fe/generic_dofs_per_object_01.output @@ -1,27 +1,27 @@ DEAL::FE_Q<1>(1) -DEAL::1 1 0 -DEAL::1 1 2 -DEAL::0 1 2 -DEAL::0 0 +DEAL::1 1 0 0 0 +DEAL::1 1 2 0 0 +DEAL::0 1 2 2 2 +DEAL::0 0 1 1 2 2 DEAL:: DEAL::FE_Q<1>(2) -DEAL::1 1 1 -DEAL::1 1 3 -DEAL::0 1 2 -DEAL::0 0 +DEAL::1 1 1 0 0 +DEAL::1 1 3 0 0 +DEAL::0 1 2 3 3 +DEAL::0 0 1 1 2 2 DEAL:: DEAL::FE_Q<2>(1) -DEAL::1 1 1 1 0 0 0 0 0 -DEAL::1 1 1 1 2 2 2 2 4 -DEAL::0 1 2 3 4 4 4 4 4 -DEAL::0 0 0 0 2 2 2 2 +DEAL::1 1 1 1 0 0 0 0 0 0 +DEAL::1 1 1 1 2 2 2 2 4 0 +DEAL::0 1 2 3 4 4 4 4 4 4 +DEAL::0 0 0 0 2 2 2 2 4 4 4 4 DEAL:: DEAL::FE_Q<2>(2) -DEAL::1 1 1 1 1 1 1 1 1 -DEAL::1 1 1 1 3 3 3 3 9 -DEAL::0 1 2 3 4 5 6 7 8 -DEAL::0 0 0 0 2 2 2 2 +DEAL::1 1 1 1 1 1 1 1 1 0 +DEAL::1 1 1 1 3 3 3 3 9 0 +DEAL::0 1 2 3 4 5 6 7 8 9 +DEAL::0 0 0 0 2 2 2 2 5 5 5 5 DEAL:: DEAL::FE_Q<3>(1) DEAL::1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 @@ -36,16 +36,16 @@ DEAL::0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 DEAL::0 0 0 0 0 0 4 4 4 4 4 4 8 8 8 8 8 8 DEAL:: DEAL::FE_SimplexP<2>(1) -DEAL::1 1 1 0 0 0 0 -DEAL::1 1 1 2 2 2 3 -DEAL::0 1 2 3 3 3 3 -DEAL::0 0 0 2 2 2 +DEAL::1 1 1 0 0 0 0 0 +DEAL::1 1 1 2 2 2 3 0 +DEAL::0 1 2 3 3 3 3 3 +DEAL::0 0 0 2 2 2 3 3 3 DEAL:: DEAL::FE_SimplexP<2>(2) -DEAL::1 1 1 1 1 1 0 -DEAL::1 1 1 3 3 3 6 -DEAL::0 1 2 3 4 5 6 -DEAL::0 0 0 2 2 2 +DEAL::1 1 1 1 1 1 0 0 +DEAL::1 1 1 3 3 3 6 0 +DEAL::0 1 2 3 4 5 6 6 +DEAL::0 0 0 2 2 2 4 4 4 DEAL:: DEAL::FE_SimplexP<3>(1) DEAL::1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 -- 2.39.5