From fd4a668a47803931c2100b4467ffe3822d27aac3 Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Tue, 26 Jan 2021 20:46:07 -0700 Subject: [PATCH] Adapt CellAccessor::neighbor_child_on_subface for triangles. Required by DoFTools::make_hanging_node_constraints. --- include/deal.II/grid/reference_cell.h | 26 + source/grid/tria_accessor.cc | 727 ++++++++++-------- tests/simplex/refinement_03.cc | 100 +++ ...finement_03.with_simplex_support=on.output | 2 + 4 files changed, 524 insertions(+), 331 deletions(-) create mode 100644 tests/simplex/refinement_03.cc create mode 100644 tests/simplex/refinement_03.with_simplex_support=on.output diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index 7bf0fadc60..b830513b07 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -1009,6 +1009,21 @@ namespace ReferenceCell return 0; } + + /** + * Indices of child cells that are adjacent to a certain face of the + * mother cell. + */ + virtual unsigned int + child_cell_on_face(const unsigned int face_n, + const unsigned int subface_n) const + { + Assert(false, ExcNotImplemented()); + (void)face_n; + (void)subface_n; + + return numbers::invalid_unsigned_int; + } }; @@ -1219,6 +1234,17 @@ namespace ReferenceCell AssertIndexRange(face_n, n_faces()); return face_n; } + + virtual unsigned int + child_cell_on_face(const unsigned int face_n, + const unsigned int subface_n) const override + { + static constexpr unsigned int subcells[3][2] = {{0, 1}, + {1, 2}, + {2, 0}}; + + return subcells[face_n][subface_n]; + } }; diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index bf6f98c58f..0611d25e3d 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -2845,368 +2845,433 @@ CellAccessor::neighbor_child_on_subface( { case 2: { - const unsigned int neighbor_neighbor = - this->neighbor_of_neighbor(face); - const unsigned int neighbor_child_index = - GeometryInfo::child_cell_on_face( - this->neighbor(face)->refinement_case(), - neighbor_neighbor, - subface); - - TriaIterator> sub_neighbor = - this->neighbor(face)->child(neighbor_child_index); - // the neighbors child can have children, - // which are not further refined along the - // face under consideration. as we are - // normally interested in one of this - // child's child, search for the right one. - while (sub_neighbor->has_children()) + if (this->reference_cell_type() == ReferenceCell::Type::Tri) { - Assert((GeometryInfo::face_refinement_case( - sub_neighbor->refinement_case(), neighbor_neighbor) == - RefinementCase::no_refinement), - ExcInternalError()); - sub_neighbor = - sub_neighbor->child(GeometryInfo::child_cell_on_face( - sub_neighbor->refinement_case(), neighbor_neighbor, 0)); + const auto neighbor_cell = this->neighbor(face); + + // only for isotropic refinement at the moment + Assert(neighbor_cell->refinement_case() == + RefinementCase<2>::isotropic_refinement, + ExcNotImplemented()); + + // determine indices for this cell's subface from the perspective + // of the neighboring cell + const unsigned int neighbor_face = + this->neighbor_of_neighbor(face); + // two neighboring cells have an opposed orientation on their + // shared face if both of them follow the same orientation type + // (i.e., standard or non-standard). + // we verify this with a XOR operation. + const unsigned int neighbor_subface = + (!(this->line_orientation(face)) != + !(neighbor_cell->line_orientation(neighbor_face))) ? + (1 - subface) : + subface; + + const auto &info = ReferenceCell::internal::Info::get_cell( + ReferenceCell::Type::Tri); + const unsigned int neighbor_child_index = + info.child_cell_on_face(neighbor_face, neighbor_subface); + TriaIterator> sub_neighbor = + neighbor_cell->child(neighbor_child_index); + + // neighbor's child is not allowed to be further refined for the + // moment + Assert(sub_neighbor->refinement_case() == + RefinementCase::no_refinement, + ExcNotImplemented()); + + return sub_neighbor; + } + else if (this->reference_cell_type() == ReferenceCell::Type::Quad) + { + const unsigned int neighbor_neighbor = + this->neighbor_of_neighbor(face); + const unsigned int neighbor_child_index = + GeometryInfo::child_cell_on_face( + this->neighbor(face)->refinement_case(), + neighbor_neighbor, + subface); + + TriaIterator> sub_neighbor = + this->neighbor(face)->child(neighbor_child_index); + // the neighbors child can have children, + // which are not further refined along the + // face under consideration. as we are + // normally interested in one of this + // child's child, search for the right one. + while (sub_neighbor->has_children()) + { + Assert((GeometryInfo::face_refinement_case( + sub_neighbor->refinement_case(), + neighbor_neighbor) == + RefinementCase::no_refinement), + ExcInternalError()); + sub_neighbor = + sub_neighbor->child(GeometryInfo::child_cell_on_face( + sub_neighbor->refinement_case(), neighbor_neighbor, 0)); + } + + return sub_neighbor; } - return sub_neighbor; + // if no reference cell type matches + Assert(false, ExcNotImplemented()); + return TriaIterator>(); } case 3: { - // this function returns the neighbor's - // child on a given face and - // subface. - - // we have to consider one other aspect here: - // The face might be refined - // anisotropically. In this case, the subface - // number refers to the following, where we - // look at the face from the current cell, - // thus the subfaces are in standard - // orientation concerning the cell - // - // for isotropic refinement - // - // *---*---* - // | 2 | 3 | - // *---*---* - // | 0 | 1 | - // *---*---* - // - // for 2*anisotropic refinement - // (first cut_y, then cut_x) - // - // *---*---* - // | 2 | 3 | - // *---*---* - // | 0 | 1 | - // *---*---* - // - // for 2*anisotropic refinement - // (first cut_x, then cut_y) - // - // *---*---* - // | 1 | 3 | - // *---*---* - // | 0 | 2 | - // *---*---* - // - // for purely anisotropic refinement: - // - // *---*---* *-------* - // | | | | 1 | - // | 0 | 1 | or *-------* - // | | | | 0 | - // *---*---* *-------* - // - // for "mixed" refinement: - // - // *---*---* *---*---* *---*---* *-------* - // | | 2 | | 1 | | | 1 | 2 | | 2 | - // | 0 *---* or *---* 2 | or *---*---* or *---*---* - // | | 1 | | 0 | | | 0 | | 0 | 1 | - // *---*---* *---*---* *-------* *---*---* - - const typename Triangulation::face_iterator - mother_face = this->face(face); - const unsigned int total_children = mother_face->number_of_children(); - AssertIndexRange(subface, total_children); - Assert(total_children <= GeometryInfo<3>::max_children_per_face, - ExcInternalError()); - - unsigned int neighbor_neighbor; - TriaIterator> neighbor_child; - const TriaIterator> neighbor = - this->neighbor(face); - - - const RefinementCase mother_face_ref_case = - mother_face->refinement_case(); - if (mother_face_ref_case == - static_cast>( - RefinementCase<2>::cut_xy)) // total_children==4 + if (this->reference_cell_type() == ReferenceCell::Type::Hex) { - // this case is quite easy. we are sure, - // that the neighbor is not coarser. + // this function returns the neighbor's + // child on a given face and + // subface. + + // we have to consider one other aspect here: + // The face might be refined + // anisotropically. In this case, the subface + // number refers to the following, where we + // look at the face from the current cell, + // thus the subfaces are in standard + // orientation concerning the cell + // + // for isotropic refinement + // + // *---*---* + // | 2 | 3 | + // *---*---* + // | 0 | 1 | + // *---*---* + // + // for 2*anisotropic refinement + // (first cut_y, then cut_x) + // + // *---*---* + // | 2 | 3 | + // *---*---* + // | 0 | 1 | + // *---*---* + // + // for 2*anisotropic refinement + // (first cut_x, then cut_y) + // + // *---*---* + // | 1 | 3 | + // *---*---* + // | 0 | 2 | + // *---*---* + // + // for purely anisotropic refinement: + // + // *---*---* *-------* + // | | | | 1 | + // | 0 | 1 | or *-------* + // | | | | 0 | + // *---*---* *-------* + // + // for "mixed" refinement: + // + // *---*---* *---*---* *---*---* *-------* + // | | 2 | | 1 | | | 1 | 2 | | 2 | + // | 0 *---* or *---* 2 | or *---*---* or *---*---* + // | | 1 | | 0 | | | 0 | | 0 | 1 | + // *---*---* *---*---* *-------* *---*---* + + const typename Triangulation::face_iterator + mother_face = this->face(face); + const unsigned int total_children = + mother_face->number_of_children(); + AssertIndexRange(subface, total_children); + Assert(total_children <= GeometryInfo<3>::max_children_per_face, + ExcInternalError()); - // get the neighbor's number for the given - // face and the neighbor - neighbor_neighbor = this->neighbor_of_neighbor(face); + unsigned int neighbor_neighbor; + TriaIterator> neighbor_child; + const TriaIterator> neighbor = + this->neighbor(face); - // now use the info provided by GeometryInfo - // to extract the neighbors child number - const unsigned int neighbor_child_index = - GeometryInfo::child_cell_on_face( - neighbor->refinement_case(), - neighbor_neighbor, - subface, - neighbor->face_orientation(neighbor_neighbor), - neighbor->face_flip(neighbor_neighbor), - neighbor->face_rotation(neighbor_neighbor)); - neighbor_child = neighbor->child(neighbor_child_index); - - // make sure that the neighbor child cell we - // have found shares the desired subface. - Assert((this->face(face)->child(subface) == - neighbor_child->face(neighbor_neighbor)), - ExcInternalError()); - } - else //-> the face is refined anisotropically - { - // first of all, we have to find the - // neighbor at one of the anisotropic - // children of the - // mother_face. determine, which of - // these we need. - unsigned int first_child_to_find; - unsigned int neighbor_child_index; - if (total_children == 2) - first_child_to_find = subface; - else + + const RefinementCase mother_face_ref_case = + mother_face->refinement_case(); + if (mother_face_ref_case == + static_cast>( + RefinementCase<2>::cut_xy)) // total_children==4 { - first_child_to_find = subface / 2; - if (total_children == 3 && subface == 1 && - !mother_face->child(0)->has_children()) - first_child_to_find = 1; + // this case is quite easy. we are sure, + // that the neighbor is not coarser. + + // get the neighbor's number for the given + // face and the neighbor + neighbor_neighbor = this->neighbor_of_neighbor(face); + + // now use the info provided by GeometryInfo + // to extract the neighbors child number + const unsigned int neighbor_child_index = + GeometryInfo::child_cell_on_face( + neighbor->refinement_case(), + neighbor_neighbor, + subface, + neighbor->face_orientation(neighbor_neighbor), + neighbor->face_flip(neighbor_neighbor), + neighbor->face_rotation(neighbor_neighbor)); + neighbor_child = neighbor->child(neighbor_child_index); + + // make sure that the neighbor child cell we + // have found shares the desired subface. + Assert((this->face(face)->child(subface) == + neighbor_child->face(neighbor_neighbor)), + ExcInternalError()); } - if (neighbor_is_coarser(face)) + else //-> the face is refined anisotropically { - std::pair indices = - neighbor_of_coarser_neighbor(face); - neighbor_neighbor = indices.first; - - - // we have to translate our - // subface_index according to the - // RefineCase and subface index of - // the coarser face (our face is an - // anisotropic child of the coarser - // face), 'a' denotes our - // subface_index 0 and 'b' denotes - // our subface_index 1, whereas 0...3 - // denote isotropic subfaces of the - // coarser face - // - // cut_x and coarser_subface_index=0 - // - // *---*---* - // |b=2| | - // | | | - // |a=0| | - // *---*---* - // - // cut_x and coarser_subface_index=1 - // - // *---*---* - // | |b=3| - // | | | - // | |a=1| - // *---*---* - // - // cut_y and coarser_subface_index=0 - // - // *-------* - // | | - // *-------* - // |a=0 b=1| - // *-------* - // - // cut_y and coarser_subface_index=1 - // - // *-------* - // |a=2 b=3| - // *-------* - // | | - // *-------* - unsigned int iso_subface; - if (neighbor->face(neighbor_neighbor)->refinement_case() == - RefinementCase<2>::cut_x) - iso_subface = 2 * first_child_to_find + indices.second; + // first of all, we have to find the + // neighbor at one of the anisotropic + // children of the + // mother_face. determine, which of + // these we need. + unsigned int first_child_to_find; + unsigned int neighbor_child_index; + if (total_children == 2) + first_child_to_find = subface; else { - Assert( - neighbor->face(neighbor_neighbor)->refinement_case() == - RefinementCase<2>::cut_y, - ExcInternalError()); - iso_subface = first_child_to_find + 2 * indices.second; + first_child_to_find = subface / 2; + if (total_children == 3 && subface == 1 && + !mother_face->child(0)->has_children()) + first_child_to_find = 1; + } + if (neighbor_is_coarser(face)) + { + std::pair indices = + neighbor_of_coarser_neighbor(face); + neighbor_neighbor = indices.first; + + + // we have to translate our + // subface_index according to the + // RefineCase and subface index of + // the coarser face (our face is an + // anisotropic child of the coarser + // face), 'a' denotes our + // subface_index 0 and 'b' denotes + // our subface_index 1, whereas 0...3 + // denote isotropic subfaces of the + // coarser face + // + // cut_x and coarser_subface_index=0 + // + // *---*---* + // |b=2| | + // | | | + // |a=0| | + // *---*---* + // + // cut_x and coarser_subface_index=1 + // + // *---*---* + // | |b=3| + // | | | + // | |a=1| + // *---*---* + // + // cut_y and coarser_subface_index=0 + // + // *-------* + // | | + // *-------* + // |a=0 b=1| + // *-------* + // + // cut_y and coarser_subface_index=1 + // + // *-------* + // |a=2 b=3| + // *-------* + // | | + // *-------* + unsigned int iso_subface; + if (neighbor->face(neighbor_neighbor) + ->refinement_case() == RefinementCase<2>::cut_x) + iso_subface = 2 * first_child_to_find + indices.second; + else + { + Assert(neighbor->face(neighbor_neighbor) + ->refinement_case() == + RefinementCase<2>::cut_y, + ExcInternalError()); + iso_subface = + first_child_to_find + 2 * indices.second; + } + neighbor_child_index = + GeometryInfo::child_cell_on_face( + neighbor->refinement_case(), + neighbor_neighbor, + iso_subface, + neighbor->face_orientation(neighbor_neighbor), + neighbor->face_flip(neighbor_neighbor), + neighbor->face_rotation(neighbor_neighbor)); + } + else // neighbor is not coarser + { + neighbor_neighbor = neighbor_of_neighbor(face); + neighbor_child_index = + GeometryInfo::child_cell_on_face( + neighbor->refinement_case(), + neighbor_neighbor, + first_child_to_find, + neighbor->face_orientation(neighbor_neighbor), + neighbor->face_flip(neighbor_neighbor), + neighbor->face_rotation(neighbor_neighbor), + mother_face_ref_case); } - neighbor_child_index = GeometryInfo::child_cell_on_face( - neighbor->refinement_case(), - neighbor_neighbor, - iso_subface, - neighbor->face_orientation(neighbor_neighbor), - neighbor->face_flip(neighbor_neighbor), - neighbor->face_rotation(neighbor_neighbor)); - } - else // neighbor is not coarser - { - neighbor_neighbor = neighbor_of_neighbor(face); - neighbor_child_index = GeometryInfo::child_cell_on_face( - neighbor->refinement_case(), - neighbor_neighbor, - first_child_to_find, - neighbor->face_orientation(neighbor_neighbor), - neighbor->face_flip(neighbor_neighbor), - neighbor->face_rotation(neighbor_neighbor), - mother_face_ref_case); - } - - neighbor_child = neighbor->child(neighbor_child_index); - // it might be, that the neighbor_child - // has children, which are not refined - // along the given subface. go down that - // list and deliver the last of those. - while (neighbor_child->has_children() && - GeometryInfo::face_refinement_case( - neighbor_child->refinement_case(), neighbor_neighbor) == - RefinementCase<2>::no_refinement) - neighbor_child = - neighbor_child->child(GeometryInfo::child_cell_on_face( - neighbor_child->refinement_case(), neighbor_neighbor, 0)); - // if there are two total subfaces, we - // are finished. if there are four we - // have to get a child of our current - // neighbor_child. If there are three, - // we have to check which of the two - // possibilities applies. - if (total_children == 3) - { - if (mother_face->child(0)->has_children()) + neighbor_child = neighbor->child(neighbor_child_index); + // it might be, that the neighbor_child + // has children, which are not refined + // along the given subface. go down that + // list and deliver the last of those. + while ( + neighbor_child->has_children() && + GeometryInfo::face_refinement_case( + neighbor_child->refinement_case(), neighbor_neighbor) == + RefinementCase<2>::no_refinement) + neighbor_child = neighbor_child->child( + GeometryInfo::child_cell_on_face( + neighbor_child->refinement_case(), + neighbor_neighbor, + 0)); + + // if there are two total subfaces, we + // are finished. if there are four we + // have to get a child of our current + // neighbor_child. If there are three, + // we have to check which of the two + // possibilities applies. + if (total_children == 3) { - if (subface < 2) - neighbor_child = neighbor_child->child( - GeometryInfo::child_cell_on_face( - neighbor_child->refinement_case(), - neighbor_neighbor, - subface, - neighbor_child->face_orientation(neighbor_neighbor), - neighbor_child->face_flip(neighbor_neighbor), - neighbor_child->face_rotation(neighbor_neighbor), - mother_face->child(0)->refinement_case())); + if (mother_face->child(0)->has_children()) + { + if (subface < 2) + neighbor_child = neighbor_child->child( + GeometryInfo::child_cell_on_face( + neighbor_child->refinement_case(), + neighbor_neighbor, + subface, + neighbor_child->face_orientation( + neighbor_neighbor), + neighbor_child->face_flip(neighbor_neighbor), + neighbor_child->face_rotation( + neighbor_neighbor), + mother_face->child(0)->refinement_case())); + } + else + { + Assert(mother_face->child(1)->has_children(), + ExcInternalError()); + if (subface > 0) + neighbor_child = neighbor_child->child( + GeometryInfo::child_cell_on_face( + neighbor_child->refinement_case(), + neighbor_neighbor, + subface - 1, + neighbor_child->face_orientation( + neighbor_neighbor), + neighbor_child->face_flip(neighbor_neighbor), + neighbor_child->face_rotation( + neighbor_neighbor), + mother_face->child(1)->refinement_case())); + } } - else + else if (total_children == 4) { - Assert(mother_face->child(1)->has_children(), - ExcInternalError()); - if (subface > 0) - neighbor_child = neighbor_child->child( - GeometryInfo::child_cell_on_face( - neighbor_child->refinement_case(), - neighbor_neighbor, - subface - 1, - neighbor_child->face_orientation(neighbor_neighbor), - neighbor_child->face_flip(neighbor_neighbor), - neighbor_child->face_rotation(neighbor_neighbor), - mother_face->child(1)->refinement_case())); + neighbor_child = neighbor_child->child( + GeometryInfo::child_cell_on_face( + neighbor_child->refinement_case(), + neighbor_neighbor, + subface % 2, + neighbor_child->face_orientation(neighbor_neighbor), + neighbor_child->face_flip(neighbor_neighbor), + neighbor_child->face_rotation(neighbor_neighbor), + mother_face->child(subface / 2)->refinement_case())); } } - else if (total_children == 4) - { - neighbor_child = - neighbor_child->child(GeometryInfo::child_cell_on_face( - neighbor_child->refinement_case(), - neighbor_neighbor, - subface % 2, - neighbor_child->face_orientation(neighbor_neighbor), - neighbor_child->face_flip(neighbor_neighbor), - neighbor_child->face_rotation(neighbor_neighbor), - mother_face->child(subface / 2)->refinement_case())); - } - } - // it might be, that the neighbor_child has - // children, which are not refined along the - // given subface. go down that list and - // deliver the last of those. - while (neighbor_child->has_children()) - neighbor_child = - neighbor_child->child(GeometryInfo::child_cell_on_face( - neighbor_child->refinement_case(), neighbor_neighbor, 0)); + // it might be, that the neighbor_child has + // children, which are not refined along the + // given subface. go down that list and + // deliver the last of those. + while (neighbor_child->has_children()) + neighbor_child = + neighbor_child->child(GeometryInfo::child_cell_on_face( + neighbor_child->refinement_case(), neighbor_neighbor, 0)); #ifdef DEBUG - // check, whether the face neighbor_child matches the requested - // subface. - typename Triangulation::face_iterator requested; - switch (this->subface_case(face)) - { - case internal::SubfaceCase<3>::case_x: - case internal::SubfaceCase<3>::case_y: - case internal::SubfaceCase<3>::case_xy: - requested = mother_face->child(subface); - break; - case internal::SubfaceCase<3>::case_x1y2y: - case internal::SubfaceCase<3>::case_y1x2x: - requested = mother_face->child(subface / 2)->child(subface % 2); - break; - - case internal::SubfaceCase<3>::case_x1y: - case internal::SubfaceCase<3>::case_y1x: - switch (subface) - { - case 0: - case 1: - requested = mother_face->child(0)->child(subface); - break; - case 2: - requested = mother_face->child(1); - break; - default: - Assert(false, ExcInternalError()); - } - break; - case internal::SubfaceCase<3>::case_x2y: - case internal::SubfaceCase<3>::case_y2x: - switch (subface) - { - case 0: - requested = mother_face->child(0); - break; - case 1: - case 2: - requested = mother_face->child(1)->child(subface - 1); - break; - default: - Assert(false, ExcInternalError()); - } - break; - default: - Assert(false, ExcInternalError()); - break; - } - Assert(requested == neighbor_child->face(neighbor_neighbor), - ExcInternalError()); + // check, whether the face neighbor_child matches the requested + // subface. + typename Triangulation::face_iterator requested; + switch (this->subface_case(face)) + { + case internal::SubfaceCase<3>::case_x: + case internal::SubfaceCase<3>::case_y: + case internal::SubfaceCase<3>::case_xy: + requested = mother_face->child(subface); + break; + case internal::SubfaceCase<3>::case_x1y2y: + case internal::SubfaceCase<3>::case_y1x2x: + requested = + mother_face->child(subface / 2)->child(subface % 2); + break; + + case internal::SubfaceCase<3>::case_x1y: + case internal::SubfaceCase<3>::case_y1x: + switch (subface) + { + case 0: + case 1: + requested = mother_face->child(0)->child(subface); + break; + case 2: + requested = mother_face->child(1); + break; + default: + Assert(false, ExcInternalError()); + } + break; + case internal::SubfaceCase<3>::case_x2y: + case internal::SubfaceCase<3>::case_y2x: + switch (subface) + { + case 0: + requested = mother_face->child(0); + break; + case 1: + case 2: + requested = mother_face->child(1)->child(subface - 1); + break; + default: + Assert(false, ExcInternalError()); + } + break; + default: + Assert(false, ExcInternalError()); + break; + } + Assert(requested == neighbor_child->face(neighbor_neighbor), + ExcInternalError()); #endif - return neighbor_child; + return neighbor_child; + } + + // if no reference cell type matches + Assert(false, ExcNotImplemented()); + return TriaIterator>(); } default: - // 1d or more than 3d + // if 1d or more than 3d Assert(false, ExcNotImplemented()); return TriaIterator>(); } diff --git a/tests/simplex/refinement_03.cc b/tests/simplex/refinement_03.cc new file mode 100644 index 0000000000..5236adfd48 --- /dev/null +++ b/tests/simplex/refinement_03.cc @@ -0,0 +1,100 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + + + +// Verify CellAccessor::neighbor_child_on_subface for a triangle mesh. + + +#include + +#include + +#include "../tests.h" + + +template +void +test() +{ + // setup + // +---+---+ + // |\ | /| + // | \|/ | + // | +---+ + // | \ | + // | \| + // +-------+ + Triangulation tria; + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); + { + auto cell = tria.begin_active(); + (++cell)->set_refine_flag(); + Assert(++cell == tria.end(), ExcInternalError()); + } + tria.execute_coarsening_and_refinement(); + + // find face index on unrefined cell to neighboring cells + const auto & unrefined_cell = tria.begin_active(0); + unsigned int unrefined_f = numbers::invalid_unsigned_int; + for (unsigned int f = 0; f < unrefined_cell->n_faces(); ++f) + if (!unrefined_cell->face(f)->at_boundary()) + unrefined_f = f; + + // verify whether unrefined cell and neighboring children have matching + // vertices on their correpsonding subface + for (unsigned int sf = 0; sf < GeometryInfo::max_children_per_face; ++sf) + { + // unrefined vertex on subface + const Point &vertex_unrefined = + unrefined_cell->face(unrefined_f)->vertex(sf); + + // child on subface [! focus of this particular test !] + const auto &neighboring_child = + unrefined_cell->neighbor_child_on_subface(unrefined_f, sf); + + // face of child + unsigned int neighboring_child_f = numbers::invalid_unsigned_int; + for (unsigned int f = 0; f < neighboring_child->n_faces(); ++f) + if (!neighboring_child->face(f)->at_boundary() && + neighboring_child->neighbor(f) == unrefined_cell) + neighboring_child_f = f; + const auto &neighboring_child_face = + neighboring_child->face(neighboring_child_f); + + // find unrefined vertex among child face vertices + bool vertex_found = false; + for (unsigned int v = 0; v < neighboring_child_face->n_vertices(); ++v) + if (neighboring_child_face->vertex(v) == vertex_unrefined) + vertex_found = true; + + Assert(vertex_found, + ExcMessage( + "Wrong assignment of neighboring children to subfaces.")); + } + + deallog << "OK" << std::endl; +} + + +int +main() +{ + initlog(); + + deallog.push("2d"); + test<2>(); + deallog.pop(); +} diff --git a/tests/simplex/refinement_03.with_simplex_support=on.output b/tests/simplex/refinement_03.with_simplex_support=on.output new file mode 100644 index 0000000000..f445833684 --- /dev/null +++ b/tests/simplex/refinement_03.with_simplex_support=on.output @@ -0,0 +1,2 @@ + +DEAL:2d::OK -- 2.39.5