From: Martin Kronbichler Date: Tue, 4 Jun 2019 14:26:44 +0000 (+0200) Subject: Write introduction X-Git-Tag: v9.2.0-rc1~1428^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a287a40fce58c135f31b2d0378a0ba8823a467b3;p=dealii.git Write introduction --- diff --git a/doc/doxygen/images/circular_mesh_boundary_cells.png b/doc/doxygen/images/circular_mesh_boundary_cells.png new file mode 100644 index 0000000000..04b546986f Binary files /dev/null and b/doc/doxygen/images/circular_mesh_boundary_cells.png differ diff --git a/doc/doxygen/images/torus_cylindrical_inner_manifold.png b/doc/doxygen/images/torus_cylindrical_inner_manifold.png new file mode 100644 index 0000000000..9dc3f642b4 Binary files /dev/null and b/doc/doxygen/images/torus_cylindrical_inner_manifold.png differ diff --git a/doc/doxygen/images/torus_no_inner_manifold.png b/doc/doxygen/images/torus_no_inner_manifold.png new file mode 100644 index 0000000000..53763a7c84 Binary files /dev/null and b/doc/doxygen/images/torus_no_inner_manifold.png differ diff --git a/doc/doxygen/images/torus_transfinite_manifold.png b/doc/doxygen/images/torus_transfinite_manifold.png new file mode 100644 index 0000000000..06f8355665 Binary files /dev/null and b/doc/doxygen/images/torus_transfinite_manifold.png differ diff --git a/examples/step-65/doc/builds-on b/examples/step-65/doc/builds-on new file mode 100644 index 0000000000..76df1d0b00 --- /dev/null +++ b/examples/step-65/doc/builds-on @@ -0,0 +1 @@ +step-49 \ No newline at end of file diff --git a/examples/step-65/doc/intro.dox b/examples/step-65/doc/intro.dox index 26b84f25d5..542d86fb95 100644 --- a/examples/step-65/doc/intro.dox +++ b/examples/step-65/doc/intro.dox @@ -8,7 +8,384 @@ This program was contributed by Martin Kronbichler.

Introduction

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

Working with manifolds

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

The class TransfiniteInterpolationManifold

+ +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 +all 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 Gordon and Hall. 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. + +

Transfinite interpolation is expensive and how to deal with it

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

The test case

+ +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: + + + +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. diff --git a/examples/step-65/doc/kind b/examples/step-65/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-65/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-65/doc/results.dox b/examples/step-65/doc/results.dox index 9738b1c95c..3480230d25 100644 --- a/examples/step-65/doc/results.dox +++ b/examples/step-65/doc/results.dox @@ -1,6 +1,8 @@

Results

+The meshes created by this program are discussed in the introduction. +

Program output

@code @@ -11,7 +13,7 @@ Scanning dependencies of target step-65 [ 66%] Built target step-65 [100%] Run step-65 with Release configuration -====== Running with the basic MappingQGeneric class ====== +====== Running with the basic MappingQGeneric class ====== Number of active cells: 6656 Number of degrees of freedom: 181609 @@ -33,7 +35,7 @@ Scanning dependencies of target step-65 | Write output | 2 | 10.3s | 18% | +---------------------------------+-----------+------------+------------+ -====== Running with the optimized MappingQCache class ====== +====== Running with the optimized MappingQCache class ====== Memory consumption cache: 22.9976 MB Number of active cells: 6656 diff --git a/examples/step-65/doc/tooltip b/examples/step-65/doc/tooltip new file mode 100644 index 0000000000..aa4fc6a1c4 --- /dev/null +++ b/examples/step-65/doc/tooltip @@ -0,0 +1 @@ +Geometry: Working efficiently with expensive manifolds.