]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Have all GridGenerator functions attach curved manifolds.
authorLuca Heltai <luca.heltai@sissa.it>
Tue, 17 Apr 2018 08:22:05 +0000 (10:22 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Tue, 17 Apr 2018 08:28:25 +0000 (10:28 +0200)
Co-authored-by: David Wells <drwells@vt.edu>
include/deal.II/grid/grid_generator.h
source/grid/grid_generator.cc

index 064ff48464604d858f6b9cdb2dd892465f2c8f69..ac3451610399ef3ae3bdb6644e34e9a9f9d4aabd 100644 (file)
@@ -33,11 +33,18 @@ DEAL_II_NAMESPACE_OPEN
  * This namespace provides a collection of functions for generating
  * triangulations for some basic geometries.
  *
- * Some of these functions receive a flag @p colorize. If this is set, parts
- * of the boundary receive different
- * @ref GlossBoundaryIndicator "boundary indicators"),
- * allowing them to be distinguished for the purpose of attaching geometry
- * objects and evaluating different boundary conditions.
+ * Some of these functions receive a flag @p colorize. If this is set, parts of
+ * the boundary receive different @ref GlossBoundaryIndicator "boundary
+ * indicators" allowing them to be distinguished for the purpose of evaluating
+ * different boundary conditions.
+ *
+ * If the domain is curved, each of the domain parts that should be
+ * refined by following an appropriate Manifold description will
+ * receive a different @ref GlossManifoldIndicator "manifold
+ * indicator", and the correct Manifold descriptor will be attached to
+ * the Triangulation. Notice that if you later tranform the
+ * triangulation, you have to make sure you attach the correct new Manifold
+ * to the triangulation.
  *
  * @ingroup grid
  */
@@ -54,7 +61,7 @@ namespace GridGenerator
    * tensor product interval $[left,right]^{\text{dim}}$ in the present number
    * of dimensions, where the limits are given as arguments. They default to
    * zero and unity, then producing the unit hypercube.
-   *
+               *
    * If the argument @p colorize is false, all boundary indicators are set to
    * zero ("not colorized") for 2d and 3d. If it is true, the boundary is
    * colorized as in hyper_rectangle(). In 1d the indicators are always
@@ -138,18 +145,14 @@ namespace GridGenerator
    * GeometryInfo class. Importantly, however, in 3d colorization does not set
    * @p boundary_ids of <i>edges</i>, but only of <i>faces</i>, because each
    * boundary edge is shared between two faces and it is not clear how the
-   * boundary id of an edge should be set in that case. This may later on lead
-   * to problems if one wants to assign boundary or manifold objects to parts
-   * of the boundary with certain boundary indicators since then the boundary
-   * object may not apply to the edges bounding the face it is meant to
-   * describe.
+   * boundary id of an edge should be set in that case.
    *
    * Additionally, if @p colorize is @p true, material ids are assigned to the
    * cells according to the octant their center is in: being in the right half
    * space for any coordinate direction <i>x<sub>i</sub></i> adds
    * 2<sup>i</sup>. For instance, a cell with center point (1,-1,1) yields a
    * material id 5, assuming that the center of the hyper rectangle lies at
-   * the origin.
+   * the origin. No manifold id is set for the cells.
    *
    * If @p dim < @p spacedim, this will create a @p dim dimensional object in
    * the first @p dim coordinate directions embedded into the @p spacedim
@@ -176,9 +179,9 @@ namespace GridGenerator
    * direction is 1.
    *
    * If the @p colorize flag is set, the @p boundary_ids of the surfaces are
-   * assigned, such that the lower one in @p x-direction is 0, the upper one
-   * is 1 (the left and the right vertical face). The indicators for the
-   * surfaces in @p y-direction are 2 and 3, the ones for @p z are 4 and 5.
+   * assigned, such that the lower one in @p x-direction is 0, the upper one is
+   * 1 (the left and the right vertical face). The indicators for the surfaces
+   * in @p y-direction are 2 and 3, the ones for @p z are 4 and 5.
    * Additionally, material ids are assigned to the cells according to the
    * octant their center is in: being in the right half plane for any
    * coordinate direction <i>x<sub>i</sub></i> adds 2<sup>i</sup>. For
@@ -260,7 +263,7 @@ namespace GridGenerator
                               const std::vector< std::vector<double> > &spacing,
                               const Point<dim>                         &p,
                               const Table<dim,types::material_id>      &material_id,
-                              const bool                                colorize=false);
+                              const bool                               colorize=false);
 
   /**
    * \brief Rectangular domain with rectangular pattern of holes
@@ -452,29 +455,36 @@ namespace GridGenerator
    * This function is declared to exist for triangulations of all space
    * dimensions, but throws an error if called in 1d.
    *
-   * You should attach a SphericalManifold to the cells and faces for correct
-   * placement of vertices upon refinement and to be able to use higher order
-   * mappings. However, it turns out that creating a mesh for a hyperball is
-   * not entirely trivial since the central cell has to be treated
-   * differently than the "cap" cells. The "Possibilities for extensions"
-   * section of step-6 has an extensive discussion of how one would construct
-   * such meshes and what one needs to do for it.
+   * By default, the manifold_id is set to 0 on the boundary faces, 1 on the
+   * the boundary cells, and types::flat_manifold_id on the central cell and on
+   * internal faces.
+   *
+   * A SphericalManifold is attached by default to the boundary faces for
+   * correct placement of boundary vertices upon refinement and to be able to
+   * use higher order mappings. However, it turns out that this strategy may
+   * not be the optimal one to create a good a mesh for a hyperball. The
+   * "Possibilities for extensions" section of step-6 has an extensive
+   * discussion of how one would construct better meshes and what one needs to
+   * do for it. Selecting the argument @p
+   * attach_spherical_manifold_on_boundary_cells to true attaches a
+   * SphericalManifold manifold also to the boundary cells, and not only to the
+   * boundary faces.
    *
    * @note The triangulation passed as argument needs to be empty when calling this function.
    */
   template <int dim>
   void hyper_ball (Triangulation<dim> &tria,
                    const Point<dim>   &center = Point<dim>(),
-                   const double        radius = 1.);
+                   const double        radius = 1.,
+                   const bool attach_spherical_manifold_on_boundary_cells=false);
 
   /**
    * Creates a hyper sphere, i.e., a surface of a ball in @p spacedim
    * dimensions. This function only exists for dim+1=spacedim in 2 and 3 space
    * dimensions. (To create a mesh of a ball, use GridGenerator::hyper_ball().)
    *
-   * You should attach a SphericalManifold to the cells and faces for correct
-   * placement of vertices upon refinement and to be able to use higher order
-   * mappings.
+   * By default, all manifold ids of the triangulation are set to zero, and a
+   * SphericalManifold is attached to the grid.
    *
    * The following pictures are generated with:
    * @code
@@ -483,9 +493,6 @@ namespace GridGenerator
    * static SphericalManifold<2,3> surface_description;
    *
    * GridGenerator::hyper_sphere(triangulation);
-   *
-   * triangulation.set_all_manifold_ids(0);
-   * triangulation.set_manifold (0, surface_description);
    * triangulation.refine_global(3);
    * @endcode
    *
@@ -509,11 +516,11 @@ namespace GridGenerator
    * relative to @p center, which contains three elements in 2d and four in 3d.
    *
    * The boundary indicators for the final triangulation are 0 for the curved
-   * boundary and 1 for the cut plane.
-   *
-   * The appropriate boundary class is HyperBallBoundary.
+   * boundary and 1 for the cut plane. The manifold id for the curved boundary is
+   * set to zero, and a SphericalManifold is attached to it.
    *
-   * @note The triangulation passed as argument needs to be empty when calling this function.
+   * @note The triangulation passed as argument needs to be empty when calling
+   * this function.
    */
   template <int dim>
   void quarter_hyper_ball (Triangulation<dim> &tria,
@@ -526,10 +533,8 @@ namespace GridGenerator
    * <i>x</i>-axis.
    *
    * The boundary indicators for the final triangulation are 0 for the curved
-   * boundary and 1 for the cut plane.
-   *
-   * The appropriate boundary class is HalfHyperBallBoundary, or
-   * HyperBallBoundary.
+   * boundary and 1 for the cut plane. The manifold id for the curved boundary is
+   * set to zero, and a SphericalManifold is attached to it.
    *
    * @note The triangulation passed as argument needs to be empty when calling this function.
    */
@@ -556,6 +561,9 @@ namespace GridGenerator
    * the GridTools::transform() function using a rotation operator as
    * argument.
    *
+   * The manifold id for the hull of the cylinder is set to zero, and a
+   * CylindricalManifold is attached to it.
+   *
    * @note The triangulation passed as argument needs to be empty when calling this function.
    */
   template <int dim>
@@ -578,11 +586,10 @@ namespace GridGenerator
    *
    * The boundaries are colored according to the following scheme: 0 for the
    * hull of the cone, 1 for the left hand face and 2 for the right hand face.
+   * Both the boundary indicators and the manifold indicators are set.
    *
-   * In three dimensions, the CylindricalManifold class is an appropriate choice
-   * for the description of the hull, with which you probably want to associate
-   * boundary indicator 0.
-   * In two dimensions the default FlatManifold is sufficient.
+   * In three dimensions, the manifold id of the hull is set to zero, and a
+   * CylindricalManifold is attached to it.
    *
    * @note The triangulation passed as argument needs to be empty when calling this function.
    *
@@ -702,11 +709,8 @@ namespace GridGenerator
    * both the faces and the edges of these boundaries. If the flag is @p
    * false, both have indicator zero.
    *
-   * You should attach a SphericalManifold to the cells and faces for correct
-   * placement of vertices upon refinement and to be able to use higher order
-   * mappings. Alternatively, it is also possible to attach a
-   * HyperShellBoundary to the inner and outer boundary. This will create
-   * inferior meshes as described below.
+   * All manifold ids are set to zero, and a SphericalManifold is attached to the
+   * the triangulation.
    *
    * In 2d, the number <tt>n_cells</tt> of elements for this initial
    * triangulation can be chosen arbitrarily. If the number of initial cells
@@ -781,6 +785,9 @@ namespace GridGenerator
    * boundary where $x=0$, get indicator 0, 1, and 2, respectively. Otherwise
    * all indicators are set to 0.
    *
+   * All manifold ids are set to zero, and a  SphericalManifold is attached
+   * to the triangulation.
+   *
    * @note The triangulation passed as argument needs to be empty when calling this function.
    */
   template <int dim>
@@ -876,9 +883,16 @@ namespace GridGenerator
 
 
   /**
-   * This class produces a square in the <i>xy</i>-plane with a circular hole
-   * in the middle. Square and circle are centered at the origin. In 3d, this
-   * geometry is extruded in $z$ direction to the interval $[0,L]$.
+   * This class produces a square in the <i>xy</i>-plane with a cylindrical
+   * hole in the middle. The square and the circle are centered at the
+   * origin. In 3d, this geometry is extruded in $z$ direction to the interval
+   * $[0,L]$.
+   *
+   * The inner boundary has a manifold id of $0$ and a boundary id of
+   * $6$. This function attaches a PolarManifold or CylindricalManifold to the
+   * interior boundary in 2d and 3d respectively. The other faces have
+   * boundary ids of $0, 1, 2, 3, 4$, or $5$ given in the standard order of
+   * faces in 2d or 3d.
    *
    * @image html cubes_hole.png
    *
index 9f461d62b3eeae372c6635221224368eaf7e1a84..b4e09ed7320f44db96ba115c1fb8956ae8732297 100644 (file)
@@ -1985,7 +1985,8 @@ namespace GridGenerator
   template <>
   void hyper_ball (Triangulation<1> &,
                    const Point<1> &,
-                   const double)
+                   const double,
+                   const bool)
   {
     Assert (false, ExcNotImplemented());
   }
@@ -2265,7 +2266,6 @@ namespace GridGenerator
         cell->face(0)->set_boundary_id(0);
         cell->face(1)->set_boundary_id(4);
         cell->face(3)->set_boundary_id(5);
-
       }
 
   }
@@ -2277,7 +2277,8 @@ namespace GridGenerator
   void
   hyper_ball (Triangulation<2> &tria,
               const Point<2>   &p,
-              const double      radius)
+              const double      radius,
+              const bool internal_manifolds)
   {
     // equilibrate cell sizes at
     // transition from the inner part
@@ -2307,12 +2308,17 @@ namespace GridGenerator
         for (unsigned int j=0; j<4; ++j)
           cells[i].vertices[j] = cell_vertices[i][j];
         cells[i].material_id = 0;
+        cells[i].manifold_id = i==2 ? 1 : numbers::flat_manifold_id;
       };
 
     tria.create_triangulation (
       std::vector<Point<2> >(std::begin(vertices), std::end(vertices)),
       cells,
       SubCellData());       // no boundary information
+    tria.set_all_manifold_ids_on_boundary(0);
+    tria.set_manifold(0, SphericalManifold<2>(p));
+    if (internal_manifolds)
+      tria.set_manifold(1, SphericalManifold<2>(p));
   }
 
 
@@ -2379,6 +2385,9 @@ namespace GridGenerator
 
     if (colorize)
       colorize_hyper_shell(tria, center, inner_radius, outer_radius);
+
+    tria.set_all_manifold_ids(0);
+    tria.set_manifold(0, SphericalManifold<2>(center));
   }
 
 
@@ -2473,6 +2482,8 @@ namespace GridGenerator
     Triangulation<dim>::cell_iterator cell = tria.begin();
     Triangulation<dim>::cell_iterator end = tria.end();
 
+    tria.set_all_manifold_ids_on_boundary(0);
+
     while (cell != end)
       {
         for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
@@ -2484,10 +2495,14 @@ namespace GridGenerator
             // component of the center, then this is part of the plane
             if (cell->face(i)->center()(0) < p(0)+1.e-5 * radius
                 || cell->face(i)->center()(1) < p(1)+1.e-5 * radius)
-              cell->face(i)->set_boundary_id(1);
+              {
+                cell->face(i)->set_boundary_id(1);
+                cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
+              }
           }
         ++cell;
       }
+    tria.set_manifold(0, SphericalManifold<2>(p));
   }
 
 
@@ -2534,6 +2549,7 @@ namespace GridGenerator
     Triangulation<2>::cell_iterator cell = tria.begin();
     Triangulation<2>::cell_iterator end = tria.end();
 
+    tria.set_all_manifold_ids_on_boundary(0);
 
     while (cell != end)
       {
@@ -2544,10 +2560,14 @@ namespace GridGenerator
 
             // If x is zero, then this is part of the plane
             if (cell->face(i)->center()(0) < p(0)+1.e-5 * radius)
-              cell->face(i)->set_boundary_id(1);
+              {
+                cell->face(i)->set_boundary_id(1);
+                cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
+              }
           }
         ++cell;
       }
+    tria.set_manifold(0, SphericalManifold<2>(p));
   }
 
 
@@ -2632,6 +2652,8 @@ namespace GridGenerator
 
         tria.last()->face(1)->set_boundary_id(2);
       }
+    tria.set_all_manifold_ids(0);
+    tria.set_manifold(0, SphericalManifold<2>(center));
   }
 
 
@@ -2712,6 +2734,9 @@ namespace GridGenerator
 
         tria.last()->face(1)->set_boundary_id(2);
       }
+
+    tria.set_all_manifold_ids(0);
+    tria.set_manifold(0, SphericalManifold<2>(center));
   }
 
 
@@ -2896,6 +2921,7 @@ namespace GridGenerator
       }
 
     triangulation.create_triangulation (vertices, cells, SubCellData ());
+    triangulation.set_all_manifold_ids_on_boundary(0);
 
     for (Triangulation<3>::cell_iterator cell = triangulation.begin ();
          cell != triangulation.end (); ++cell)
@@ -2903,6 +2929,7 @@ namespace GridGenerator
         if (cell->vertex (0) (0) == -half_length)
           {
             cell->face (4)->set_boundary_id (1);
+            cell->face (4)->set_manifold_id(numbers::flat_manifold_id);
 
             for (unsigned int i = 0; i < 4; ++i)
               cell->line (i)->set_boundary_id (0);
@@ -2911,6 +2938,7 @@ namespace GridGenerator
         if (cell->vertex (4) (0) == half_length)
           {
             cell->face (5)->set_boundary_id (2);
+            cell->face (5)->set_manifold_id(numbers::flat_manifold_id);
 
             for (unsigned int i = 4; i < 8; ++i)
               cell->line (i)->set_boundary_id (0);
@@ -2919,6 +2947,8 @@ namespace GridGenerator
         for (unsigned int i = 0; i < 4; ++i)
           cell->face (i)->set_boundary_id (0);
       }
+
+    triangulation.set_manifold(0, CylindricalManifold<3>());
   }
 
 
@@ -3002,7 +3032,8 @@ namespace GridGenerator
   void
   hyper_ball (Triangulation<3> &tria,
               const Point<3>   &p,
-              const double radius)
+              const double radius,
+              const bool internal_manifold)
   {
     const double a = 1./(1+std::sqrt(3.0)); // equilibrate cell sizes at transition
     // from the inner part to the radial
@@ -3052,12 +3083,17 @@ namespace GridGenerator
         for (unsigned int j=0; j<GeometryInfo<3>::vertices_per_cell; ++j)
           cells[i].vertices[j] = cell_vertices[i][j];
         cells[i].material_id = 0;
+        cells[i].manifold_id = i == 0 ? numbers::flat_manifold_id : 1;
       };
 
     tria.create_triangulation (
       std::vector<Point<3> >(std::begin(vertices), std::end(vertices)),
       cells,
       SubCellData());       // no boundary information
+    tria.set_all_manifold_ids_on_boundary(0);
+    tria.set_manifold(0, SphericalManifold<3>(p));
+    if (internal_manifold)
+      tria.set_manifold(1, SphericalManifold<3>(p));
   }
 
 
@@ -3074,6 +3110,8 @@ namespace GridGenerator
     boundary_ids.insert (0);
     GridGenerator::extract_boundary_mesh (volume_mesh, tria,
                                           boundary_ids);
+    tria.set_all_manifold_ids(0);
+    tria.set_manifold(0, SphericalManifold<spacedim-1, spacedim>(p));
   }
 
 
@@ -3163,6 +3201,8 @@ namespace GridGenerator
     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=0; i<GeometryInfo<3>::faces_per_cell; ++i)
         if (cell->at_boundary(i))
@@ -3170,26 +3210,35 @@ namespace GridGenerator
             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_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_boundary_id(1);
+                      cell->face(i)->line(e)->set_manifold_id(numbers::flat_manifold_id);
+                    }
               }
           }
+    tria.set_manifold(0, CylindricalManifold<3>());
   }
 
 
@@ -3245,6 +3294,7 @@ namespace GridGenerator
     Triangulation<dim>::cell_iterator cell = tria.begin();
     Triangulation<dim>::cell_iterator end = tria.end();
 
+    tria.set_all_manifold_ids_on_boundary(0);
     while (cell != end)
       {
         for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
@@ -3258,6 +3308,7 @@ namespace GridGenerator
                 || cell->face(i)->center()(2) < center(2)+1.e-5 * radius)
               {
                 cell->face(i)->set_boundary_id(1);
+                cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
                 // also set the boundary indicators of the bounding lines,
                 // unless both vertices are on the perimeter
                 for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
@@ -3271,13 +3322,17 @@ namespace GridGenerator
                         ||
                         (std::fabs(line_vertices[1].distance(center)-radius) >
                          1e-5*radius))
-                      cell->face(i)->line(j)->set_boundary_id(1);
+                      {
+                        cell->face(i)->line(j)->set_boundary_id(1);
+                        cell->face(i)->line(j)->set_manifold_id(numbers::flat_manifold_id);
+                      }
                   }
 
               }
           }
         ++cell;
       }
+    tria.set_manifold(0, SphericalManifold<3>(center));
   }
 
 
@@ -3346,6 +3401,8 @@ namespace GridGenerator
     Triangulation<3>::cell_iterator cell = tria.begin();
     Triangulation<3>::cell_iterator end = tria.end();
 
+    tria.set_all_manifold_ids_on_boundary(0);
+
     // go over all faces. for the ones on the flat face, set boundary
     // indicator for face and edges to one; the rest will remain at
     // zero but we have to pay attention to those edges that are
@@ -3364,6 +3421,7 @@ namespace GridGenerator
             if (cell->face(i)->center()(0) < center(0)+1.e-5*radius)
               {
                 cell->face(i)->set_boundary_id(1);
+                cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
                 for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
                   {
                     const Point<3> line_vertices[2]
@@ -3375,12 +3433,16 @@ namespace GridGenerator
                         ||
                         (std::fabs(line_vertices[1].distance(center)-radius) >
                          1e-5*radius))
-                      cell->face(i)->line(j)->set_boundary_id(1);
+                      {
+                        cell->face(i)->line(j)->set_boundary_id(1);
+                        cell->face(i)->line(j)->set_manifold_id(numbers::flat_manifold_id);
+                      }
                   }
               }
           }
         ++cell;
       }
+    tria.set_manifold(0, SphericalManifold<3>(center));
   }
 
 
@@ -3608,6 +3670,8 @@ namespace GridGenerator
 
     if (colorize)
       colorize_hyper_shell(tria, p, inner_radius, outer_radius);
+    tria.set_all_manifold_ids(0);
+    tria.set_manifold(0, SphericalManifold<3>(p));
   }
 
 
@@ -3719,6 +3783,8 @@ namespace GridGenerator
                   }
               }
       }
+    tria.set_all_manifold_ids(0);
+    tria.set_manifold(0, SphericalManifold<3>(center));
   }
 
 
@@ -3783,6 +3849,9 @@ namespace GridGenerator
 
     if (colorize)
       colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius);
+
+    tria.set_all_manifold_ids(0);
+    tria.set_manifold(0, SphericalManifold<3>(center));
   }
 
 
@@ -3865,6 +3934,8 @@ namespace GridGenerator
 
     tria.create_triangulation (
       vertices_3d, cells, SubCellData());
+    tria.set_all_manifold_ids(0);
+    tria.set_manifold(0, CylindricalManifold<3>(2));
   }
 
 
@@ -4263,18 +4334,25 @@ namespace GridGenerator
                   else if (std::abs(dy - outer_radius) < eps)
                     cell->face(f)->set_boundary_id(3);
                   else
-                    cell->face(f)->set_boundary_id(4);
+                    {
+                      cell->face(f)->set_boundary_id(4);
+                      cell->face(f)->set_manifold_id(0);
+                    }
                 }
               else
                 {
                   double d = (cell->face(f)->center() - center).norm();
                   if (d-inner_radius < 0)
-                    cell->face(f)->set_boundary_id(1);
+                    {
+                      cell->face(f)->set_boundary_id(1);
+                      cell->face(f)->set_manifold_id(0);
+                    }
                   else
                     cell->face(f)->set_boundary_id(0);
                 }
             }
       }
+    triangulation.set_manifold(0, PolarManifold<2>(center));
   }
 
 
@@ -4373,9 +4451,8 @@ namespace GridGenerator
 
                   else
                     {
-                      cell->face(f)->set_boundary_id(6);
-                      for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
-                        cell->face(f)->line(l)->set_boundary_id(6);
+                      cell->face(f)->set_all_boundary_ids(6);
+                      cell->face(f)->set_all_manifold_ids(0);
                     }
 
                 }
@@ -4386,15 +4463,15 @@ namespace GridGenerator
                   double d = c.norm();
                   if (d-inner_radius < 0)
                     {
-                      cell->face(f)->set_boundary_id(1);
-                      for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
-                        cell->face(f)->line(l)->set_boundary_id(1);
+                      cell->face(f)->set_all_boundary_ids(1);
+                      cell->face(f)->set_all_manifold_ids(0);
                     }
                   else
                     cell->face(f)->set_boundary_id(0);
                 }
             }
       }
+    triangulation.set_manifold(0, CylindricalManifold<3>(2));
   }
 
   template <int dim, int spacedim1, int spacedim2>

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.