Adapt CellAccessor::neighbor_child_on_subface for triangles.

Required by DoFTools::make_hanging_node_constraints.

diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index 7bf0fadc60d6ed17311f114ba1bac834e9713a84..b830513b07bba7cea533e5de5e3bbc25c4fc2595 100644 (file)
--- 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 bf6f98c58fa054eb4c94d5ffc553a811ad6e492f..0611d25e3d2f08e8aa817101991b279b6fabac3b 100644 (file)
--- a/source/grid/tria_accessor.cc
+++ b/source/grid/tria_accessor.cc
@@ -2845,368 +2845,433 @@ CellAccessor<dim, spacedim>::neighbor_child_on_subface(
     {
       case 2:
         {
-          const unsigned int neighbor_neighbor =
-            this->neighbor_of_neighbor(face);
-          const unsigned int neighbor_child_index =
-            GeometryInfo<dim>::child_cell_on_face(
-              this->neighbor(face)->refinement_case(),
-              neighbor_neighbor,
-              subface);
-
-          TriaIterator<CellAccessor<dim, spacedim>> 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<dim>::face_refinement_case(
-                        sub_neighbor->refinement_case(), neighbor_neighbor) ==
-                      RefinementCase<dim>::no_refinement),
-                     ExcInternalError());
-              sub_neighbor =
-                sub_neighbor->child(GeometryInfo<dim>::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<CellAccessor<dim, spacedim>> 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<dim>::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<dim>::child_cell_on_face(
+                  this->neighbor(face)->refinement_case(),
+                  neighbor_neighbor,
+                  subface);
+
+              TriaIterator<CellAccessor<dim, spacedim>> 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<dim>::face_refinement_case(
+                            sub_neighbor->refinement_case(),
+                            neighbor_neighbor) ==
+                          RefinementCase<dim>::no_refinement),
+                         ExcInternalError());
+                  sub_neighbor =
+                    sub_neighbor->child(GeometryInfo<dim>::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<CellAccessor<dim, spacedim>>();
         }
 
 
       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<dim, spacedim>::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<CellAccessor<dim, spacedim>>       neighbor_child;
-          const TriaIterator<CellAccessor<dim, spacedim>> neighbor =
-            this->neighbor(face);
-
-
-          const RefinementCase<dim - 1> mother_face_ref_case =
-            mother_face->refinement_case();
-          if (mother_face_ref_case ==
-              static_cast<RefinementCase<dim - 1>>(
-                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<dim, spacedim>::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<CellAccessor<dim, spacedim>>       neighbor_child;
+              const TriaIterator<CellAccessor<dim, spacedim>> neighbor =
+                this->neighbor(face);
 
-              // now use the info provided by GeometryInfo
-              // to extract the neighbors child number
-              const unsigned int neighbor_child_index =
-                GeometryInfo<dim>::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<dim - 1> mother_face_ref_case =
+                mother_face->refinement_case();
+              if (mother_face_ref_case ==
+                  static_cast<RefinementCase<dim - 1>>(
+                    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<dim>::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<unsigned int, unsigned int> 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<unsigned int, unsigned int> 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<dim>::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<dim>::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<dim>::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<dim>::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<dim>::face_refinement_case(
-                       neighbor_child->refinement_case(), neighbor_neighbor) ==
-                       RefinementCase<2>::no_refinement)
-                neighbor_child =
-                  neighbor_child->child(GeometryInfo<dim>::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<dim>::face_refinement_case(
+                      neighbor_child->refinement_case(), neighbor_neighbor) ==
+                      RefinementCase<2>::no_refinement)
+                    neighbor_child = neighbor_child->child(
+                      GeometryInfo<dim>::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<dim>::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<dim>::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<dim>::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<dim>::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<dim>::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<dim>::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<dim>::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<dim>::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<dim, spacedim>::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<dim, spacedim>::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<CellAccessor<dim, spacedim>>();
         }
 
       default:
-        // 1d or more than 3d
+        // if 1d or more than 3d
         Assert(false, ExcNotImplemented());
         return TriaIterator<CellAccessor<dim, spacedim>>();
     }

diff --git a/tests/simplex/refinement_03.cc b/tests/simplex/refinement_03.cc
new file mode 100644 (file)
index 0000000..5236adf
--- /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 <deal.II/grid/tria.h>
+
+#include <deal.II/simplex/grid_generator.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test()
+{
+  // setup
+  //   +---+---+
+  //   |\  |  /|
+  //   |  \|/  |
+  //   |   +---+
+  //   |    \  |
+  //   |      \|
+  //   +-------+
+  Triangulation<dim> 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<dim>::max_children_per_face; ++sf)
+    {
+      // unrefined vertex on subface
+      const Point<dim> &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 (file)
index 0000000..f445833
--- /dev/null
+++ b/tests/simplex/refinement_03.with_simplex_support=on.output
@@ -0,0 +1,2 @@
+
+DEAL:2d::OK