<a name="Intro"></a>
<h1>Introduction</h1>
-This tutorial program presents a way to handle expensive manifold classes,
-exemplified for the TransfiniteInterpolationManifold, a special manifold class
-that can blend a curved approximation on a boundary with straight-sided
-descriptions elsewhere.
+This tutorial program presents an advanced manifold class,
+TransfiniteInterpolationManifold, and how to work around its main
+disadvantage, the relatively high cost.
+
+<h3>Working with manifolds</h3>
+
+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
+library have been introduced. Given a base mesh, deal.II is then able to
+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
+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
+overcome this limitation by using CAD surfaces in terms of the OpenCASCADE
+library.
+
+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
+attaching any manifolds, we would obtain the following mesh:
+
+@image html hypershell-nothing.png ""
+
+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
+
+@image html hypershell-boundary-only.png ""
+
+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,
+and one would rather like to use the following mesh:
+
+@image html hypershell-all.png ""
+
+In this last (optimal) case, which is also the default produced by
+GridGenerator::hyper_shell(), the curved manifold description (in this case a
+polar manifold description) is applied not only to the boundary faces, but to
+the whole domain. Whenever the triangulation requests a new point, e.g., the
+mid point of the edges or the cells when it refines a cell into four children,
+it will place them along the respective mid points in the polar coordinate
+system. By contrast, the case above where only the boundary was subject to the
+polar manifold, only mid points along the boundary would be placed along the
+curved description, whereas mid points in the interior would be computed by
+suitable averages of the surrounding points in the Cartesian coordinate system
+(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
+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
+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
+refinement. If the mesh is refined once, we would obtain the following mesh,
+shown with the upper half of the mesh clipped away:
+
+@image html torus_no_inner_manifold.png ""
+
+This is clearly sub-optimal, and 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
+direction follows a cylindrical coordindate system, so attaching a
+TorusManifold to the surface combined with CylindricalManifold with
+appropriate periodicity in toroidal direction applied to all interior entities
+would produce a high-quality mesh as follows, now shown with two top cells
+hidden:
+
+@image html torus_cylindrical_inner_manifold.png ""
+
+This mesh is pretty good, but obviously it is linked to a good description of
+the volume, which we lack in other cases. Actually, there is an imperfection
+also in this case, as we can see some unnatural kinks of two adjacent cells in
+the interior of the domain which are hidden by the top two boundary cells, as
+opposed to the following setup (the default manifolds applied by
+GridGenerator::torus() and using the TransfiniteInterpolationManifold):
+
+@image html torus_transfinite_manifold.png ""
+
+<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
+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
+polar manifold only for the boundary of the disk, we get the following mesh
+upon four global refinements
+
+@image html circular_mesh_only_boundary_manifold.png ""
+
+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:
+
+@image html circular_mesh_boundary_cells.png ""
+
+The picture above highlights those cells on the disk that are touching the
+boundary and where boundary information could in principle be taken into
+account when only looking at a single cell at the time. Clearly, the area
+where some curvature can be taken into account gets more limited as the mesh
+is refined, thus creating the seemingly irregular spots in the mesh: When
+computing the center of any one of the boundary cells in the leftmost picture,
+the ideal position is the mid point between the outer circle and the cell in
+the middle. This is exactly what is used for the first refinement step in the
+Triangulation class. However, for the second refinement all interior edges as
+well as the interior cell layers can only add points according to a flat
+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
+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
+the highlighted cells of the coarse grid in the leftmost panel of the figure
+above. This class adheres to the general manifold interfaces, i.e., given any
+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
+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
+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.)
+
+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
+coordinate system $(\xi,\eta)\in (0,1)^2$ where we compute averages between
+certain points. If we were to use a simple bilinear map spanned by four
+vertices $(x_0,y_0), (x_1,y_1), (x_2,y_2), (x_3, y_3)$, the image of a point
+$(\xi, \eta)\in (0,1)^2$ would be
+@f{align*}{
+(x,y) = (1-\xi)(1-\eta) (x_0,y_0) + \xi(1-\eta) (x_1,y_1) +
+ (1-\xi)\eta (x_2,y_2) + \xi\eta (x_3,y_3).
+@f}
+
+For the case of the curved surface, we want to modify this formula. For the
+top cell of the coarse mesh of the disk, we can assume that the points
+$(x_0,y_0)$ and $(x_1,y_1)$ sit along the straight line at the lower end and
+the points $(x_2,y_2)$ and $(x_3,y_3)$ are connected by a quarter circle along
+the top. We would then map a point $(\xi, \eta)$ as
+@f{align*}{
+(x,y) = (1-\eta) \big[(1-\xi) (x_0,y_0) + \xi (x_1,y_1)\big] +
+ \eta \mathbf{c}_3(\xi),
+@f}
+where $\mathbf{c}_3(\xi)$ is a curve that describes the $(x,y)$ coordinates of
+the quarter circle in terms of an arclength parameter $\xi\in (0,1)$. This
+represents a linear interpolation between the straight lower edge and the
+curved upper edge of the cell, and is the basis for the picture shown above.
+
+This formula is easily generalized to the case where all four edges are
+described by a curve rather than a straight line. We call the four functions,
+parameterized by a single coordinate $\xi$ or $\eta$ in the horizontal and
+vertical directions, $\mathbf{c}_0, \mathbf{c}_1, \mathbf{c}_2,
+\mathbf{c}_3$ for the left, right, lower, and upper edge of a
+quadrilateral, respectively. The interpolation then reads
+@f{align*}{
+(x,y) =& (1-\xi)\mathbf{c}_0(\eta) + \xi \mathbf{c}_1(\eta)
+ +(1-\eta)\mathbf{c}_2(\xi) + \eta \mathbf{c}_3(\xi)\\
+ &-\big[(1-\xi)(1-\eta) (x_0,y_0) + \xi(1-\eta) (x_1,y_1) +
+ (1-\xi)\eta (x_2,y_2) + \xi\eta (x_3,y_3)\big].
+@f}
+
+This formula assumes that the boundary curves match and coincide with the
+vertices $(x_0,y_0), (x_1,y_1), (x_2,y_2), (x_3, y_3)$, e.g. $\mathbf{c}_0(0)
+= (x_0,y_0)$ or $\mathbf{c}_0(1) = (x_2,y_2)$. The subtraction of the bilinear
+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
+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.
+
+This formula, called transfinite interpolation, was introduced in 1973 by <a
+href="https://doi.org/10.1002%2Fnme.1620070405">Gordon and Hall</a>. Even
+though transfinite interpolation essentially only represents a linear blending
+of the bounding curves, the interpolation exactly follows the boundary curves
+for each real number $\xi\in (0,1)$ or $\eta\in (0,1)$, i.e., it interpolates
+in an infinite number of points, which was the original motivation to label
+this variant of interpolation a transfinite one by Gordon and Hall. Another
+interpretation is that the transfinite interpolation interpolates from the
+left and right and the top and bottom linearly, from which we need to subtract
+the bilinear interpolation to ensure a unit weight in the interior of the
+domain.
+
+The transfinite interpolation is easily generalized to three spatial
+dimensions. In that case, the interpolation allows to blend 6 different
+surface descriptions for any of the quads of a three-dimensional cell and 12
+edge descriptions for the lines of a cell. Again, to ensure a consistent map,
+it is necessary to subtract the contribution of edges and add the contribution
+of vertices again to make the curves follow the prescribed surface or edge
+description. In the three-dimensional case, it is also possible to use a
+transfinite interpolation from a curved edge both into the adjacent faces and
+the adjacent cells.
+
+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
+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
+deal.II library is to enable it during adaptive refinement and coarsening of
+the mesh, and for creating higher-degree mappings that use manifolds to insert
+additional points beyond the mesh vertices.
+
+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.
+
+<h3>Transfinite interpolation is expensive and how to deal with it</h3>
+
+A mesh with a transfinite manifold description is typically set up in two
+steps. The first step is to create a mesh (or read it in from a file) and to
+attach a curved manifold to some of the mesh entities. For the above example
+of the disk, we attach a polar manifold to the faces along the outer circle
+(this is done automatically by GridGenerator::hyper_ball()). Before we start
+refining the mesh, we then assign a TransfiniteInterpolationManifold to all
+interior cells and edges of the mesh, which of course needs to be based on
+some manifold id that we have assigned to those entities (everything except
+the circle on the boundary). It does not matter whether we also assign a
+TransfiniteInterpolationManifold to the inner square of the disk or not
+because the transfinite interpolation becomes a flat representation to cells
+where all surrounding objects are also flat (or a transfinite interpolation
+with flat sub-entities).
+
+Later, when the mesh is refined or when a higher-order mapping is set up based
+on this mesh, the cells will query the underlying manifold object for new
+points. This process takes a set of surrounding points, for example the four
+vertices of a two-dimensional cell, and a set of weights to each of these
+points, for definition a new point. For the mid point of a cell, each of the
+four vertices would get weight 0.25. For the transfinite interpolation
+manifold, the process of building weighted sums requires some serious work. By
+construction, we want to combine the points in terms of the reference
+coordinates $\xi$ and $\eta$ (or $\xi, \eta, \zeta$ in 3D) of the surrounding
+points. However, the interface of the manifold classes in deal.II does not get
+the reference coordinates of the surrounding points (as they are not stored
+globally) but rather the physical coordinates only. Thus, the first step the
+transfinite interpolation manifold has to do is to invert the mapping and find
+the reference coordinates within one of the coarse cells of the transfinite
+interpolation (e.g. one of the four shaded coarse-grid cells of the disk mesh
+above). This inversion is done by a Newton iteration (or rather,
+finite-difference based Newton scheme combined with Broyden's method) and
+queries the transfinite interpolation according to the formula above several
+times. Each of these queries in turn might call an expensive manifold, e.g. a
+spherical description of a ball, and be expensive on its own. Since the
+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
+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
+the physical space. The collection of the charts of the coarse grid cells is
+an atlas, and as usual, the first thing one does when looking up something in
+an atlas is to find the right chart.
+
+Once the reference coordinates of the surrounding points have been found, a
+new point in the reference coordinate system is computed by a simple weighted
+sum. Finally, the reference point is inserted into the formula for the
+transfinite interpolation, which gives the desired new point.
+
+In a number of cases, the curved manifold is not only used during mesh
+refinement, but also to ensure a curved representation of boundaries within
+the cells of the computational domain. This is a necessity to guarantee
+high-order convergence for high-order polynomials on complex geometries
+anyway, but sometimes an accurate geometry is also desired with linear shape
+functions. This is often done by polynomial descriptions of the cells and
+called the isoparametric concept if the polynomial degree to represent the
+curved mesh elements is the same as the degree of the polynomials for the
+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
+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
+manifold. In case the transfinite interpolation manifold is used, we can
+imagine that going through the pull-back into reference coordinates of some
+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
+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
+FEValues object. However, also the interpolation of boundary values, the
+computation of numerical errors, writing the output, and evaluation of error
+estimators must involve the same mapping to ensure a consistent interpretation
+of the solution vectors. Thus, even a linear stationary problem that is solved
+once will evaluate the points of the mapping several times. For the cubic case
+in 3D mentioned above, this means computing 117 points per cell by an
+expensive algorithm many times. The situation is more pressing for nonlinear
+or time-dependent problems where those operations are done over and over
+again.
+
+As the manifold description via a transfinite interpolation can easily be
+hundreds of times more expensive than a similar query on a flat manifold, it
+makes sense to compute the additional points only once and use them in all
+subsequent calls. The deal.II library provides the class MappingQCache for
+exactly this purpose. The cache is typically not overly big compared to the
+memory consumed by a system matrix, as will become clear when looking at the
+results of this tutorial program. The usage of MappingQCache is simple: Once
+the mesh has been set up (or changed during refinement), we call
+MappingQCache::initialize() with the desired triangulation as well as a
+desired mapping as arguments. The initialization then goes through all cells
+of the mesh and queries the given mapping for its additional points. Those get
+stored for an identifier of the cell so that they can later be returned
+whenever the mapping computes some quantities related to the cell (like the
+Jacobians of the map between the reference and physical coordinates).
+
+As a final note, we mention that the TransfiniteInterpolationManifold also
+makes the refinement of the mesh more expensive. In this case, 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.
+
+<h3>The test case</h3>
+
+In this tutorial program, the usage of TransfiniteInterpolationManifold is
+exemplified in combination with MappingQCache. The test case is relatively
+simple and takes up the solution stages involved in many typical programs,
+e.g., the step-6 tutorial program. As a geometry, we select one prototype use
+of TransfiniteInterpolationManifold, namely a setup involving a spherical ball
+that is in turn surrounded by a cube. Such a setup is used for example in case
+that a material interface is located at the boundary of the ball within the
+computational domain that should be tracked by an element interface. A
+visualization of the grid is given here:
+
+<img src="https://www.dealii.org/images/steps/developer/step-65-mesh.png" alt="">
+
+For this case, we want to attach a spherical description to the surface inside
+the domain and use the transfinite interpolation to smoothly switch to the
+straight lines of the outer cube and the cube at the center of the ball.
+
+Within the program, we will follow a typical flow in finite element programs,
+starting from the setup of DoFHandler and sparsity patterns, the assembly of a
+linear system for solving the Poisson equation with a jumping coefficient, its
+solution with a simple iterative method, computation of some numerical error
+with VectorTools::integrate_difference() as well as an error estimator. We
+record timings for each section and run the code twice. In the first run, we
+hand a MappingQGeneric object to each stage of the program separately, where
+points get re-computed over and over again. In the second run, we use
+MappingQCache instead.
projects / dealii.git / commitdiff
