]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Redo the step-49 extensions to use Manifolds. 5860/head
authorDavid Wells <wellsd2@rpi.edu>
Sun, 4 Feb 2018 21:11:23 +0000 (16:11 -0500)
committerDavid Wells <wellsd2@rpi.edu>
Sat, 10 Feb 2018 12:40:03 +0000 (07:40 -0500)
The old approach used the deprecated Boundary classes.

examples/step-49/doc/results.dox

index a392ccc3aa805206284704daae91943480d279b7..f1b381bf1a4d35fc3cdcc1ec0f3a7a52c086b581 100644 (file)
@@ -4,41 +4,55 @@ The program produces a series of <code>.eps</code> files of the
 triangulations. The methods are discussed above.
 
 
-<h3>Next steps: Curved boundaries</h3>
-
-As mentioned in the introduction,
-creating a coarse mesh using the methods discussed here is only the first
-step. In order to refine a mesh, the Triangulation needs to know where to put
-new vertices on the mid-points of edges and faces. By default, these new
-points will be placed at the centers of the old edge but this isn't what you
-want if you need curved boundaries that aren't already adequately resolved by
-the coarse mesh. Several of the meshes shown in the introduction section fall
-into this category. For example, for this mesh the central hole is supposed to
-be round:
+<h3>Next steps: Curved Cells</h3>
+
+As mentioned in the introduction, creating a coarse mesh using the methods
+discussed here is only the first step. In order to refine a mesh, the
+Triangulation needs to know where to put new vertices on the mid-points of
+edges, faces, and cells. By default, these new points will be placed at the
+arithmetic mean of the surrounding points, but this isn't what you want if you
+need curved boundaries that aren't already adequately resolved by the coarse
+mesh. Several of the meshes shown in the introduction section fall into this
+category. For example, for this mesh the central hole is supposed to be round:
 
 <img src="https://www.dealii.org/images/steps/developer/step-49.grid-2a.png" alt="" height="200px">
 
-On the other hand, if you simply refine it, the Triangulation class can not
-know whether you wanted the hole to be round or to be an octagon. The default
-is to place new points along existing edges. After two mesh refinement steps,
-this would yield the following mesh, which is not what we wanted:
+If you simply refine it, the Triangulation class can not know whether you wanted
+the hole to be round or to be an octagon. The default is to place new points
+along existing straight lines. After two mesh refinement steps, this would yield
+the following mesh, which is not what we wanted:
 
 <img src="https://www.dealii.org/images/steps/developer/step-49.grid-2d-refined.png" alt="" height="200px">
 
 What needs to happen is that you tell the triangulation that you in fact want
-to use a curved boundary. The way to do this requires three steps:
-- Create an object that describes the boundary in terms that allow the
-  Triangulation::execute_coarsening_and_refinement function to ask the
-  boundary description where a new point should be located upon mesh
-  refinement.
-- Tell the triangulation object that you want this object to be used for all
-  boundaries with boundary indicates equal to a particular value (for more
-  information on boundary indicators, see the
-  @ref GlossBoundaryIndicator "glossary entry on this topic".)
-- Mark those parts of the boundary of the domain for which you want the
-  boundary to be so treated with the value of the boundary indicator used in
-  the previous step. (The order of this step and the previous one does not
-  matter.)
+to use a curved geometry. The way to do this requires three steps:
+- Create an object that describes the desired geometry. This object will be
+  queried when refining the Triangulation for new point placement. It will also
+  be used to calculate shape function values if a high degree mapping, like
+  MappingQ or MappingQGeneric, is used during system assembly.
+  In deal.II the Manifold class and classes inheriting from it (e.g.,
+  PolarManifold and FlatManifold) perform these calculations.
+- Notify the Triangulation object which Manifold classes to use. By default, a
+  Triangulation uses FlatManifold to do all geometric calculations,
+  which assumes that all cell edges are straight lines and all quadrilaterals
+  are flat. You can attach Manifold classes to a Triangulation by calling
+  Triangulation::set_manifold function, which associates a
+  <code>manifold_id</code> with a Manifold object. For more information on this
+  see the @ref GlossManifoldIndicator "glossary entry on this topic".
+- Finally, you must mark cells and cell faces with the correct
+  <code>manifold_id</code>. For example, you could get an annular domain with
+  curved cells in Cartesian coordinates (but rectangles in polar coordinates)
+  by doing the following:
+  @code
+  PolarManifold<2> polar_manifold;
+  Triangulation<2> triangulation;
+  const types::manifold_id polar_id = 42;
+  GridGenerator::hyper_shell(triangulation, Point<2>(), 0.5, 1.0, 10);
+  triangulation.set_manifold(polar_id, polar_manifold);
+  triangulation.set_all_manifold_ids(polar_id);
+  @endcode
+  Now, when the grid is refined, all cell splitting calculations will be done in
+  polar coordinates.
 
 To illustrate this process in more detail, let us consider an example created
 by Yuhan Zhou as part of a 2013 semester project at Texas A&amp;M University.
@@ -52,10 +66,10 @@ In the following, we will walk you through the entire process of creating a
 mesh for this geometry, including a number of common pitfalls by showing the
 things that can go wrong.
 
-The first step in getting there was to create a coarse mesh, which was done
-by creating a 2d coarse mesh for each of the two cross section, extruding them
-into the third direction, and gluing them together. The following code does
-this, using the techniques previously described:
+The first step in getting there was to create a coarse mesh, which was done by
+creating a 2d coarse mesh for each of cross sections, extruding them into the
+third direction, and gluing them together. The following code does this, using
+the techniques previously described:
 
 @code
 // Given a list of points and how vertices connect to cells,
@@ -225,158 +239,109 @@ void create_3d_grid (Triangulation<3> &triangulation)
 }
 @endcode
 
-With this code, you get a mesh that looks like this:
-
-<img src="https://www.dealii.org/images/steps/developer/step-49.yuhan.2.png" alt="">
-
-The next step is to teach each of the top surfaces that they should be
-curved. We can do this by creating CylinderBoundary objects that
-describe this. A first attempt looks like this:
+This creates the following mesh:
+
+<img src="https://www.dealii.org/images/steps/developer/step-49.yuhan.8.png"
+     alt="" width="400" height="355">
+
+This mesh has the right general shape, but the top cells are now polygonal: their
+edges are no longer along circles and we do not have a very accurate
+representation of the original geometry. The next step is to teach the top part
+of the domain that it should be curved. Put another way, all calculations done
+on the top boundary cells should be done in cylindrical coordinates rather than
+Cartesian coordinates. We can do this by creating a CylindricalManifold object
+and associating it with the cells above $y = 3$. This way, when we refine the
+cells on top, we will place new points along concentric circles instead of
+straight lines.
+
+In deal.II we describe all geometries with classes that inherit from
+Manifold. The default geometry is Cartesian and is implemented in the
+FlatManifold class. As the name suggests, Manifold and its inheriting classes
+provide a way to describe curves and curved cells in a general way with ideas
+and terminology from differential geometry: for example, CylindricalManifold
+inherits from ChartManifold, which describes a geometry through pull backs
+and push forwards. In general, one should think that the Triangulation class
+describes the topology of a domain (in addition, of course, to storing the
+locations of the vertices) while the Manifold classes describe the geometry of a
+domain (e.g., whether or not a pair of vertices lie along a circular arc or a
+straight line). A Triangulation will refine cells by doing computations with the
+Manifold associated with that cell regardless of whether or not the cell is on
+the boundary. Put another way: the Manifold classes do not need any information
+about where the boundary of the Triangulation actually is: it is up to the
+Triangulation to query the right Manifold for calculations on a cell. Most
+Manifold functions (e.g., Manifold::get_intermediate_point) know nothing about
+the domain itself and just assume that the points given to it lie along a
+geodesic. In this case, with the CylindricalManifold constructed below, the
+geodesics are arcs along circles orthogonal to the $z$-axis centered along the
+line $(0, 3, z)$.
+
+Since all three top parts of the domain use the same geodesics, we will
+mark all cells with centers above the $y = 3$ line as being cylindrical in
+nature:
 
 @code
+  const Tensor<1, 3> axis({0.0, 0.0, 1.0});
+  const Point<3> axial_point(0, 3.0, 0.0);
+  const CylindricalManifold<3> cylinder(axis, axial_point);
+  const types::manifold_id cylinder_id = 8;
+
   Triangulation<3> triangulation;
-  create_3d_grid (triangulation);
-
-  // Create the objects that describe the boundaries and attach them
-  // to the triangulation as the ones to use on boundaries marked
-  // with boundary indicators 8 and 9
-  const double inner_radius = 1.5;
-  const double outer_radius = 2.5;
-
-  static const CylinderBoundary<3> inner_cylinder(inner_radius, 2);
-  static const CylinderBoundary<3> outer_cylinder(outer_radius, 2);
-
-  triangulation.set_boundary (8, inner_cylinder);
-  triangulation.set_boundary (9, outer_cylinder);
-
-  // Then loop over all faces of the domain and, if for the position
-  // of the center of a face the following holds then set boundary
-  // indicators:
-  // - if y>3 and z<=2.5 or z>=5 then use boundary indicator 8
-  // - if y>3 and 2.5<=z<=5 then use boundary indicator 9
-  typename Triangulation<3>::active_cell_iterator
-    cell = triangulation.begin_active(),
-    endc = triangulation.end();
-  for (; cell!=endc; ++cell)
-    for (unsigned int f=0;
-        f < GeometryInfo<3>::faces_per_cell;
-        ++f)
-      {
-       const Point<3> face_center = cell->face(f)->center();
-
-       if (cell->face(f)->at_boundary())
-         {
-           if ((face_center[2] <= 2.5 || face_center[2] >= 5) &&
-               face_center[1] >= 3)
-             cell->face(f)->set_boundary_id(8);
-
-           if (face_center[2] >= 2.5 &&
-               face_center[2] <= 5
-               && face_center[1] >= 3)
-             cell->face(f)->set_boundary_id(9);
-         }
-      }
+  create_3d_grid(triangulation);
+  triangulation.set_manifold (cylinder_id, cylinder);
+
+  for (auto cell : triangulation.active_cell_iterators())
+    if (cell->center()[1] >= 3.0)
+      cell->set_all_manifold_ids(cylinder_id);
 
-  // Then refine the mesh once
   triangulation.refine_global (1);
 @endcode
 
 With this code, we get a mesh that looks like this:
 
-<img src="https://www.dealii.org/images/steps/developer/step-49.yuhan.3.png" alt="">
+<img src="https://www.dealii.org/images/steps/developer/step-49.yuhan.9.png"
+     alt="" width="400" height="355">
 
-This is clearly not correct: The new vertices that have been entered at
-mid-edge and mid-face points are not where they should have been. Upon some
-reflection, it turns out that while the radii of the cylinders are correct,
-the axes of the two cylinder objects should not have been along coordinate
-axes but shifted. This can be corrected by creating them as follows, the
-two points given as arguments indicating the direction and a point on the
-axis:
+This change fixes the boundary but creates a new problem: the cells adjacent to
+the cylinder's axis are badly distorted. We should use Cartesian coordinates for
+calculations on these central cells to avoid this issue. The cells along the
+center line all have a face that touches the line $(0, 3, z)$ so, to implement
+this, we go back and overwrite the <code>manifold_id</code>s on these cells to
+be zero (which is the default):
 
 @code
-  static const CylinderBoundary<3> inner_cylinder (inner_radius,
-                                                   Point<3>(0,0,1),
-                                                   Point<3>(0,3,0));
-  static const CylinderBoundary<3> outer_cylinder (outer_radius,
-                                                   Point<3>(0,0,1),
-                                                   Point<3>(0,3,0));
-  triangulation.set_boundary (9, outer_cylinder);
-@endcode
+  const Tensor<1, 3> axis({0.0, 0.0, 1.0});
+  const Point<3> axial_point(0, 3.0, 0.0);
+  const CylindricalManifold<3> cylinder(axis, axial_point);
+  const types::manifold_id cylinder_id = 8;
 
-This yields an improvement, though it is still not quite correct:
-
-<img src="https://www.dealii.org/images/steps/developer/step-49.yuhan.4.png" alt="">
-
-Looking closely at this mesh, we realize that the new points on mid-face
-vertices are where they should be, though the new vertices inserted at
-mid-edge points are in the wrong place (you see this by comparing the
-picture with the one of the coarse mesh). What is happening is that we
-are only telling the triangulation to use these geometry objects for
-the <i>faces</i> but not for the adjacent <i>edges</i> as well. This is
-easily fixed by using the function TriaAccessor::set_all_boundary_ids()
-instead of TriaAccessor::set_boundary_id() used above. With this change,
-the grid now looks like this:
-
-<img src="https://www.dealii.org/images/steps/developer/step-49.yuhan.5.png" alt="">
-
-This is already better. However, something is still going wrong on the
-front left face. On second look, we can also see that the faces where
-the geometry widens have been refined at the bottom, that there is one
-transition face that looks wrong because it has a triangle rather than
-a quadrilateral, and that finally the transition faces in the cylindrical
-region appear not to have been refined at all in radial direction.
-
-This due to the fact that we have (erroneously) marked all boundary faces
-between $0\le z\le 2.5$ with the boundary indicator for the small cylinder
-and similarly for the other regions. This condition includes the faces parallel
-to the x-y plane. To fix it, we need to exclude faces whose center points have
-$z$ values equal to (or at least close to, since we should not compare
-for equality in floating point arithmetic) 0, 2.5, 5 or 7.5. This replacement
-code does the trick:
+  Triangulation<3> triangulation;
+  create_3d_grid(triangulation);
+  triangulation.set_manifold (cylinder_id, cylinder);
 
-@code
-  // Then loop over all faces of the domain and, if for the position
-  // of the center of a face the following holds then set boundary
-  // indicators:
-  // - if y>3 and z<2.5 or z>5 then use boundary indicator 8
-  // - if y>3 and 2.5<z<5 then use boundary indicator 9
-  // In this process, exclude faces whose z-coordinates are
-  // within a small distance of z=0, z=2.5, z=5 or z=7.5.
-  typename Triangulation<3>::active_cell_iterator
-    cell = triangulation.begin_active(),
-    endc = triangulation.end();
-  for (; cell!=endc; ++cell)
-    for (unsigned int f=0;
-         f < GeometryInfo<3>::faces_per_cell;
-         ++f)
+  for (auto cell : triangulation.active_cell_iterators())
+    if (cell->center()[1] >= 3.0)
+      cell->set_all_manifold_ids(cylinder_id);
+
+  for (auto cell : triangulation.active_cell_iterators())
+    for (unsigned int face_n = 0; face_n < GeometryInfo<3>::faces_per_cell; ++face_n)
       {
-        const Point<3> face_center = cell->face(f)->center();
-
-        if (cell->face(f)->at_boundary())
-          if ((face_center[2]>1e-6) &&
-              (face_center[2]<7.5-1e-6) &&
-              (std::fabs(face_center[2]-2.5)>1e-6) &&
-              (std::fabs(face_center[2]-5.0)>1e-6))
-            {
-              if ((face_center[2] < 2.5 || face_center[2] > 5)
-                  && face_center[1] >= 3)
-                cell->face(f)->set_all_boundary_ids(8);
-
-              if (face_center[2] > 2.5 && face_center[2] < 5
-                  && face_center[1] >= 3)
-                cell->face(f)->set_all_boundary_ids(9);
-            }
+        const Point<3> face_center = cell->face(face_n)->center();
+        if (std::abs(face_center[0]) < 1.0e-5 &&
+            std::abs(face_center[1] - 3.0) < 1.0e-5)
+          cell->set_all_manifold_ids(numbers::flat_manifold_id);
       }
-}
-@endcode
 
-With this, we finally get a mesh that looks good:
+  triangulation.refine_global (1);
+@endcode
 
-<img src="https://www.dealii.org/images/steps/developer/step-49.yuhan.6.png" alt="">
+This gives us the following grid:
 
-We can then refine the mesh two more times to see in more detail what
-happens to the curved part of the boundary:
+<img src="https://www.dealii.org/images/steps/developer/step-49.yuhan.10.png"
+     alt="" width="400" height="355">
 
- <img src="https://www.dealii.org/images/steps/developer/step-49.yuhan.7.png" alt="">
+This gives us a good mesh, where cells at the center of each circle are still
+Cartesian and cells around the boundary lie along a circle. We can really see
+the nice detail of the boundary fitted mesh if we refine two more times:
 
- So, yes!, this is finally what we were looking for!
+<img src="https://www.dealii.org/images/steps/developer/step-49.yuhan.11.png"
+     alt="" width="400" height="355">

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.