From 14a1c5b661542d1fb2ad1d93de2f5d671f8524bf Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Thu, 18 Nov 2021 22:57:20 +0100 Subject: [PATCH] Refactor HangingNodes::setup_constraints() --- .../matrix_free/cuda_matrix_free.templates.h | 2 +- include/deal.II/matrix_free/dof_info.h | 6 +- .../deal.II/matrix_free/dof_info.templates.h | 6 +- .../matrix_free/hanging_nodes_internal.h | 784 +++++++++--------- .../matrix_free/matrix_free.templates.h | 40 +- source/matrix_free/dof_info.cc | 6 +- .../matrix_free/mixed_mesh_hn_algorithm_01.cc | 165 ++++ .../mixed_mesh_hn_algorithm_01.output | 32 + 8 files changed, 641 insertions(+), 400 deletions(-) create mode 100644 tests/matrix_free/mixed_mesh_hn_algorithm_01.cc create mode 100644 tests/matrix_free/mixed_mesh_hn_algorithm_01.output diff --git a/include/deal.II/matrix_free/cuda_matrix_free.templates.h b/include/deal.II/matrix_free/cuda_matrix_free.templates.h index 4c64f02188..9efc73c550 100644 --- a/include/deal.II/matrix_free/cuda_matrix_free.templates.h +++ b/include/deal.II/matrix_free/cuda_matrix_free.templates.h @@ -387,7 +387,7 @@ namespace CUDAWrappers hanging_nodes.setup_constraints(cell, partitioner, - lexicographic_inv, + {lexicographic_inv}, lexicographic_dof_indices, cell_id_view); diff --git a/include/deal.II/matrix_free/dof_info.h b/include/deal.II/matrix_free/dof_info.h index d4b69bce4c..08178f7eab 100644 --- a/include/deal.II/matrix_free/dof_info.h +++ b/include/deal.II/matrix_free/dof_info.h @@ -222,9 +222,9 @@ namespace internal template bool process_hanging_node_constraints( - const HangingNodes & hanging_nodes, - const std::vector &lexicographic_mapping, - const unsigned int cell_number, + const HangingNodes & hanging_nodes, + const std::vector> &lexicographic_mapping, + const unsigned int cell_number, const TriaIterator> &cell, std::vector & dof_indices); diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h index ffb59ed4c3..814988d4d7 100644 --- a/include/deal.II/matrix_free/dof_info.templates.h +++ b/include/deal.II/matrix_free/dof_info.templates.h @@ -277,9 +277,9 @@ namespace internal template bool DoFInfo::process_hanging_node_constraints( - const HangingNodes & hanging_nodes, - const std::vector &lexicographic_mapping, - const unsigned int cell_number, + const HangingNodes & hanging_nodes, + const std::vector> &lexicographic_mapping, + const unsigned int cell_number, const TriaIterator> &cell, std::vector & dof_indices) { diff --git a/include/deal.II/matrix_free/hanging_nodes_internal.h b/include/deal.II/matrix_free/hanging_nodes_internal.h index a5bb9c42eb..2c3313292b 100644 --- a/include/deal.II/matrix_free/hanging_nodes_internal.h +++ b/include/deal.II/matrix_free/hanging_nodes_internal.h @@ -18,6 +18,7 @@ #include +#include #include #include @@ -208,9 +209,38 @@ namespace internal setup_constraints( const CellIterator & cell, const std::shared_ptr &partitioner, - const std::vector & lexicographic_mapping, - std::vector &dof_indices, - const ArrayView & mask) const; + const std::vector> &lexicographic_mapping, + std::vector & dof_indices, + const ArrayView & mask) const; + + /** + * Compute the supported components of all entries of the given + * hp::FECollection object. + */ + std::vector> + compute_supported_components( + const dealii::hp::FECollection &fe) const; + + /** + * Determine the refinement configuration of the given cell. + */ + template + ConstraintKinds + compute_refinement_configuration(const CellIterator &cell) const; + + /** + * Update the DoF indices of a given cell for the given refinement + * configuration and for the given components. + */ + template + void + update_dof_indices( + const CellIterator & cell, + const std::shared_ptr &partitioner, + const std::vector> &lexicographic_mapping, + const std::vector> & component_mask, + const ConstraintKinds & refinement_configuration, + std::vector & dof_indices) const; private: /** @@ -242,6 +272,11 @@ namespace internal std::vector::cell_iterator, unsigned int>>> line_to_cells; + + const dealii::ndarray local_lines = { + {{{{{7, 3}}, {{6, 2}}}}, + {{{{5, 1}}, {{4, 0}}}}, + {{{{11, 9}}, {{10, 8}}}}}}; }; @@ -345,24 +380,142 @@ namespace internal + template + inline std::vector> + HangingNodes::compute_supported_components( + const dealii::hp::FECollection &fe_collection) const + { + std::vector> supported_components( + fe_collection.size(), + std::vector(fe_collection.n_components(), false)); + + for (unsigned int i = 0; i < fe_collection.size(); ++i) + { + for (unsigned int base_element_index = 0, comp = 0; + base_element_index < fe_collection[i].n_base_elements(); + ++base_element_index) + for (unsigned int c = 0; + c < fe_collection[i].element_multiplicity(base_element_index); + ++c, ++comp) + if (dim == 1 || dynamic_cast *>( + &fe_collection[i].base_element( + base_element_index)) == nullptr) + supported_components[i][comp] = false; + else + supported_components[i][comp] = true; + } + + return supported_components; + } + + + template template - inline bool - HangingNodes::setup_constraints( + inline ConstraintKinds + HangingNodes::compute_refinement_configuration( + const CellIterator &cell) const + { + // TODO: for simplex or mixed meshes: nothing to do + if ((dim == 3 && line_to_cells.size() == 0) || + (cell->reference_cell().is_hyper_cube() == false)) + return ConstraintKinds::unconstrained; + + if (cell->level() == 0) + return ConstraintKinds::unconstrained; + + const std::uint16_t subcell = + cell->parent()->child_iterator_to_index(cell); + const std::uint16_t subcell_x = (subcell >> 0) & 1; + const std::uint16_t subcell_y = (subcell >> 1) & 1; + const std::uint16_t subcell_z = (subcell >> 2) & 1; + + std::uint16_t face = 0; + std::uint16_t edge = 0; + + for (unsigned int direction = 0; direction < dim; ++direction) + { + const auto side = (subcell >> direction) & 1U; + const auto face_no = direction * 2 + side; + + // ignore if at boundary + if (cell->at_boundary(face_no)) + continue; + + const auto &neighbor = cell->neighbor(face_no); + + // ignore neighbors that are artificial or have the same level or + // have children + if (neighbor->has_children() || neighbor->is_artificial() || + neighbor->level() == cell->level()) + continue; + + face |= 1 << direction; + } + + if (dim == 3) + for (unsigned int direction = 0; direction < dim; ++direction) + if (face == 0 || face == (1 << direction)) + { + const unsigned int line_no = + direction == 0 ? + (local_lines[0][subcell_y == 0][subcell_z == 0]) : + (direction == 1 ? + (local_lines[1][subcell_x == 0][subcell_z == 0]) : + (local_lines[2][subcell_x == 0][subcell_y == 0])); + + const unsigned int line_index = cell->line(line_no)->index(); + + const auto edge_neighbor = + std::find_if(line_to_cells[line_index].begin(), + line_to_cells[line_index].end(), + [&cell](const auto &edge_neighbor) { + return edge_neighbor.first->is_artificial() == + false && + edge_neighbor.first->level() < + cell->level(); + }); + + if (edge_neighbor == line_to_cells[line_index].end()) + continue; + + edge |= 1 << direction; + } + + if ((face == 0) && (edge == 0)) + return ConstraintKinds::unconstrained; + + const std::uint16_t inverted_subcell = (subcell ^ (dim == 2 ? 3 : 7)); + + const auto refinement_configuration = static_cast( + inverted_subcell + (face << 3) + (edge << 6)); + Assert(check(refinement_configuration, dim), ExcInternalError()); + return refinement_configuration; + } + + + + template + template + inline void + HangingNodes::update_dof_indices( const CellIterator & cell, const std::shared_ptr &partitioner, - const std::vector & lexicographic_mapping, - std::vector &dof_indices, - const ArrayView & masks) const + const std::vector> &lexicographic_mapping, + const std::vector> & supported_components, + const ConstraintKinds & refinement_configuration, + std::vector & dof_indices) const { - bool cell_has_hanging_node_constraints = false; - - // for simplex or mixed meshes: nothing to do - if (dim == 3 && line_to_cells.size() == 0) - return cell_has_hanging_node_constraints; + if (std::find(supported_components[cell->active_fe_index()].begin(), + supported_components[cell->active_fe_index()].end(), + true) == + supported_components[cell->active_fe_index()].end()) + return; const auto &fe = cell->get_fe(); + AssertDimension(fe.n_unique_faces(), 1); + std::vector> component_to_system_index_face_array(fe.n_components()); @@ -373,7 +526,6 @@ namespace internal std::vector idx_offset = {0}; - for (unsigned int base_element_index = 0; base_element_index < cell->get_fe().n_base_elements(); ++base_element_index) @@ -384,386 +536,252 @@ namespace internal idx_offset.back() + cell->get_fe().base_element(base_element_index).n_dofs_per_cell()); - for (unsigned int base_element_index = 0, comp = 0; - base_element_index < cell->get_fe().n_base_elements(); - ++base_element_index) - for (unsigned int c = 0; - c < cell->get_fe().element_multiplicity(base_element_index); - ++c, ++comp) + std::vector neighbor_dofs_all(idx_offset.back()); + std::vector neighbor_dofs_all_temp( + idx_offset.back()); + + const auto get_face_idx = [](const auto n_dofs_1d, + const auto face_no, + const auto i, + const auto j) -> unsigned int { + const auto direction = face_no / 2; + const auto side = face_no % 2; + const auto offset = (side == 1) ? (n_dofs_1d - 1) : 0; + + if (dim == 2) + return (direction == 0) ? (n_dofs_1d * i + offset) : + (n_dofs_1d * offset + i); + else if (dim == 3) + switch (direction) + { + case 0: + return n_dofs_1d * n_dofs_1d * i + n_dofs_1d * j + offset; + case 1: + return n_dofs_1d * n_dofs_1d * j + n_dofs_1d * offset + i; + case 2: + return n_dofs_1d * n_dofs_1d * offset + n_dofs_1d * i + j; + default: + Assert(false, ExcNotImplemented()); + } + + Assert(false, ExcNotImplemented()); + + return 0; + }; + + const std::uint16_t kind = + static_cast(refinement_configuration); + const std::uint16_t subcell = (kind >> 0) & 7; + const std::uint16_t subcell_x = (subcell >> 0) & 1; + const std::uint16_t subcell_y = (subcell >> 1) & 1; + const std::uint16_t subcell_z = (subcell >> 2) & 1; + const std::uint16_t face = (kind >> 3) & 7; + const std::uint16_t edge = (kind >> 6) & 7; + + for (unsigned int direction = 0; direction < dim; ++direction) + if ((face >> direction) & 1U) { - auto &mask = masks[comp]; - mask = ConstraintKinds::unconstrained; + const auto side = ((subcell >> direction) & 1U) == 0; + const auto face_no = direction * 2 + side; + + // read DoFs of parent of face, ... + cell->neighbor(face_no) + ->face(cell->neighbor_face_no(face_no)) + ->get_dof_indices(neighbor_dofs_all, + cell->neighbor(face_no)->active_fe_index()); + + // ... convert the global DoFs to serial ones, and ... + if (partitioner) + for (auto &index : neighbor_dofs_all) + index = partitioner->global_to_local(index); + + for (unsigned int base_element_index = 0, comp = 0; + base_element_index < cell->get_fe().n_base_elements(); + ++base_element_index) + for (unsigned int c = 0; + c < cell->get_fe().element_multiplicity(base_element_index); + ++c, ++comp) + { + if (supported_components[cell->active_fe_index()][comp] == + false) + continue; + + const unsigned int n_dofs_1d = + cell->get_fe() + .base_element(base_element_index) + .tensor_degree() + + 1; + const unsigned int dofs_per_face = + Utilities::pow(n_dofs_1d, dim - 1); + std::vector neighbor_dofs( + dofs_per_face); + const auto lex_face_mapping = + FETools::lexicographic_to_hierarchic_numbering( + n_dofs_1d - 1); + + // ... extract the DoFs of the current component + for (unsigned int i = 0; i < dofs_per_face; ++i) + neighbor_dofs[i] = neighbor_dofs_all + [component_to_system_index_face_array[comp][i]]; + + // fix DoFs depending on orientation, flip, and rotation + if (dim == 2) + { + // TODO: for mixed meshes we need to take care of + // orientation here + Assert(cell->face_orientation(face_no), + ExcNotImplemented()); + } + else if (dim == 3) + { + int rotate = 0; // TODO + if (cell->face_rotation(face_no)) // + rotate -= 1; // + if (cell->face_flip(face_no)) // + rotate -= 2; // + + rotate_face(rotate, n_dofs_1d, neighbor_dofs); + + if (cell->face_orientation(face_no) == false) + transpose_face(n_dofs_1d - 1, neighbor_dofs); + } + else + { + Assert(false, ExcNotImplemented()); + } + + // update DoF map + for (unsigned int i = 0, k = 0; i < n_dofs_1d; ++i) + for (unsigned int j = 0; j < (dim == 2 ? 1 : n_dofs_1d); + ++j, ++k) + dof_indices[get_face_idx(n_dofs_1d, face_no, i, j) + + idx_offset[comp]] = + neighbor_dofs[lex_face_mapping[k]]; + } + } - const auto &fe_base = - cell->get_fe().base_element(base_element_index); + if (dim == 3) + for (unsigned int direction = 0; direction < dim; ++direction) + if ((edge >> direction) & 1U) + { + const unsigned int line_no = + direction == 0 ? + (local_lines[0][subcell_y][subcell_z]) : + (direction == 1 ? (local_lines[1][subcell_x][subcell_z]) : + (local_lines[2][subcell_x][subcell_y])); + + const unsigned int line_index = cell->line(line_no)->index(); + + const auto edge_neighbor = + std::find_if(line_to_cells[line_index].begin(), + line_to_cells[line_index].end(), + [&cell](const auto &edge_neighbor) { + return edge_neighbor.first->is_artificial() == + false && + edge_neighbor.first->level() < + cell->level(); + }); + + if (edge_neighbor == line_to_cells[line_index].end()) + continue; + + const auto neighbor_cell = edge_neighbor->first; + const auto local_line_neighbor = edge_neighbor->second; + + DoFCellAccessor( + &neighbor_cell->get_triangulation(), + neighbor_cell->level(), + neighbor_cell->index(), + &cell->get_dof_handler()) + .get_dof_indices(neighbor_dofs_all); + + if (partitioner) + for (auto &index : neighbor_dofs_all) + index = partitioner->global_to_local(index); + + for (unsigned int i = 0; i < neighbor_dofs_all_temp.size(); ++i) + neighbor_dofs_all_temp[i] = neighbor_dofs_all + [lexicographic_mapping[cell->active_fe_index()][i]]; + + const bool flipped = + cell->line_orientation(line_no) != + neighbor_cell->line_orientation(local_line_neighbor); + + for (unsigned int base_element_index = 0, comp = 0; + base_element_index < cell->get_fe().n_base_elements(); + ++base_element_index) + for (unsigned int c = 0; + c < + cell->get_fe().element_multiplicity(base_element_index); + ++c, ++comp) + { + if (supported_components[cell->active_fe_index()][comp] == + false) + continue; - if (dim == 1 || - dynamic_cast *>(&fe_base) == nullptr) - continue; + const unsigned int n_dofs_1d = + cell->get_fe() + .base_element(base_element_index) + .tensor_degree() + + 1; + + for (unsigned int i = 0; i < n_dofs_1d; ++i) + dof_indices[line_dof_idx(line_no, i, n_dofs_1d) + + idx_offset[comp]] = neighbor_dofs_all_temp + [line_dof_idx(local_line_neighbor, + flipped ? (n_dofs_1d - 1 - i) : i, + n_dofs_1d) + + idx_offset[comp]]; + } + } + } - const unsigned int fe_degree = fe_base.tensor_degree(); - const unsigned int n_dofs_1d = fe_degree + 1; - const unsigned int dofs_per_face = - Utilities::fixed_power(n_dofs_1d); - std::vector neighbor_dofs_all( - idx_offset.back()); - std::vector neighbor_dofs_all_temp( - idx_offset.back()); - std::vector neighbor_dofs(dofs_per_face); + template + template + inline bool + HangingNodes::setup_constraints( + const CellIterator & cell, + const std::shared_ptr &partitioner, + const std::vector> &lexicographic_mapping, + std::vector & dof_indices, + const ArrayView & masks) const + { + // 1) check if finite elements support fast hanging-node algorithm + const auto supported_components = compute_supported_components( + cell->get_dof_handler().get_fe_collection()); + + if ([](const auto &supported_components) { + return std::none_of(supported_components.begin(), + supported_components.end(), + [](const auto &a) { + return *std::max_element(a.begin(), a.end()); + }); + }(supported_components)) + return false; - const auto lex_face_mapping = - FETools::lexicographic_to_hierarchic_numbering( - fe_degree); + // 2) determine the refinement configuration of the cell + const auto refinement_configuration = + compute_refinement_configuration(cell); - for (const unsigned int face : GeometryInfo::face_indices()) - { - if ((!cell->at_boundary(face)) && - (cell->neighbor(face)->has_children() == false)) - { - const auto &neighbor = cell->neighbor(face); + if (refinement_configuration == ConstraintKinds::unconstrained) + return false; - if (neighbor->is_artificial()) - continue; + // 3) update DoF indices of cell for specified components + update_dof_indices(cell, + partitioner, + lexicographic_mapping, + supported_components, + refinement_configuration, + dof_indices); - // Neighbor is coarser than us, i.e., face is constrained - if (neighbor->level() < cell->level()) - { - const unsigned int neighbor_face = - cell->neighbor_face_no(face); - - // Find position of face on neighbor - unsigned int subface = 0; - for (; - subface < GeometryInfo::max_children_per_face; - ++subface) - if (neighbor->neighbor_child_on_subface(neighbor_face, - subface) == - cell) - break; - - // Get indices to read - DoFAccessor( - &neighbor->face(neighbor_face)->get_triangulation(), - neighbor->face(neighbor_face)->level(), - neighbor->face(neighbor_face)->index(), - &cell->get_dof_handler()) - .get_dof_indices(neighbor_dofs_all); - - for (unsigned int i = 0; i < dofs_per_face; ++i) - neighbor_dofs[i] = neighbor_dofs_all - [component_to_system_index_face_array[comp][i]]; - - // If the vector is distributed, we need to transform - // the global indices to local ones. - if (partitioner) - for (auto &index : neighbor_dofs) - index = partitioner->global_to_local(index); - - if (dim == 2) - { - if (face < 2) - { - mask |= ConstraintKinds::face_x; - if (face == 0) - mask |= ConstraintKinds::subcell_x; - if (subface == 0) - mask |= ConstraintKinds::subcell_y; - } - else - { - mask |= ConstraintKinds::face_y; - if (face == 2) - mask |= ConstraintKinds::subcell_y; - if (subface == 0) - mask |= ConstraintKinds::subcell_x; - } - - // Reorder neighbor_dofs and copy into faceth face - // of dof_indices - - // Offset if upper/right face - unsigned int offset = - (face % 2 == 1) ? fe_degree : 0; - - for (unsigned int i = 0; i < n_dofs_1d; ++i) - { - unsigned int idx = 0; - // If X-line, i.e., if y = 0 or y = fe_degree - if (face > 1) - idx = n_dofs_1d * offset + i; - // If Y-line, i.e., if x = 0 or x = fe_degree - else - idx = n_dofs_1d * i + offset; - - dof_indices[idx + idx_offset[comp]] = - neighbor_dofs[lex_face_mapping[i]]; - } - } - else if (dim == 3) - { - const bool transpose = - !(cell->face_orientation(face)); - - int rotate = 0; - - if (cell->face_rotation(face)) - rotate -= 1; - if (cell->face_flip(face)) - rotate -= 2; - - rotate_face(rotate, n_dofs_1d, neighbor_dofs); - rotate_subface_index(rotate, subface); - - if (transpose) - { - transpose_face(fe_degree, neighbor_dofs); - transpose_subface_index(subface); - } - - // YZ-plane - if (face < 2) - { - mask |= ConstraintKinds::face_x; - if (face == 0) - mask |= ConstraintKinds::subcell_x; - if (subface % 2 == 0) - mask |= ConstraintKinds::subcell_y; - if (subface / 2 == 0) - mask |= ConstraintKinds::subcell_z; - } - // XZ-plane - else if (face < 4) - { - mask |= ConstraintKinds::face_y; - if (face == 2) - mask |= ConstraintKinds::subcell_y; - if (subface % 2 == 0) - mask |= ConstraintKinds::subcell_z; - if (subface / 2 == 0) - mask |= ConstraintKinds::subcell_x; - } - // XY-plane - else - { - mask |= ConstraintKinds::face_z; - if (face == 4) - mask |= ConstraintKinds::subcell_z; - if (subface % 2 == 0) - mask |= ConstraintKinds::subcell_x; - if (subface / 2 == 0) - mask |= ConstraintKinds::subcell_y; - } - - // Offset if upper/right/back face - unsigned int offset = - (face % 2 == 1) ? fe_degree : 0; - - for (unsigned int i = 0; i < n_dofs_1d; ++i) - { - for (unsigned int j = 0; j < n_dofs_1d; ++j) - { - unsigned int idx = 0; - // If YZ-plane, i.e., if x = 0 or x = - // fe_degree, and orientation standard - if (face < 2) - idx = n_dofs_1d * n_dofs_1d * i + - n_dofs_1d * j + offset; - // If XZ-plane, i.e., if y = 0 or y = - // fe_degree, and orientation standard - else if (face < 4) - idx = n_dofs_1d * n_dofs_1d * j + - n_dofs_1d * offset + i; - // If XY-plane, i.e., if z = 0 or z = - // fe_degree, and orientation standard - else - idx = n_dofs_1d * n_dofs_1d * offset + - n_dofs_1d * i + j; - - dof_indices[idx + idx_offset[comp]] = - neighbor_dofs - [lex_face_mapping[n_dofs_1d * i + j]]; - } - } - } - else - ExcNotImplemented(); - } - } - } - - // In 3D we can have a situation where only DoFs on an edge are - // constrained. Append these here. - if (dim == 3) - { - // For each line on cell, which faces does it belong to, what is - // the edge mask, what is the types of the faces it belong to, - // and what is the type along the edge. - const ConstraintKinds line_to_edge[12][4] = { - {ConstraintKinds::face_x | ConstraintKinds::face_z, - ConstraintKinds::edge_y, - ConstraintKinds::subcell_x | ConstraintKinds::subcell_z, - ConstraintKinds::subcell_y}, - {ConstraintKinds::face_x | ConstraintKinds::face_z, - ConstraintKinds::edge_y, - ConstraintKinds::subcell_z, - ConstraintKinds::subcell_y}, - {ConstraintKinds::face_y | ConstraintKinds::face_z, - ConstraintKinds::edge_x, - ConstraintKinds::subcell_y | ConstraintKinds::subcell_z, - ConstraintKinds::subcell_x}, - {ConstraintKinds::face_y | ConstraintKinds::face_z, - ConstraintKinds::edge_x, - ConstraintKinds::subcell_z, - ConstraintKinds::subcell_x}, - {ConstraintKinds::face_x | ConstraintKinds::face_z, - ConstraintKinds::edge_y, - ConstraintKinds::subcell_x, - ConstraintKinds::subcell_y}, - {ConstraintKinds::face_x | ConstraintKinds::face_z, - ConstraintKinds::edge_y, - ConstraintKinds::unconstrained, - ConstraintKinds::subcell_y}, - {ConstraintKinds::face_y | ConstraintKinds::face_z, - ConstraintKinds::edge_x, - ConstraintKinds::subcell_y, - ConstraintKinds::subcell_x}, - {ConstraintKinds::face_y | ConstraintKinds::face_z, - ConstraintKinds::edge_x, - ConstraintKinds::unconstrained, - ConstraintKinds::subcell_x}, - {ConstraintKinds::face_x | ConstraintKinds::face_y, - ConstraintKinds::edge_z, - ConstraintKinds::subcell_x | ConstraintKinds::subcell_y, - ConstraintKinds::subcell_z}, - {ConstraintKinds::face_x | ConstraintKinds::face_y, - ConstraintKinds::edge_z, - ConstraintKinds::subcell_y, - ConstraintKinds::subcell_z}, - {ConstraintKinds::face_x | ConstraintKinds::face_y, - ConstraintKinds::edge_z, - ConstraintKinds::subcell_x, - ConstraintKinds::subcell_z}, - {ConstraintKinds::face_x | ConstraintKinds::face_y, - ConstraintKinds::edge_z, - ConstraintKinds::unconstrained, - ConstraintKinds::subcell_z}}; - - for (unsigned int local_line = 0; - local_line < GeometryInfo::lines_per_cell; - ++local_line) - { - // If we don't already have a constraint for as part of a - // face - if ((mask & line_to_edge[local_line][0]) == - ConstraintKinds::unconstrained) - { - // For each cell which share that edge - const unsigned int line = - cell->line(local_line)->index(); - for (const auto &edge_neighbor : line_to_cells[line]) - { - // If one of them is coarser than us - const auto neighbor_cell = edge_neighbor.first; - - if (neighbor_cell->is_artificial()) - continue; - - if (neighbor_cell->level() < cell->level()) - { - const unsigned int local_line_neighbor = - edge_neighbor.second; - mask |= line_to_edge[local_line][1] | - line_to_edge[local_line][2]; - - bool flipped = false; - if (cell->line(local_line)->vertex_index(0) == - neighbor_cell->line(local_line_neighbor) - ->vertex_index(0)) - { - // Assuming line directions match axes - // directions, we have an unflipped edge of - // first type - mask |= line_to_edge[local_line][3]; - } - else if (cell->line(local_line) - ->vertex_index(1) == - neighbor_cell - ->line(local_line_neighbor) - ->vertex_index(1)) - { - // We have an unflipped edge of second type - } - else if (cell->line(local_line) - ->vertex_index(1) == - neighbor_cell - ->line(local_line_neighbor) - ->vertex_index(0)) - { - // We have a flipped edge of second type - flipped = true; - } - else if (cell->line(local_line) - ->vertex_index(0) == - neighbor_cell - ->line(local_line_neighbor) - ->vertex_index(1)) - { - // We have a flipped edge of first type - mask |= line_to_edge[local_line][3]; - flipped = true; - } - else - ExcInternalError(); - - // Copy the unconstrained values - DoFCellAccessor( - &neighbor_cell->get_triangulation(), - neighbor_cell->level(), - neighbor_cell->index(), - &cell->get_dof_handler()) - .get_dof_indices(neighbor_dofs_all); - // If the vector is distributed, we need to - // transform the global indices to local ones. - if (partitioner) - for (auto &index : neighbor_dofs_all) - index = partitioner->global_to_local(index); - - for (unsigned int i = 0; - i < neighbor_dofs_all_temp.size(); - ++i) - neighbor_dofs_all_temp[i] = - neighbor_dofs_all[lexicographic_mapping[i]]; - - for (unsigned int i = 0; i < n_dofs_1d; ++i) - { - // Get local dof index along line - const unsigned int idx = - line_dof_idx(local_line, i, n_dofs_1d); - - dof_indices[idx + idx_offset[comp]] = - neighbor_dofs_all_temp - [line_dof_idx(local_line_neighbor, - flipped ? fe_degree - i : - i, - n_dofs_1d) + - idx_offset[comp]]; - } - - // Stop looping over edge neighbors - break; - } - } - } - } - } - Assert(check(mask, dim), ExcInternalError()); + // 4) TODO: copy refinement configuration to all components + for (unsigned int c = 0; c < supported_components[0].size(); ++c) + if (supported_components[cell->active_fe_index()][c]) + masks[c] = refinement_configuration; - cell_has_hanging_node_constraints |= - mask != ConstraintKinds::unconstrained; - } - return cell_has_hanging_node_constraints; + return true; } diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 3a7dc80fc1..7b5c1911eb 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -1023,7 +1023,7 @@ namespace internal MatrixFreeFunctions::FaceSetup & face_setup, MatrixFreeFunctions::ConstraintValues & constraint_values, const bool use_vector_data_exchanger_full, - const bool use_fast_hanging_node_algorithm) + const bool use_fast_hanging_node_algorithm_in) { if (do_face_integrals) face_setup.initialize(dof_handler[0]->get_triangulation(), @@ -1049,6 +1049,34 @@ namespace internal bool cell_categorization_enabled = !cell_vectorization_category.empty(); + bool use_fast_hanging_node_algorithm = use_fast_hanging_node_algorithm_in; + + if (use_fast_hanging_node_algorithm) + { + const auto &reference_cells = + dof_handler[0]->get_triangulation().get_reference_cells(); + use_fast_hanging_node_algorithm = + std::all_of(reference_cells.begin(), + reference_cells.end(), + [](const auto &r) { + return r.is_hyper_cube() || r.is_simplex(); + }); + } + + if (use_fast_hanging_node_algorithm) + for (unsigned int no = 0; no < n_dof_handlers; ++no) + { + const dealii::hp::FECollection &fes = + dof_handler[no]->get_fe_collection(); + + use_fast_hanging_node_algorithm &= + std::all_of(fes.begin(), fes.end(), [&fes](const auto &fe) { + return fes[0].compare_for_domination(fe) == + FiniteElementDomination::Domination:: + either_element_can_dominate; + }); + } + for (unsigned int no = 0; no < n_dof_handlers; ++no) { const dealii::hp::FECollection &fes = @@ -1171,9 +1199,8 @@ namespace internal hanging_nodes = std::make_unique< dealii::internal::MatrixFreeFunctions::HangingNodes>(tria); for (unsigned int no = 0; no < n_dof_handlers; ++no) - if (dof_handler[no]->get_fe_collection().size() == 1) - dof_info[no].hanging_node_constraint_masks.resize( - n_active_cells * dof_handler[no]->get_fe().n_components()); + dof_info[no].hanging_node_constraint_masks.resize( + n_active_cells * dof_handler[no]->get_fe().n_components()); } for (unsigned int counter = 0; counter < n_active_cells; ++counter) @@ -1212,15 +1239,14 @@ namespace internal bool cell_has_hanging_node_constraints = false; - if (dim > 1 && use_fast_hanging_node_algorithm && - dofh->get_fe_collection().size() == 1) + if (dim > 1 && use_fast_hanging_node_algorithm) { local_dof_indices_resolved = local_dof_indices; cell_has_hanging_node_constraints = dof_info[no].process_hanging_node_constraints( *hanging_nodes, - lexicographic[no][0], + lexicographic[no], counter, cell_it, local_dof_indices_resolved); diff --git a/source/matrix_free/dof_info.cc b/source/matrix_free/dof_info.cc index 8e0eaf5a73..6bd7669f11 100644 --- a/source/matrix_free/dof_info.cc +++ b/source/matrix_free/dof_info.cc @@ -1519,21 +1519,21 @@ namespace internal template bool DoFInfo::process_hanging_node_constraints<1>( const HangingNodes<1> & hanging_nodes, - const std::vector & lexicographic_mapping, + const std::vector> & lexicographic_mapping, const unsigned int cell_number, const TriaIterator> &cell, std::vector & dof_indices); template bool DoFInfo::process_hanging_node_constraints<2>( const HangingNodes<2> & hanging_nodes, - const std::vector & lexicographic_mapping, + const std::vector> & lexicographic_mapping, const unsigned int cell_number, const TriaIterator> &cell, std::vector & dof_indices); template bool DoFInfo::process_hanging_node_constraints<3>( const HangingNodes<3> & hanging_nodes, - const std::vector & lexicographic_mapping, + const std::vector> & lexicographic_mapping, const unsigned int cell_number, const TriaIterator> &cell, std::vector & dof_indices); diff --git a/tests/matrix_free/mixed_mesh_hn_algorithm_01.cc b/tests/matrix_free/mixed_mesh_hn_algorithm_01.cc new file mode 100644 index 0000000000..61903625b0 --- /dev/null +++ b/tests/matrix_free/mixed_mesh_hn_algorithm_01.cc @@ -0,0 +1,165 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// Test that the fast matrix-free hanging-node algorithm is also working on +// adaptively refined 2D mixed meshes. + +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include "../tests.h" + + + +int +main() +{ + initlog(); + + const unsigned int dim = 2; + const unsigned int degree = 1; + + Triangulation tria_0, tria_1, tria_2, tria_3, tria; + + GridGenerator::subdivided_hyper_rectangle(tria_0, + {1, 1}, + {0.0, 0.0}, + {0.5, 0.5}); + GridGenerator::subdivided_hyper_rectangle_with_simplices(tria_1, + {1, 1}, + {0.5, 0.0}, + {1.0, 0.5}); + GridGenerator::subdivided_hyper_rectangle_with_simplices(tria_2, + {1, 1}, + {0.0, 0.5}, + {0.5, 1.0}); + GridGenerator::subdivided_hyper_rectangle_with_simplices(tria_3, + {1, 1}, + {0.5, 0.5}, + {1.0, 1.0}); + + GridGenerator::merge_triangulations({&tria_0, &tria_1, &tria_2, &tria_3}, + tria); + + auto cell = tria.begin(); + + cell->set_refine_flag(); + cell++; + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + if (false) + { + GridOut grid_out; + std::ofstream out("mesh.vtk"); + grid_out.write_vtk(tria, out); + } + + DoFHandler dof_handler(tria); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->reference_cell() == ReferenceCells::Triangle) + cell->set_active_fe_index(0); + else if (cell->reference_cell() == ReferenceCells::Quadrilateral) + cell->set_active_fe_index(1); + else + Assert(false, ExcNotImplemented()); + } + + const hp::MappingCollection<2> mapping(MappingFE<2>(FE_SimplexP<2>(1)), + MappingFE<2>(FE_Q<2>(1))); + const hp::FECollection<2> fe(FE_SimplexP<2>{degree}, FE_Q<2>{degree}); + const hp::QCollection<2> quadrature_formula(QGaussSimplex<2>(degree + 1), + QGauss<2>(degree + 1)); + + dof_handler.distribute_dofs(fe); + + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + const auto print = [](const auto &label, const auto &matrix_free) { + deallog << label << std::endl; + + for (unsigned int c = 0; + c < matrix_free.get_dof_info(0).row_starts.size() - 1; + ++c) + { + deallog + << std::setw(3) + << (matrix_free.get_dof_info(0) + .hanging_node_constraint_masks.size() == 0 ? + 0 : + static_cast( + matrix_free.get_dof_info(0).hanging_node_constraint_masks[c])) + << " : "; + + for (unsigned int i = matrix_free.get_dof_info(0).row_starts[c].first; + i < matrix_free.get_dof_info(0).row_starts[c + 1].first; + ++i) + deallog << std::setw(3) << matrix_free.get_dof_info(0).dof_indices[i] + << " "; + deallog << std::endl; + } + deallog << std::endl; + }; + + { + typename MatrixFree>::AdditionalData + additional_data; + additional_data.mapping_update_flags = update_gradients | update_values; + + MatrixFree> matrix_free; + matrix_free.reinit( + mapping, dof_handler, constraints, quadrature_formula, additional_data); + + print("use_fast_hanging_node_algorithm = true", matrix_free); + } + + { + typename MatrixFree>::AdditionalData + additional_data; + additional_data.mapping_update_flags = update_gradients | update_values; + additional_data.use_fast_hanging_node_algorithm = false; + + MatrixFree> matrix_free; + matrix_free.reinit( + mapping, dof_handler, constraints, quadrature_formula, additional_data); + + for (const auto &i : + matrix_free.get_dof_info(0).hanging_node_constraint_masks) + deallog << static_cast(i) << std::endl; + deallog << std::endl; + + print("use_fast_hanging_node_algorithm = false", matrix_free); + } +} diff --git a/tests/matrix_free/mixed_mesh_hn_algorithm_01.output b/tests/matrix_free/mixed_mesh_hn_algorithm_01.output new file mode 100644 index 0000000000..384befad90 --- /dev/null +++ b/tests/matrix_free/mixed_mesh_hn_algorithm_01.output @@ -0,0 +1,32 @@ + +DEAL::use_fast_hanging_node_algorithm = true +DEAL::0 : 12 13 14 +DEAL::0 : 13 2 2 1 +DEAL::0 : 14 2 1 1 +DEAL::0 : 13 2 1 14 +DEAL::0 : 0 1 2 +DEAL::0 : 3 1 4 +DEAL::0 : 5 4 1 +DEAL::0 : 1 0 5 +DEAL::0 : 6 5 0 +DEAL::0 : 7 8 9 10 +DEAL::0 : 8 12 10 14 +DEAL::17 : 9 10 3 1 +DEAL::16 : 10 14 3 1 +DEAL:: +DEAL:: +DEAL::use_fast_hanging_node_algorithm = false +DEAL::0 : 12 13 14 +DEAL::0 : 13 2 2 1 +DEAL::0 : 14 2 1 1 +DEAL::0 : 13 2 1 14 +DEAL::0 : 0 1 2 +DEAL::0 : 3 1 4 +DEAL::0 : 5 4 1 +DEAL::0 : 1 0 5 +DEAL::0 : 6 5 0 +DEAL::0 : 7 8 9 10 +DEAL::0 : 8 12 10 14 +DEAL::0 : 9 10 3 3 1 +DEAL::0 : 10 14 3 1 1 +DEAL:: -- 2.39.5