From 6e83e2948fb1ad62a10afe3116677e8d9ae8bf92 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 29 Sep 2023 21:21:55 +0200 Subject: [PATCH] Reduce compile time by using more common code --- source/dofs/dof_tools_sparsity.cc | 816 +++++++++--------------------- 1 file changed, 227 insertions(+), 589 deletions(-) diff --git a/source/dofs/dof_tools_sparsity.cc b/source/dofs/dof_tools_sparsity.cc index 940d24a824..bde01c8d47 100644 --- a/source/dofs/dof_tools_sparsity.cc +++ b/source/dofs/dof_tools_sparsity.cc @@ -763,6 +763,92 @@ namespace DoFTools { namespace { + // helper function + template + void + add_cell_entries( + const Iterator &cell, + const unsigned int face_no, + const Iterator2 &neighbor, + const unsigned int neighbor_face_no, + const Table<2, Coupling> &flux_mask, + const std::vector &dofs_on_this_cell, + std::vector &dofs_on_other_cell, + std::vector> &cell_entries) + { + dofs_on_other_cell.resize(neighbor->get_fe().n_dofs_per_cell()); + neighbor->get_dof_indices(dofs_on_other_cell); + + // Keep expensive data structures in separate vectors for inner j loop + // in separate vectors + boost::container::small_vector + component_indices_neighbor(neighbor->get_fe().n_dofs_per_cell()); + boost::container::small_vector support_on_face_i( + neighbor->get_fe().n_dofs_per_cell()); + boost::container::small_vector support_on_face_e( + neighbor->get_fe().n_dofs_per_cell()); + for (unsigned int j = 0; j < neighbor->get_fe().n_dofs_per_cell(); ++j) + { + component_indices_neighbor[j] = + (neighbor->get_fe().is_primitive(j) ? + neighbor->get_fe().system_to_component_index(j).first : + neighbor->get_fe() + .get_nonzero_components(j) + .first_selected_component()); + support_on_face_i[j] = + neighbor->get_fe().has_support_on_face(j, face_no); + support_on_face_e[j] = + neighbor->get_fe().has_support_on_face(j, neighbor_face_no); + } + + // For the parallel setting, must include also the off-diagonal + // coupling, otherwise we will include those on the other cell + for (unsigned int f = 0; f < (neighbor->is_locally_owned() ? 1 : 2); + ++f) + { + const auto &fe = (f == 0) ? cell->get_fe() : neighbor->get_fe(); + const auto &dofs_i = + (f == 0) ? dofs_on_this_cell : dofs_on_other_cell; + for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) + { + const unsigned int ii = + (fe.is_primitive(i) ? + fe.system_to_component_index(i).first : + fe.get_nonzero_components(i).first_selected_component()); + + Assert(ii < fe.n_components(), ExcInternalError()); + const bool i_non_zero_i = + fe.has_support_on_face(i, + (f == 0 ? face_no : neighbor_face_no)); + + for (unsigned int j = 0; + j < neighbor->get_fe().n_dofs_per_cell(); + ++j) + { + const bool j_non_zero_e = support_on_face_e[j]; + const unsigned int jj = component_indices_neighbor[j]; + + Assert(jj < neighbor->get_fe().n_components(), + ExcInternalError()); + + if ((flux_mask(ii, jj) == always) || + (flux_mask(ii, jj) == nonzero && i_non_zero_i && + j_non_zero_e)) + cell_entries.emplace_back(dofs_i[i], + dofs_on_other_cell[j]); + if ((flux_mask(jj, ii) == always) || + (flux_mask(jj, ii) == nonzero && j_non_zero_e && + i_non_zero_i)) + cell_entries.emplace_back(dofs_on_other_cell[j], + dofs_i[i]); + } + } + } + } + + + // implementation of the same function in namespace DoFTools for // non-hp-DoFHandlers template @@ -784,607 +870,159 @@ namespace DoFTools SparsityPatternBase::size_type>> cell_entries; - if (dof.has_hp_capabilities() == false) + const dealii::hp::FECollection &fe = + dof.get_fe_collection(); + + std::vector dofs_on_this_cell( + dof.get_fe_collection().max_dofs_per_cell()); + std::vector dofs_on_other_cell( + dof.get_fe_collection().max_dofs_per_cell()); + + const unsigned int n_components = fe.n_components(); + AssertDimension(int_mask.size(0), n_components); + AssertDimension(int_mask.size(1), n_components); + AssertDimension(flux_mask.size(0), n_components); + AssertDimension(flux_mask.size(1), n_components); + + // note that we also need to set the respective entries if flux_mask + // says so. this is necessary since we need to consider all degrees + // of freedom on a cell for interior faces. + Table<2, Coupling> int_and_flux_mask(n_components, n_components); + for (unsigned int c1 = 0; c1 < n_components; ++c1) + for (unsigned int c2 = 0; c2 < n_components; ++c2) + if (int_mask(c1, c2) != none || flux_mask(c1, c2) != none) + int_and_flux_mask(c1, c2) = always; + + // Convert the int_dof_mask to bool_int_dof_mask so we can pass it + // to constraints.add_entries_local_to_global() + std::vector> int_and_flux_dof_mask = + dof_couplings_from_component_couplings(fe, int_and_flux_mask); + + // Convert int_and_flux_dof_mask to bool_int_and_flux_dof_mask so we + // can pass it to constraints.add_entries_local_to_global() + std::vector> bool_int_and_flux_dof_mask(fe.size()); + for (unsigned int f = 0; f < fe.size(); ++f) { - const FiniteElement &fe = dof.get_fe(); - - std::vector dofs_on_this_cell( - fe.n_dofs_per_cell()); - std::vector dofs_on_other_cell( - fe.n_dofs_per_cell()); - - const Table<2, Coupling> - int_dof_mask = - dof_couplings_from_component_couplings(fe, int_mask), - flux_dof_mask = - dof_couplings_from_component_couplings(fe, flux_mask); - - Table<2, bool> support_on_face(fe.n_dofs_per_cell(), - GeometryInfo::faces_per_cell); - for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) - for (const unsigned int f : GeometryInfo::face_indices()) - support_on_face(i, f) = fe.has_support_on_face(i, f); - - // Convert the int_dof_mask to bool_int_dof_mask so we can pass it - // to constraints.add_entries_local_to_global() - Table<2, bool> bool_int_dof_mask(fe.n_dofs_per_cell(), - fe.n_dofs_per_cell()); - bool_int_dof_mask.fill(false); - for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) - for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) - if (int_dof_mask(i, j) != none) - bool_int_dof_mask(i, j) = true; - - for (const auto &cell : dof.active_cell_iterators()) - if (((subdomain_id == numbers::invalid_subdomain_id) || - (subdomain_id == cell->subdomain_id())) && - cell->is_locally_owned()) - { - cell->get_dof_indices(dofs_on_this_cell); - // make sparsity pattern for this cell - constraints.add_entries_local_to_global(dofs_on_this_cell, - sparsity, - keep_constrained_dofs, - bool_int_dof_mask); - // Loop over all interior neighbors - for (const unsigned int face_n : cell->face_indices()) - { - const typename DoFHandler::face_iterator - cell_face = cell->face(face_n); - - const bool periodic_neighbor = - cell->has_periodic_neighbor(face_n); - - if (cell->at_boundary(face_n) && (!periodic_neighbor)) - { - for (unsigned int i = 0; i < fe.n_dofs_per_cell(); - ++i) - { - const bool i_non_zero_i = - support_on_face(i, face_n); - for (unsigned int j = 0; j < fe.n_dofs_per_cell(); - ++j) - { - const bool j_non_zero_i = - support_on_face(j, face_n); - - if (flux_dof_mask(i, j) == always || - (flux_dof_mask(i, j) == nonzero && - i_non_zero_i && j_non_zero_i)) - cell_entries.emplace_back( - dofs_on_this_cell[i], - dofs_on_this_cell[j]); - } - } - sparsity.add_entries(make_array_view(cell_entries)); - cell_entries.clear(); - } - else - { - typename DoFHandler::level_cell_iterator - neighbor = - cell->neighbor_or_periodic_neighbor(face_n); - // If the cells are on the same level (and both are - // active, locally-owned cells) then only add to the - // sparsity pattern if the current cell is 'greater' - // in the total ordering. - if (neighbor->level() == cell->level() && - neighbor->index() > cell->index() && - neighbor->is_active() && - neighbor->is_locally_owned()) - continue; - // If we are more refined then the neighbor, then we - // will automatically find the active neighbor cell - // when we call 'neighbor (face_n)' above. The - // opposite is not true; if the neighbor is more - // refined then the call 'neighbor (face_n)' will - // *not* return an active cell. Hence, only add things - // to the sparsity pattern if (when the levels are - // different) the neighbor is coarser than the current - // cell. - // - // Like above, do not use this optimization if the - // neighbor is not locally owned. - if (neighbor->level() != cell->level() && - ((!periodic_neighbor && - !cell->neighbor_is_coarser(face_n)) || - (periodic_neighbor && - !cell->periodic_neighbor_is_coarser(face_n))) && - neighbor->is_locally_owned()) - continue; // (the neighbor is finer) - - if (!face_has_flux_coupling(cell, face_n)) - continue; - - - const unsigned int neighbor_face_n = - periodic_neighbor ? - cell->periodic_neighbor_face_no(face_n) : - cell->neighbor_face_no(face_n); - - - // In 1d, go straight to the cell behind this - // particular cell's most terminal cell. This makes us - // skip the if (neighbor->has_children()) section - // below. We need to do this since we otherwise - // iterate over the children of the face, which are - // always 0 in 1d. - if (dim == 1) - while (neighbor->has_children()) - neighbor = neighbor->child(face_n == 0 ? 1 : 0); - - if (neighbor->has_children()) - { - for (unsigned int sub_nr = 0; - sub_nr != cell_face->n_children(); - ++sub_nr) - { - const typename DoFHandler:: - level_cell_iterator sub_neighbor = - periodic_neighbor ? - cell - ->periodic_neighbor_child_on_subface( - face_n, sub_nr) : - cell->neighbor_child_on_subface(face_n, - sub_nr); - - sub_neighbor->get_dof_indices( - dofs_on_other_cell); - for (unsigned int i = 0; - i < fe.n_dofs_per_cell(); - ++i) - { - const bool i_non_zero_i = - support_on_face(i, face_n); - const bool i_non_zero_e = - support_on_face(i, neighbor_face_n); - for (unsigned int j = 0; - j < fe.n_dofs_per_cell(); - ++j) - { - const bool j_non_zero_i = - support_on_face(j, face_n); - const bool j_non_zero_e = - support_on_face(j, neighbor_face_n); - - if (flux_dof_mask(i, j) == always) - { - cell_entries.emplace_back( - dofs_on_this_cell[i], - dofs_on_other_cell[j]); - cell_entries.emplace_back( - dofs_on_other_cell[i], - dofs_on_this_cell[j]); - cell_entries.emplace_back( - dofs_on_this_cell[i], - dofs_on_this_cell[j]); - cell_entries.emplace_back( - dofs_on_other_cell[i], - dofs_on_other_cell[j]); - } - else if (flux_dof_mask(i, j) == - nonzero) - { - if (i_non_zero_i && j_non_zero_e) - cell_entries.emplace_back( - dofs_on_this_cell[i], - dofs_on_other_cell[j]); - if (i_non_zero_e && j_non_zero_i) - cell_entries.emplace_back( - dofs_on_other_cell[i], - dofs_on_this_cell[j]); - if (i_non_zero_i && j_non_zero_i) - cell_entries.emplace_back( - dofs_on_this_cell[i], - dofs_on_this_cell[j]); - if (i_non_zero_e && j_non_zero_e) - cell_entries.emplace_back( - dofs_on_other_cell[i], - dofs_on_other_cell[j]); - } - - if (flux_dof_mask(j, i) == always) - { - cell_entries.emplace_back( - dofs_on_this_cell[j], - dofs_on_other_cell[i]); - cell_entries.emplace_back( - dofs_on_other_cell[j], - dofs_on_this_cell[i]); - cell_entries.emplace_back( - dofs_on_this_cell[j], - dofs_on_this_cell[i]); - cell_entries.emplace_back( - dofs_on_other_cell[j], - dofs_on_other_cell[i]); - } - else if (flux_dof_mask(j, i) == - nonzero) - { - if (j_non_zero_i && i_non_zero_e) - cell_entries.emplace_back( - dofs_on_this_cell[j], - dofs_on_other_cell[i]); - if (j_non_zero_e && i_non_zero_i) - cell_entries.emplace_back( - dofs_on_other_cell[j], - dofs_on_this_cell[i]); - if (j_non_zero_i && i_non_zero_i) - cell_entries.emplace_back( - dofs_on_this_cell[j], - dofs_on_this_cell[i]); - if (j_non_zero_e && i_non_zero_e) - cell_entries.emplace_back( - dofs_on_other_cell[j], - dofs_on_other_cell[i]); - } - } - } - } - sparsity.add_entries( - make_array_view(cell_entries)); - cell_entries.clear(); - } - else - { - neighbor->get_dof_indices(dofs_on_other_cell); - for (unsigned int i = 0; i < fe.n_dofs_per_cell(); - ++i) - { - const bool i_non_zero_i = - support_on_face(i, face_n); - const bool i_non_zero_e = - support_on_face(i, neighbor_face_n); - for (unsigned int j = 0; - j < fe.n_dofs_per_cell(); - ++j) - { - const bool j_non_zero_i = - support_on_face(j, face_n); - const bool j_non_zero_e = - support_on_face(j, neighbor_face_n); - if (flux_dof_mask(i, j) == always) - { - cell_entries.emplace_back( - dofs_on_this_cell[i], - dofs_on_other_cell[j]); - cell_entries.emplace_back( - dofs_on_other_cell[i], - dofs_on_this_cell[j]); - cell_entries.emplace_back( - dofs_on_this_cell[i], - dofs_on_this_cell[j]); - cell_entries.emplace_back( - dofs_on_other_cell[i], - dofs_on_other_cell[j]); - } - if (flux_dof_mask(i, j) == nonzero) - { - if (i_non_zero_i && j_non_zero_e) - cell_entries.emplace_back( - dofs_on_this_cell[i], - dofs_on_other_cell[j]); - if (i_non_zero_e && j_non_zero_i) - cell_entries.emplace_back( - dofs_on_other_cell[i], - dofs_on_this_cell[j]); - if (i_non_zero_i && j_non_zero_i) - cell_entries.emplace_back( - dofs_on_this_cell[i], - dofs_on_this_cell[j]); - if (i_non_zero_e && j_non_zero_e) - cell_entries.emplace_back( - dofs_on_other_cell[i], - dofs_on_other_cell[j]); - } - - if (flux_dof_mask(j, i) == always) - { - cell_entries.emplace_back( - dofs_on_this_cell[j], - dofs_on_other_cell[i]); - cell_entries.emplace_back( - dofs_on_other_cell[j], - dofs_on_this_cell[i]); - cell_entries.emplace_back( - dofs_on_this_cell[j], - dofs_on_this_cell[i]); - cell_entries.emplace_back( - dofs_on_other_cell[j], - dofs_on_other_cell[i]); - } - if (flux_dof_mask(j, i) == nonzero) - { - if (j_non_zero_i && i_non_zero_e) - cell_entries.emplace_back( - dofs_on_this_cell[j], - dofs_on_other_cell[i]); - if (j_non_zero_e && i_non_zero_i) - cell_entries.emplace_back( - dofs_on_other_cell[j], - dofs_on_this_cell[i]); - if (j_non_zero_i && i_non_zero_i) - cell_entries.emplace_back( - dofs_on_this_cell[j], - dofs_on_this_cell[i]); - if (j_non_zero_e && i_non_zero_e) - cell_entries.emplace_back( - dofs_on_other_cell[j], - dofs_on_other_cell[i]); - } - } - } - sparsity.add_entries( - make_array_view(cell_entries)); - cell_entries.clear(); - } - } - } - } + bool_int_and_flux_dof_mask[f].reinit( + TableIndices<2>(fe[f].n_dofs_per_cell(), + fe[f].n_dofs_per_cell())); + bool_int_and_flux_dof_mask[f].fill(false); + for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i) + for (unsigned int j = 0; j < fe[f].n_dofs_per_cell(); ++j) + if (int_and_flux_dof_mask[f](i, j) != none) + bool_int_and_flux_dof_mask[f](i, j) = true; } - else - { - // while the implementation above is quite optimized and caches a - // lot of data (see e.g. the int/flux_dof_mask tables), this is no - // longer practical for the hp-version since we would have to have - // it for all combinations of elements in the hp::FECollection. - // consequently, the implementation here is simpler and probably - // less efficient but at least readable... - - const dealii::hp::FECollection &fe = - dof.get_fe_collection(); - - std::vector dofs_on_this_cell( - dof.get_fe_collection().max_dofs_per_cell()); - std::vector dofs_on_other_cell( - dof.get_fe_collection().max_dofs_per_cell()); - - const unsigned int n_components = fe.n_components(); - AssertDimension(int_mask.size(0), n_components); - AssertDimension(int_mask.size(1), n_components); - AssertDimension(flux_mask.size(0), n_components); - AssertDimension(flux_mask.size(1), n_components); - - // note that we also need to set the respective entries if flux_mask - // says so. this is necessary since we need to consider all degrees - // of freedom on a cell for interior faces. - Table<2, Coupling> int_and_flux_mask(n_components, n_components); - for (unsigned int c1 = 0; c1 < n_components; ++c1) - for (unsigned int c2 = 0; c2 < n_components; ++c2) - if (int_mask(c1, c2) != none || flux_mask(c1, c2) != none) - int_and_flux_mask(c1, c2) = always; - - std::vector> int_and_flux_dof_mask = - dof_couplings_from_component_couplings(fe, int_and_flux_mask); - - // Convert int_and_flux_dof_mask to bool_int_and_flux_dof_mask so we - // can pass it to constraints.add_entries_local_to_global() - std::vector> bool_int_and_flux_dof_mask(fe.size()); - for (unsigned int f = 0; f < fe.size(); ++f) - { - bool_int_and_flux_dof_mask[f].reinit( - TableIndices<2>(fe[f].n_dofs_per_cell(), - fe[f].n_dofs_per_cell())); - bool_int_and_flux_dof_mask[f].fill(false); - for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i) - for (unsigned int j = 0; j < fe[f].n_dofs_per_cell(); ++j) - if (int_and_flux_dof_mask[f](i, j) != none) - bool_int_and_flux_dof_mask[f](i, j) = true; - } - for (const auto &cell : dof.active_cell_iterators()) - if (((subdomain_id == numbers::invalid_subdomain_id) || - (subdomain_id == cell->subdomain_id())) && - cell->is_locally_owned()) + for (const auto &cell : dof.active_cell_iterators()) + if (((subdomain_id == numbers::invalid_subdomain_id) || + (subdomain_id == cell->subdomain_id())) && + cell->is_locally_owned()) + { + dofs_on_this_cell.resize(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(dofs_on_this_cell); + + // make sparsity pattern for this cell also taking into + // account the couplings due to face contributions on the same + // cell + constraints.add_entries_local_to_global( + dofs_on_this_cell, + sparsity, + keep_constrained_dofs, + bool_int_and_flux_dof_mask[cell->active_fe_index()]); + + // Loop over interior faces + for (const unsigned int face : cell->face_indices()) { - dofs_on_this_cell.resize(cell->get_fe().n_dofs_per_cell()); - cell->get_dof_indices(dofs_on_this_cell); - - // make sparsity pattern for this cell also taking into - // account the couplings due to face contributions on the same - // cell - constraints.add_entries_local_to_global( - dofs_on_this_cell, - sparsity, - keep_constrained_dofs, - bool_int_and_flux_dof_mask[cell->active_fe_index()]); - - // Loop over interior faces - for (const unsigned int face : cell->face_indices()) - { - const typename DoFHandler::face_iterator - cell_face = cell->face(face); + const bool periodic_neighbor = + cell->has_periodic_neighbor(face); - const bool periodic_neighbor = - cell->has_periodic_neighbor(face); - - if ((!cell->at_boundary(face)) || periodic_neighbor) + if ((!cell->at_boundary(face)) || periodic_neighbor) + { + typename DoFHandler::level_cell_iterator + neighbor = cell->neighbor_or_periodic_neighbor(face); + + // If the cells are on the same level (and both are + // active, locally-owned cells) then only add to the + // sparsity pattern if the current cell is 'greater' in + // the total ordering. + if (neighbor->level() == cell->level() && + neighbor->index() > cell->index() && + neighbor->is_active() && neighbor->is_locally_owned()) + continue; + + // If we are more refined then the neighbor, then we + // will automatically find the active neighbor cell when + // we call 'neighbor (face)' above. The opposite is not + // true; if the neighbor is more refined then the call + // 'neighbor (face)' will *not* return an active + // cell. Hence, only add things to the sparsity pattern + // if (when the levels are different) the neighbor is + // coarser than the current cell, except in the case + // when the neighbor is not locally owned. + if (neighbor->level() != cell->level() && + ((!periodic_neighbor && + !cell->neighbor_is_coarser(face)) || + (periodic_neighbor && + !cell->periodic_neighbor_is_coarser(face))) && + neighbor->is_locally_owned()) + continue; // (the neighbor is finer) + + if (!face_has_flux_coupling(cell, face)) + continue; + + const unsigned int neighbor_face_no = + periodic_neighbor ? + cell->periodic_neighbor_face_no(face) : + cell->neighbor_face_no(face); + + // In 1d, go straight to the cell behind this + // particular cell's most terminal cell. This makes us + // skip the if (neighbor->has_children()) section + // below. We need to do this since we otherwise + // iterate over the children of the face, which are + // always 0 in 1d. + if (dim == 1) + while (neighbor->has_children()) + neighbor = neighbor->child(face == 0 ? 1 : 0); + + if (neighbor->has_children()) { - typename DoFHandler::level_cell_iterator - neighbor = - cell->neighbor_or_periodic_neighbor(face); - - // Like the non-hp-case: If the cells are on the same - // level (and both are active, locally-owned cells) - // then only add to the sparsity pattern if the - // current cell is 'greater' in the total ordering. - if (neighbor->level() == cell->level() && - neighbor->index() > cell->index() && - neighbor->is_active() && - neighbor->is_locally_owned()) - continue; - // Again, like the non-hp-case: If we are more refined - // then the neighbor, then we will automatically find - // the active neighbor cell when we call 'neighbor - // (face)' above. The opposite is not true; if the - // neighbor is more refined then the call 'neighbor - // (face)' will *not* return an active cell. Hence, - // only add things to the sparsity pattern if (when - // the levels are different) the neighbor is coarser - // than the current cell. - // - // Like above, do not use this optimization if the - // neighbor is not locally owned. - if (neighbor->level() != cell->level() && - ((!periodic_neighbor && - !cell->neighbor_is_coarser(face)) || - (periodic_neighbor && - !cell->periodic_neighbor_is_coarser(face))) && - neighbor->is_locally_owned()) - continue; // (the neighbor is finer) - - - if (!face_has_flux_coupling(cell, face)) - continue; - - // In 1d, go straight to the cell behind this - // particular cell's most terminal cell. This makes us - // skip the if (neighbor->has_children()) section - // below. We need to do this since we otherwise - // iterate over the children of the face, which are - // always 0 in 1d. - if (dim == 1) - while (neighbor->has_children()) - neighbor = neighbor->child(face == 0 ? 1 : 0); - - if (neighbor->has_children()) + for (unsigned int sub_nr = 0; + sub_nr != cell->face(face)->n_children(); + ++sub_nr) { - for (unsigned int sub_nr = 0; - sub_nr != cell_face->n_children(); - ++sub_nr) - { - const typename DoFHandler:: - level_cell_iterator sub_neighbor = - periodic_neighbor ? - cell - ->periodic_neighbor_child_on_subface( - face, sub_nr) : - cell->neighbor_child_on_subface(face, - sub_nr); - - dofs_on_other_cell.resize( - sub_neighbor->get_fe().n_dofs_per_cell()); - sub_neighbor->get_dof_indices( - dofs_on_other_cell); - for (unsigned int i = 0; - i < cell->get_fe().n_dofs_per_cell(); - ++i) - { - const unsigned int ii = - (cell->get_fe().is_primitive(i) ? - cell->get_fe() - .system_to_component_index(i) - .first : - cell->get_fe() - .get_nonzero_components(i) - .first_selected_component()); - - Assert(ii < cell->get_fe().n_components(), - ExcInternalError()); - - for (unsigned int j = 0; - j < sub_neighbor->get_fe() - .n_dofs_per_cell(); - ++j) - { - const unsigned int jj = - (sub_neighbor->get_fe() - .is_primitive(j) ? - sub_neighbor->get_fe() - .system_to_component_index(j) - .first : - sub_neighbor->get_fe() - .get_nonzero_components(j) - .first_selected_component()); - - Assert(jj < sub_neighbor->get_fe() - .n_components(), - ExcInternalError()); - - if ((flux_mask(ii, jj) == always) || - (flux_mask(ii, jj) == nonzero)) - { - cell_entries.emplace_back( - dofs_on_this_cell[i], - dofs_on_other_cell[j]); - } - - if ((flux_mask(jj, ii) == always) || - (flux_mask(jj, ii) == nonzero)) - { - cell_entries.emplace_back( - dofs_on_other_cell[j], - dofs_on_this_cell[i]); - } - } - } - } - } - else - { - dofs_on_other_cell.resize( - neighbor->get_fe().n_dofs_per_cell()); - neighbor->get_dof_indices(dofs_on_other_cell); - for (unsigned int i = 0; - i < cell->get_fe().n_dofs_per_cell(); - ++i) - { - const unsigned int ii = - (cell->get_fe().is_primitive(i) ? - cell->get_fe() - .system_to_component_index(i) - .first : - cell->get_fe() - .get_nonzero_components(i) - .first_selected_component()); - - Assert(ii < cell->get_fe().n_components(), - ExcInternalError()); - - for (unsigned int j = 0; - j < neighbor->get_fe().n_dofs_per_cell(); - ++j) - { - const unsigned int jj = - (neighbor->get_fe().is_primitive(j) ? - neighbor->get_fe() - .system_to_component_index(j) - .first : - neighbor->get_fe() - .get_nonzero_components(j) - .first_selected_component()); - - Assert( - jj < neighbor->get_fe().n_components(), - ExcInternalError()); - - if ((flux_mask(ii, jj) == always) || - (flux_mask(ii, jj) == nonzero)) - { - cell_entries.emplace_back( - dofs_on_this_cell[i], - dofs_on_other_cell[j]); - } - - if ((flux_mask(jj, ii) == always) || - (flux_mask(jj, ii) == nonzero)) - { - cell_entries.emplace_back( - dofs_on_other_cell[j], - dofs_on_this_cell[i]); - } - } - } + const typename DoFHandler:: + level_cell_iterator sub_neighbor = + periodic_neighbor ? + cell->periodic_neighbor_child_on_subface( + face, sub_nr) : + cell->neighbor_child_on_subface(face, + sub_nr); + add_cell_entries(cell, + face, + sub_neighbor, + neighbor_face_no, + flux_mask, + dofs_on_this_cell, + dofs_on_other_cell, + cell_entries); } } + else + add_cell_entries(cell, + face, + neighbor, + neighbor_face_no, + flux_mask, + dofs_on_this_cell, + dofs_on_other_cell, + cell_entries); } - sparsity.add_entries(make_array_view(cell_entries)); - cell_entries.clear(); } - } + sparsity.add_entries(make_array_view(cell_entries)); + cell_entries.clear(); + } } } // namespace -- 2.39.5