Fix vertex orientations in torus mesh

doc/news/changes/minor/231122KronbichlerAnselmannMunch

@@ -0,0 +1,5 @@
+Fixed: Meshes created with GridGenerator::torus() would previously run
+into an assertion in p4est when using parallel::distributed::Triangulation. This
+is now fixed.
+<br>
+(Martin Kronbichler, Mathias Anselmann, Peter Munch, 2023/11/22)

source/grid/grid_generator.cc

@@ -2002,9 +2002,9 @@ namespace GridGenerator
            ExcMessage("Outer radius R must be greater than the inner "
                       "radius r."));
     Assert(r > 0.0, ExcMessage("The inner radius r must be positive."));
-    Assert(n_cells_toroidal > 2,
+    Assert(n_cells_toroidal > static_cast<unsigned int>(phi / numbers::PI),
            ExcMessage("Number of cells in toroidal direction has "
-                      "to be at least 3."));
+                      "to be at least 3 for a torus of polar extent 2*pi."));
     AssertThrow(phi > 0.0 && phi < 2.0 * numbers::PI + 1.0e-15,
                 ExcMessage("Invalid angle phi specified."));
 
@@ -2063,20 +2063,20 @@ namespace GridGenerator
             cells[5 * c + 1].vertices[2 + j * 4] = offset + 4;
             cells[5 * c + 1].vertices[3 + j * 4] = offset + 5;
             // cell 2 in x-y-plane
-            cells[5 * c + 2].vertices[0 + j * 4] = offset + 4;
-            cells[5 * c + 2].vertices[1 + j * 4] = offset + 5;
-            cells[5 * c + 2].vertices[2 + j * 4] = offset + 6;
-            cells[5 * c + 2].vertices[3 + j * 4] = offset + 7;
+            cells[5 * c + 2].vertices[0 + j * 4] = offset + 6;
+            cells[5 * c + 2].vertices[1 + j * 4] = offset + 4;
+            cells[5 * c + 2].vertices[2 + j * 4] = offset + 7;
+            cells[5 * c + 2].vertices[3 + j * 4] = offset + 5;
             // cell 3 in x-y-plane
             cells[5 * c + 3].vertices[0 + j * 4] = offset + 0;
             cells[5 * c + 3].vertices[1 + j * 4] = offset + 2;
             cells[5 * c + 3].vertices[2 + j * 4] = offset + 6;
             cells[5 * c + 3].vertices[3 + j * 4] = offset + 4;
             // cell 4 in x-y-plane
-            cells[5 * c + 4].vertices[0 + j * 4] = offset + 3;
-            cells[5 * c + 4].vertices[1 + j * 4] = offset + 1;
-            cells[5 * c + 4].vertices[2 + j * 4] = offset + 5;
-            cells[5 * c + 4].vertices[3 + j * 4] = offset + 7;
+            cells[5 * c + 4].vertices[0 + j * 4] = offset + 1;
+            cells[5 * c + 4].vertices[1 + j * 4] = offset + 7;
+            cells[5 * c + 4].vertices[2 + j * 4] = offset + 3;
+            cells[5 * c + 4].vertices[3 + j * 4] = offset + 5;
           }
 
         cells[5 * c].material_id = 0;

tests/distributed_grids/create_torus_mesh.cc

@@ -0,0 +1,59 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 - 2023 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 GridGenerator::torus, which used to run into an assertion with p4est
+// due to bad orientations before 2023, see
+// https://github.com/dealii/dealii/issues/16272
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include "../tests.h"
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  initlog();
+
+  // create mesh
+  const int dim = 3;
+
+  parallel::distributed::Triangulation<dim> triangulation(
+    MPI_COMM_WORLD,
+    typename Triangulation<dim>::MeshSmoothing(
+      Triangulation<dim>::smoothing_on_refinement |
+      Triangulation<dim>::smoothing_on_coarsening));
+
+  const double       R                = 2.;
+  const double       r                = 1.;
+  const unsigned int n_cells_toroidal = 1;
+  const double       phi              = 0.5 * dealii::numbers::PI;
+  dealii::GridGenerator::torus<dim, dim>(
+    triangulation, R, r, n_cells_toroidal, phi);
+  triangulation.reset_all_manifolds();
+  triangulation.refine_global(1);
+
+  deallog << "Tria info: " << triangulation.n_vertices() << " "
+          << triangulation.n_cells() << " " << triangulation.n_faces() << " "
+          << triangulation.n_raw_lines() << std::endl;
+
+  return 0;
+}

tests/distributed_grids/create_torus_mesh.output

@@ -0,0 +1,2 @@
+
+DEAL::Tria info: 75 45 170 214