Regular cylinder uses subdivided_cylinder

diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index 3f2fb6382e6212adfde21cb31ce1c8454dc2b2f1..160681e6ef79dcc4b28bb43299b3c5cd33db2f4b 100644 (file)
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -5461,122 +5461,7 @@ namespace GridGenerator
                 const double      radius,
                 const double      half_length)
   {
-    // Copy the base from hyper_ball<3>
-    // and transform it to yz
-    const double d            = radius / std::sqrt(2.0);
-    const double a            = d / (1 + std::sqrt(2.0));
-    Point<3>     vertices[24] = {
-      Point<3>(-d, -half_length, -d),
-      Point<3>(d, -half_length, -d),
-      Point<3>(-a, -half_length, -a),
-      Point<3>(a, -half_length, -a),
-      Point<3>(-a, -half_length, a),
-      Point<3>(a, -half_length, a),
-      Point<3>(-d, -half_length, d),
-      Point<3>(d, -half_length, d),
-      Point<3>(-d, 0, -d),
-      Point<3>(d, 0, -d),
-      Point<3>(-a, 0, -a),
-      Point<3>(a, 0, -a),
-      Point<3>(-a, 0, a),
-      Point<3>(a, 0, a),
-      Point<3>(-d, 0, d),
-      Point<3>(d, 0, d),
-      Point<3>(-d, half_length, -d),
-      Point<3>(d, half_length, -d),
-      Point<3>(-a, half_length, -a),
-      Point<3>(a, half_length, -a),
-      Point<3>(-a, half_length, a),
-      Point<3>(a, half_length, a),
-      Point<3>(-d, half_length, d),
-      Point<3>(d, half_length, d),
-    };
-    // Turn cylinder such that y->x
-    for (auto &vertex : vertices)
-      {
-        const double h = vertex(1);
-        vertex(1)      = -vertex(0);
-        vertex(0)      = h;
-      }
-
-    int cell_vertices[10][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
-                                {0, 2, 8, 10, 6, 4, 14, 12},
-                                {2, 3, 10, 11, 4, 5, 12, 13},
-                                {1, 7, 9, 15, 3, 5, 11, 13},
-                                {6, 4, 14, 12, 7, 5, 15, 13}};
-    for (unsigned int i = 0; i < 5; ++i)
-      for (unsigned int j = 0; j < 8; ++j)
-        cell_vertices[i + 5][j] = cell_vertices[i][j] + 8;
-
-    std::vector<CellData<3>> cells(10, CellData<3>());
-
-    for (unsigned int i = 0; i < 10; ++i)
-      {
-        for (unsigned int j = 0; j < 8; ++j)
-          cells[i].vertices[j] = cell_vertices[i][j];
-        cells[i].material_id = 0;
-      }
-
-    tria.create_triangulation(std::vector<Point<3>>(std::begin(vertices),
-                                                    std::end(vertices)),
-                              cells,
-                              SubCellData()); // no boundary information
-
-    // set boundary indicators for the
-    // faces at the ends to 1 and 2,
-    // respectively. note that we also
-    // have to deal with those lines
-    // that are purely in the interior
-    // of the ends. we determine whether
-    // an edge is purely in the
-    // interior if one of its vertices
-    // is at coordinates '+-a' as set
-    // above
-    Triangulation<3>::cell_iterator cell = tria.begin();
-    Triangulation<3>::cell_iterator end  = tria.end();
-
-    tria.set_all_manifold_ids_on_boundary(0);
-
-    for (; cell != end; ++cell)
-      for (unsigned int i : GeometryInfo<3>::face_indices())
-        if (cell->at_boundary(i))
-          {
-            if (cell->face(i)->center()(0) > half_length - 1.e-5)
-              {
-                cell->face(i)->set_boundary_id(2);
-                cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
-
-                for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
-                     ++e)
-                  if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
-                      (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
-                      (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
-                      (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
-                    {
-                      cell->face(i)->line(e)->set_boundary_id(2);
-                      cell->face(i)->line(e)->set_manifold_id(
-                        numbers::flat_manifold_id);
-                    }
-              }
-            else if (cell->face(i)->center()(0) < -half_length + 1.e-5)
-              {
-                cell->face(i)->set_boundary_id(1);
-                cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
-
-                for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
-                     ++e)
-                  if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
-                      (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
-                      (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
-                      (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
-                    {
-                      cell->face(i)->line(e)->set_boundary_id(1);
-                      cell->face(i)->line(e)->set_manifold_id(
-                        numbers::flat_manifold_id);
-                    }
-              }
-          }
-    tria.set_manifold(0, CylindricalManifold<3>());
+    subdivided_cylinder(tria, 2, radius, half_length);
   }
 
   // Implementation for 3D only