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&M University.
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,
}
@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">
projects / dealii.git / commitdiff