From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Wed, 13 Dec 2023 06:17:04 +0000 (-0700)
Subject: Disallow using direction flags for spacedim>dim+1.
X-Git-Tag: relicensing~246^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1bedb961b24b036487b00f701c85ade3f816a905;p=dealii.git

Disallow using direction flags for spacedim>dim+1.
---

diff --git a/doc/doxygen/headers/glossary.h b/doc/doxygen/headers/glossary.h
index 10d7a91672..62767ab8f3 100644
--- a/doc/doxygen/headers/glossary.h
+++ b/doc/doxygen/headers/glossary.h
@@ -827,7 +827,7 @@
  * can change all direction flags of a triangulation using the
  * Triangulation::flip_all_direction_flags() function.
  *
- * The flag is necessary to make cases like this work: assume we have a
+ * The flag is necessary to make cases like this work: Assume we have a
  * one-dimensional mesh embedded in a two-dimensional space,
  *
  *   @image html direction_flag.png "One dimensional mesh in two dimensions"
@@ -863,6 +863,14 @@
  * dimensions. We note that it would not be possible to find consistent
  * direction flags if the two-dimensional manifold is not orientable; such
  * manifolds are not currently supported by deal.II.
+ *
+ * Finally, the direction flag cannot be used for triangulations where
+ * `spacedim>dim+1`, such as for meshes with one-dimensional cells in 3d.
+ * In these cases, the normal vector to a cell does not simply point to
+ * one side or the other of a cell, but must lie in a two-dimensional sub-space
+ * perpendicular to the cell. As a consequence, we cannot make normal vectors
+ * consistent between adjacent cells simply by flipping it from one side to
+ * the other.
  * </dd>
  *
  *
diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h
index 40b326e65e..48066c61b7 100644
--- a/include/deal.II/grid/tria.h
+++ b/include/deal.II/grid/tria.h
@@ -2098,10 +2098,12 @@ public:
       &construction_data);
 
   /**
-   * Revert or flip the direction_flags of a dim<spacedim triangulation, see
+   * Revert or flip the direction flags of a triangulation with
+   * `dim==spacedim-1`, see
    * @ref GlossDirectionFlag.
    *
-   * This function throws an exception if dim equals spacedim.
+   * This function throws an exception if `dim==spacedim` or if
+   * `dim<spacedim-1`.
    */
   void
   flip_all_direction_flags();
diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h
index 651415465d..e5bc974b3b 100644
--- a/include/deal.II/grid/tria_accessor.h
+++ b/include/deal.II/grid/tria_accessor.h
@@ -3858,7 +3858,10 @@ public:
    */
 
   /**
-   * Return the orientation of this cell.
+   * Return the orientation of this cell. This function always returns
+   * `true` if `dim==spacedim`. It can return `true` or `false` if
+   * `dim==spacedim-1`. The function cannot be called (and will abort
+   * with an error) if called for `dim<spacedim-1`.
    *
    * For the meaning of this flag, see
    * @ref GlossDirectionFlag.
@@ -4172,7 +4175,11 @@ private:
   set_parent(const unsigned int parent_index);
 
   /**
-   * Set the orientation of this cell.
+   * Set the orientation of this cell. This function can only be
+   * called if the argument is `true` if `dim==spacedim`. It can be
+   * called with either `true` or `false` if `dim==spacedim-1`. The
+   * function cannot be called (and will abort with an error) if called
+   * for `dim<spacedim-1`.
    *
    * For the meaning of this flag, see
    * @ref GlossDirectionFlag.
diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index 085fa32595..28011c22f7 100644
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -4947,7 +4947,7 @@ namespace internal
         // already cleared at the
         // beginning of this function
 
-        if (dim < spacedim)
+        if (dim == spacedim - 1)
           for (unsigned int c = 0; c < n_children; ++c)
             cell->child(c)->set_direction_flag(cell->direction_flag());
       }
@@ -5414,7 +5414,7 @@ namespace internal
 
           cell->set_refinement_case(ref_case);
 
-          if (dim < spacedim)
+          if (dim == spacedim - 1)
             for (unsigned int c = 0; c < n_children; ++c)
               cell->child(c)->set_direction_flag(cell->direction_flag());
         };
@@ -5595,7 +5595,8 @@ namespace internal
                   first_child->set_material_id(cell->material_id());
                   first_child->set_manifold_id(cell->manifold_id());
                   first_child->set_subdomain_id(subdomainid);
-                  first_child->set_direction_flag(cell->direction_flag());
+                  if (dim == spacedim - 1)
+                    first_child->set_direction_flag(cell->direction_flag());
 
                   first_child->set_parent(cell->index());
 
@@ -5646,7 +5647,8 @@ namespace internal
                   second_child->set_material_id(cell->material_id());
                   second_child->set_manifold_id(cell->manifold_id());
                   second_child->set_subdomain_id(subdomainid);
-                  second_child->set_direction_flag(cell->direction_flag());
+                  if (dim == spacedim - 1)
+                    second_child->set_direction_flag(cell->direction_flag());
 
                   if (cell->neighbor(1).state() != IteratorState::valid)
                     second_child->set_neighbor(1, cell->neighbor(1));
@@ -12929,7 +12931,8 @@ DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
 void Triangulation<dim, spacedim>::flip_all_direction_flags()
 {
   AssertThrow(dim + 1 == spacedim,
-              ExcMessage("Only works for dim == spacedim-1"));
+              ExcMessage(
+                "This function can only be called if dim == spacedim-1."));
   for (const auto &cell : this->active_cell_iterators())
     cell->set_direction_flag(!cell->direction_flag());
 }
diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc
index 0836ac904c..4c084ea65f 100644
--- a/source/grid/tria_accessor.cc
+++ b/source/grid/tria_accessor.cc
@@ -2225,11 +2225,17 @@ bool
 CellAccessor<dim, spacedim>::direction_flag() const
 {
   Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
-  if (dim == spacedim)
+  if constexpr (dim == spacedim)
     return true;
-  else
+  else if constexpr (dim == spacedim - 1)
     return this->tria->levels[this->present_level]
       ->direction_flags[this->present_index];
+  else
+    {
+      Assert(false,
+             ExcMessage("This function cannot be called if dim<spacedim-1."));
+      return true;
+    }
 }
 
 
@@ -2240,13 +2246,17 @@ CellAccessor<dim, spacedim>::set_direction_flag(
   const bool new_direction_flag) const
 {
   Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
-  if (dim < spacedim)
+  if constexpr (dim == spacedim)
+    Assert(new_direction_flag == true,
+           ExcMessage("If dim==spacedim, direction flags are always true and "
+                      "can not be set to anything else."));
+  else if constexpr (dim == spacedim - 1)
     this->tria->levels[this->present_level]
       ->direction_flags[this->present_index] = new_direction_flag;
   else
     Assert(new_direction_flag == true,
-           ExcMessage("If dim==spacedim, direction flags are always true and "
-                      "can not be set to anything else."));
+           ExcMessage("If dim<spacedim-1, then this function can be called "
+                      "only if the argument is 'true'."));
 }