<h3>Working with manifolds</h3>
+<h4>What we want</h4>
+
In many applications, the finite element mesh must be able to represent a
relatively complex geometry. In the step-1, step-49, and step-53 tutorial
programs, some techniques to generate grids available within the deal.II
create a finer mesh by subdividing the cells into children, either uniformly
or only in selected parts of the computational domain. Besides the basic
meshing capabilities collected in the GridGenerator namespace, deal.II also
-comes with interfaces to read in meshes generated by (hex-only) mesh
-generators, as for example demonstrated in step-5. A fundamental limitation of
+comes with interfaces to read in meshes generated by (quad- and hex-only) mesh
+generators using the functions of namespace GridIn, as for example
+demonstrated in step-5. A fundamental limitation of
externally generated meshes is that the information provided by the generated
cells in the mesh only consists of the position of the vertices and their
connectivity, without the context of the underlying geometry that used to be
available in the mesh generator that originally created this mesh. This
becomes problematic once the mesh is refined within deal.II and additional
-points need to be placed. The step-54 tutorial program has shown how to
+points need to be placed. The step-54 tutorial program shows how to
overcome this limitation by using CAD surfaces in terms of the OpenCASCADE
-library.
+library, and step-53 by providing the same kind of information
+programmatically from within the source code.
Within deal.II, the placement of new points during mesh refinement or for the
definition of higher order mappings is controlled by manifold objects, see the
@ref manifold "manifold module"
for details.
-
To give an example, consider the following situation of a two-dimensional
annulus (with pictures taken from the manifold module). If we start with an
initial mesh of 10 cells and refine the mesh three times globally without
@image html hypershell-nothing.png ""
+The picture looks like this because, by default, deal.II only knows
+where to put the vertices of child cells by averaging the locations of
+the vertices of the parent cell. This yields a polygonal domain whose
+faces are the edges of the original (coarse mesh) cells.
Obviously, we must attach a curved description to the boundary faces of the
triangulation to reproduce the circular shape upon mesh refinement, like in
-the following picture
+the following picture:
@image html hypershell-boundary-only.png ""
+This is better: At least the inner and outer boundaries are now
+approaching real circles if we continue to refine the mesh.
However, the mesh in this picture is still not optimal for an annulus in the
-sense that the lines from one cell to the next have kinks at certain vertices,
+sense that the <i>interior</i> lines from one cell to the next have kinks at certain vertices,
and one would rather like to use the following mesh:
@image html hypershell-all.png ""
(see the @ref manifold "manifold module" for more details).
At this point, one might assume that curved volume descriptions are the way to
-go. However, this becomes impossible for as simple a case as the
-two-dimensional disk because the polar manifold degenerates in the origin and
-would not produce reasonable new points. A similar thing happens at the origin
-of the three-dimensional ball when one tries to attach a spherical manifold to
-the whole volume &endash; in this case, the computation of new manifold points
-would abort with an exception. These two simple examples make it clear that
+go. This is generally not wrong, though it is sometimes not so easy to
+describe how exactly this should work. Here are couple examples:
+
+- Imagine that the mesh above had actually been a disk, not just a ring.
+ In that case the polar manifold degenerates at the origin and
+ would not produce reasonable new points. In fact, defining a
+ manifold description for things that are supposed "to look round"
+ but might have points at or close to the origin is surprisingly very
+ difficult.
+- A similar thing happens at the origin
+ of the three-dimensional ball when one tries to attach a spherical manifold to
+ the whole volume &endash; in this case, the computation of new manifold points
+ would abort with an exception.
+- CAD geometries often only describe the boundary of the domain, in a
+ similar way to how we only attached a manifold to the boundary in
+ the second picture above. Similarly, step-54 only uses the CAD
+ geometry to generate a surface mesh (maybe because that is what is
+ needed to solve the problem in question), but if one wanted to solve
+ a problem in the water or the air around the ship described there,
+ we would need to have a volume mesh. The question is then how
+ exactly we should describe what is supposed to happen in the
+ interior of the domain.
+
+These simple examples make it clear that
for many interesting cases we must step back from the desire to have an
-analytic curved description for the full volume. This is particularly true if
-the boundary description is provided by some CAD files like in the step-54
-tutorial program, which are intrinsically surface-only.
-
-Yet, a curved boundary description alone is sometimes not enough. Consider the
+analytic curved description for the full volume: There will need to be
+<i>some</i> kind of information that leads to curvature also in the
+interior, but it must be possible to do this without actually writing
+down an explicit formula that describes the kind of geometry.
+
+So what happens if we don't do anything at all in the interior and
+only describe the surface as a manifold? Sometimes, as in the ring
+shown above, the result is not terrible. But sometimes it is. Consider the
case of a torus (e.g. generated with GridGenerator::torus()) with a
TorusManifold object attached to the surface only, no additional manifolds on
the interior cells and faces, and with six cells in toroidal direction before
@image html torus_no_inner_manifold.png ""
-This is clearly sub-optimal, and the mapping actually inverts in some regions
+This is clearly sub-optimal. Indeed, if we had started with fewer than
+the six cells shown above in toroidal direction, the mapping actually
+inverts in some regions
because the new points placed along interior cells intersect with the boundary
as they are not following the circular shape along the toroidal direction. The
simple case of a torus can still be fixed because we know that the toroidal
<h3>The class TransfiniteInterpolationManifold</h3>
In order to find a better strategy, let us look at the two-dimensional disk
-again (that is also the base entity revoluted along the torioidal direction in
+again (that is also the base entity rotated along the toroidal direction in
the torus). As we learned above, we can only apply the curved polar
description to the boundary (or a rim of cells sufficiently far away from the
origin) but must eventually transition to a straight description towards the
-disk's center. If we use a flat manifold in the interior of the cells and a
+disk's center. If we use a flat manifold in the interior of the cells
+(i.e., one in which new vertices are created by averaging of the
+adjacent existing ones) and a
polar manifold only for the boundary of the disk, we get the following mesh
-upon four global refinements
+upon four global refinements:
@image html circular_mesh_only_boundary_manifold.png ""
+If you know that the original coarse mesh consisted of a single square
+in the middle, with four caps around it, then it's not hard to see
+every refinement step that happened to this mesh to get the picture
+above.
+
While the triangulation class of deal.II tries to propagate information from
the boundary into the interior when creating new points, the reach of this
algorithm is limited:
manifold description.
At this point, we realize what would be needed to create a better mesh: For
-<b>all</b> new points in any child cell that is created within the red shaded
+<i>all</i> new points in any child cell that is created within the red shaded
layer on the leftmost picture, we want to compute the interpolation with
respect to the curvature in the area covered by the respective coarse
cell. This is achieved by adding the class TransfiniteInterpolationManifold to
set of points within its domain of definition, it can compute weighted
averages conforming to the manifold (using a formula that will be given in a
minute). These weighted averages are used whenever the mesh is refined, or
-when a higher order mapping, MappingQGeneric, is evaluated on a given cell
+when a higher order mapping (such as MappingQGeneric or MappingC1)
+is evaluated on a given cell
subject to this manifold. Using this manifold on the shaded cells of the
coarse grid with of the disk produces the following mesh upon four global
steps of refinement:
@image html circular_mesh_transfinite_interpolation.png ""
-Given a straight-sided central cell, this representation is the best possible
+There are still some kinks in the lines of this mesh, but they are
+restricted to the faces between coarse mesh cells, whereas the rest of
+the mesh is about as smooth as one would like. Indeed,
+given a straight-sided central cell, this representation is the best possible
one as all mesh cells follow a smooth transition from the straight sides in
the square block in the interior to the circular shape on the boundary. (One
could possibly do a bit better by allowing some curvature also in the central
square block, that eventually vanishes as the center is approached.)
+
+<h4>How it works</h4>
+
In the simple case of a disk with one curved and three straight edges, we can
explicitly write down how to achieve the blending of the shapes. For this, it
is useful to map the physical cell, like the top one, back to the reference
interpolation in the second line of the formula makes sure that the prescribed
curves are followed exactly on the boundary: Along each of the four edges, we
need to subtract the contribution of the two adjacent edges evaluated in the
-corners, which is nothing else than a vertex position. It is easy to check
+corners, which is then simply a vertex position. It is easy to check
that the formula for the circle above is reproduced if three of the four
curves $\mathbf{c}_i$ are straight and thus coincide with the bilinear
interpolation.
The interpolation of the transfinite interpolation in deal.II is general in
the sense that it can deal with arbitrary curves. It will evaluate the curves
in terms of their original coordinates of the $d$-dimensional space but with
-one (or two for edges in 3D) coordinate held fixed at $0$ or $1$ to ensure
+one (or two, in the case of edges in 3D) coordinate held fixed at $0$ or $1$ to ensure
that any other manifold class, including CAD files if desired, can be applied
out of the box. Transfinite interpolation is a standard ingredient in mesh
generators, so the main strength of the integration of this feature within the
As a final remark on transfinite interpolation, we mention that the mesh
refinement strategies in deal.II in absence of a volume manifold description
are also based on the weights of the transfinite interpolation and optimal in
-that sense. As mentioned above, this is however limited to operations on those
-cells touching the curved manifolds.
+that sense. The difference is that the default algorithm sees only one
+cell at a time, and so will apply the optimal algorithm only on those
+cells touching the curved manifolds. In contrast, using the
+transfinite mapping on entire <i>patches</i> of cells (originating
+from one coarser cell) allows to use the transfinite interpolation
+method in a way that propagates information from the boundary to cells
+far away.
+
<h3>Transfinite interpolation is expensive and how to deal with it</h3>
Manifold interface class of deal.II only provides a set of points, the
transfinite interpolation initially does not even know to which coarse grid
cell the set of surrounding points belong to and needs to search among several
-cells based on some heuristics. In terms of charts, one could describe the
-implementation of the transfinite interpolation as an atlas-based
+cells based on some heuristics. In terms of <a
+href="https://en.wikipedia.org/wiki/Atlas_(topology)#Charts">charts</a>,
+one could describe the
+implementation of the transfinite interpolation as an <a
+href="https://en.wikipedia.org/wiki/Atlas_(topology)">atlas</a>-based
implementation: Each cell of the initial coarse grid of the triangulation
represents a chart with its own reference space, and the surrounding manifolds
provide a way to transform from the chart space (i.e., the reference cell) to
numerical solution. If the degree of the geometry is higher or lower than the
solution, one calls that a super- or sub-parametric geometry representation,
respectively. In deal.II, the standard class for polynomial representation is
-MappingQGeneric. If this class is used with polynomial degree $4$ in 3D, a
+MappingQGeneric. If, for example, this class is used with polynomial degree $4$ in 3D, a
total of 125 ($=(4+1)^3$) points are needed for the tri-cubic
interpolation. Among these points, 8 are the mesh vertices and already
available from the mesh, but the other 117 need to be provided by the
yet to be determined coarse cell, followed by subsequent push-forward on each
of the 117 points, is a lot of work and can be very time consuming.
-What makes things worse is that the structure of many programs is that the
+What makes things worse is that the structure of many programs is such
+that the
mapping is queried several times independently for the same cell. Its primary
use is in the assembly of the linear system, i.e., the computation of the
system matrix and the right hand side, via the `mapping` argument of the
MappingQCache does not help and there currently does not exist a more
efficient mechanism in deal.II. However, the mesh refinement contains many
other expensive steps as well, so it is not as big as an issue compared to the
-rest of the computation.
+rest of the computation. It also only happens at most once per time
+step or nonlinear iteration.
<h3>The test case</h3>
projects / dealii.git / commitdiff