From 5e6d645a3e8b15510aad55ec63edd1d48a6248ae Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 23 Oct 2023 17:30:54 -0600 Subject: [PATCH] Simplify code logic. --- source/fe/fe_system.cc | 143 ++++++++++++++++++----------------------- 1 file changed, 63 insertions(+), 80 deletions(-) diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index da84c7809f..4df28e839c 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -905,9 +905,6 @@ FESystem::get_restriction_matrix( this->n_dofs_per_cell()) return this->restriction[refinement_case - 1][child]; - // Check if some of the matrices of the base elements are void. - bool do_restriction = true; - // shortcut for accessing local restrictions further down std::vector *> base_matrices( this->n_base_elements()); @@ -916,57 +913,50 @@ FESystem::get_restriction_matrix( { base_matrices[i] = &base_element(i).get_restriction_matrix(child, refinement_case); - if (base_matrices[i]->n() != base_element(i).n_dofs_per_cell()) - do_restriction = false; - } - Assert(do_restriction, - (typename FiniteElement::ExcProjectionVoid())); - // if we did not encounter void matrices, initialize the matrix sizes - if (do_restriction) - { - FullMatrix restriction(this->n_dofs_per_cell(), - this->n_dofs_per_cell()); + Assert(base_matrices[i]->n() == base_element(i).n_dofs_per_cell(), + (typename FiniteElement::ExcProjectionVoid())); + } - // distribute the matrices of the base finite elements to the - // matrices of this object. for this, loop over all degrees of - // freedom and take the respective entry of the underlying base - // element. - // - // note that we by definition of a base element, they are - // independent, i.e. do not couple. only DoFs that belong to the - // same instance of a base element may couple - for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i) - for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j) - { - // first find out to which base element indices i and j - // belong, and which instance thereof in case the base element - // has a multiplicity greater than one. if they should not - // happen to belong to the same instance of a base element, - // then they cannot couple, so go on with the next index - if (this->system_to_base_table[i].first != - this->system_to_base_table[j].first) - continue; - - // so get the common base element and the indices therein: - const unsigned int base = - this->system_to_base_table[i].first.first; - - const unsigned int base_index_i = - this->system_to_base_table[i].second, - base_index_j = - this->system_to_base_table[j].second; - - // if we are sure that DoFs i and j may couple, then copy - // entries of the matrices: - restriction(i, j) = - (*base_matrices[base])(base_index_i, base_index_j); - } + FullMatrix restriction(this->n_dofs_per_cell(), + this->n_dofs_per_cell()); + + // distribute the matrices of the base finite elements to the + // matrices of this object. for this, loop over all degrees of + // freedom and take the respective entry of the underlying base + // element. + // + // note that we by definition of a base element, they are + // independent, i.e. do not couple. only DoFs that belong to the + // same instance of a base element may couple + for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i) + for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j) + { + // first find out to which base element indices i and j + // belong, and which instance thereof in case the base element + // has a multiplicity greater than one. if they should not + // happen to belong to the same instance of a base element, + // then they cannot couple, so go on with the next index + if (this->system_to_base_table[i].first != + this->system_to_base_table[j].first) + continue; + + // so get the common base element and the indices therein: + const unsigned int base = this->system_to_base_table[i].first.first; + + const unsigned int base_index_i = + this->system_to_base_table[i].second, + base_index_j = + this->system_to_base_table[j].second; + + // if we are sure that DoFs i and j may couple, then copy + // entries of the matrices: + restriction(i, j) = + (*base_matrices[base])(base_index_i, base_index_j); + } - const_cast &>( - this->restriction[refinement_case - 1][child]) = - std::move(restriction); - } + const_cast &>( + this->restriction[refinement_case - 1][child]) = std::move(restriction); } return this->restriction[refinement_case - 1][child]; @@ -997,45 +987,38 @@ FESystem::get_prolongation_matrix( this->n_dofs_per_cell()) return this->prolongation[refinement_case - 1][child]; - bool do_prolongation = true; std::vector *> base_matrices( this->n_base_elements()); for (unsigned int i = 0; i < this->n_base_elements(); ++i) { base_matrices[i] = &base_element(i).get_prolongation_matrix(child, refinement_case); - if (base_matrices[i]->n() != base_element(i).n_dofs_per_cell()) - do_prolongation = false; + + Assert(base_matrices[i]->n() == base_element(i).n_dofs_per_cell(), + (typename FiniteElement::ExcEmbeddingVoid())); } - Assert(do_prolongation, - (typename FiniteElement::ExcEmbeddingVoid())); - if (do_prolongation) - { - FullMatrix prolongate(this->n_dofs_per_cell(), - this->n_dofs_per_cell()); + FullMatrix prolongate(this->n_dofs_per_cell(), + this->n_dofs_per_cell()); - for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i) - for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j) - { - if (this->system_to_base_table[i].first != - this->system_to_base_table[j].first) - continue; - const unsigned int base = - this->system_to_base_table[i].first.first; - - const unsigned int base_index_i = - this->system_to_base_table[i].second, - base_index_j = - this->system_to_base_table[j].second; - prolongate(i, j) = - (*base_matrices[base])(base_index_i, base_index_j); - } + for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i) + for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j) + { + if (this->system_to_base_table[i].first != + this->system_to_base_table[j].first) + continue; + const unsigned int base = this->system_to_base_table[i].first.first; + + const unsigned int base_index_i = + this->system_to_base_table[i].second, + base_index_j = + this->system_to_base_table[j].second; + prolongate(i, j) = + (*base_matrices[base])(base_index_i, base_index_j); + } - const_cast &>( - this->prolongation[refinement_case - 1][child]) = - std::move(prolongate); - } + const_cast &>( + this->prolongation[refinement_case - 1][child]) = std::move(prolongate); } return this->prolongation[refinement_case - 1][child]; -- 2.39.5