From: David Wells Date: Sun, 4 Feb 2018 21:11:23 +0000 (-0500) Subject: Redo the step-49 extensions to use Manifolds. X-Git-Tag: v9.0.0-rc1~450^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b3360284f188ff0f9cac8a7f968f7fb233329e80;p=dealii.git Redo the step-49 extensions to use Manifolds. The old approach used the deprecated Boundary classes. --- diff --git a/examples/step-49/doc/results.dox b/examples/step-49/doc/results.dox index a392ccc3aa..f1b381bf1a 100644 --- a/examples/step-49/doc/results.dox +++ b/examples/step-49/doc/results.dox @@ -4,41 +4,55 @@ The program produces a series of .eps files of the triangulations. The methods are discussed above. -

Next steps: Curved boundaries

- -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: +

Next steps: Curved Cells

+ +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: -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: 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 + manifold_id 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 + manifold_id. 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. @@ -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: - - - -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: + + + +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: - + -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 manifold_ids 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: - - - -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 faces but not for the adjacent edges 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: - - - -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::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 - +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: + - +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! +