 |
Reference documentation for deal.II version 9.2.0
|
|
Reference documentation for deal.II version 9.2.0
|
This tutorial depends on step-7.
| Table of contents | |
|---|---|
This is a rather short example which only shows some aspects of using higher order mappings. By mapping we mean the transformation between the unit cell (i.e. the unit line, square, or cube) to the cells in real space. In all the previous examples, we have implicitly used linear or d-linear mappings; you will not have noticed this at all, since this is what happens if you do not do anything special. However, if your domain has curved boundaries, there are cases where the piecewise linear approximation of the boundary (i.e. by straight line segments) is not sufficient, and you want that your computational domain is an approximation to the real domain using curved boundaries as well. If the boundary approximation uses piecewise quadratic parabolas to approximate the true boundary, then we say that this is a quadratic or \(Q_2\) approximation. If we use piecewise graphs of cubic polynomials, then this is a \(Q_3\) approximation, and so on.
For some differential equations, it is known that piecewise linear approximations of the boundary, i.e. \(Q_1\) mappings, are not sufficient if the boundary of the exact domain is curved. Examples are the biharmonic equation using \(C^1\) elements, or the Euler equations of gas dynamics on domains with curved reflective boundaries. In these cases, it is necessary to compute the integrals using a higher order mapping. If we do not use such a higher order mapping, the order of approximation of the boundary dominates the order of convergence of the entire numerical scheme, irrespective of the order of convergence of the discretization in the interior of the domain.
Rather than demonstrating the use of higher order mappings with one of these more complicated examples, we do only a brief computation: calculating the value of \(\pi=3.141592653589793238462643\ldots\) by two different methods.
The first method uses a triangulated approximation of the circle with unit radius and integrates the function that is constant one over it. Of course, if the domain were the exact unit circle, then the area would be \(\pi\), but since we only use an approximation by piecewise polynomial segments, the value of the area we integrate over is not exactly \(\pi\). However, it is known that as we refine the triangulation, a \(Q_p\) mapping approximates the boundary with an order \(h^{p+1}\), where \(h\) is the mesh size. We will check the values of the computed area of the circle and their convergence towards \(\pi\) under mesh refinement for different mappings. We will also find a convergence behavior that is surprising at first, but has a good explanation.
The second method works similarly, but this time does not use the area of the triangulated unit circle, but rather its perimeter. \(\pi\) is then approximated by half of the perimeter, as we choose the radius equal to one.
The first of the following include files are probably well-known by now and need no further explanation.
This include file is new. Even if we are not solving a PDE in this tutorial, we want to use a dummy finite element with zero degrees of freedoms provided by the FE_Nothing class.
The following header file is also new: in it, we declare the MappingQ class which we will use for polynomial mappings of arbitrary order:
And this again is C++:
The last step is as in previous programs:
Now, as we want to compute the value of \(\pi\), we have to compare to something. These are the first few digits of \(\pi\), which we define beforehand for later use. Since we would like to compute the difference between two numbers which are quite accurate, with the accuracy of the computed approximation to \(\pi\) being in the range of the number of digits which a double variable can hold, we rather declare the reference value as a long double and give it a number of extra digits:
Then, the first task will be to generate some output. Since this program is so small, we do not employ object oriented techniques in it and do not declare classes (although, of course, we use the object oriented features of the library). Rather, we just pack the functionality into separate functions. We make these functions templates on the number of space dimensions to conform to usual practice when using deal.II, although we will only use them for two space dimensions.
The first of these functions just generates a triangulation of a circle (hyperball) and outputs the \(Q_p\) mapping of its cells for different values of p. Then, we refine the grid once and do so again.
So first generate a coarse triangulation of the circle and associate a suitable boundary description to it. By default, GridGenerator::hyper_ball attaches a SphericalManifold to the boundary (and uses FlatManifold for the interior) so we simply call that function and move on:
Then alternate between generating output on the current mesh for \(Q_1\), \(Q_2\), and \(Q_3\) mappings, and (at the end of the loop body) refining the mesh once globally.
For this, first set up an object describing the mapping. This is done using the MappingQ class, which takes as argument to the constructor the polynomial degree which it shall use.
As a side note, for a piecewise linear mapping, you could give a value of 1 to the constructor of MappingQ, but there is also a class MappingQ1 that achieves the same effect. Historically, it did a lot of things in a simpler way than MappingQ but is today just a wrapper around the latter. It is, however, still the class that is used implicitly in many places of the library if you do not specify another mapping explicitly.
In order to actually write out the present grid with this mapping, we set up an object which we will use for output. We will generate Gnuplot output, which consists of a set of lines describing the mapped triangulation. By default, only one line is drawn for each face of the triangulation, but since we want to explicitly see the effect of the mapping, we want to have the faces in more detail. This can be done by passing the output object a structure which contains some flags. In the present case, since Gnuplot can only draw straight lines, we output a number of additional points on the faces so that each face is drawn by 30 small lines instead of only one. This is sufficient to give us the impression of seeing a curved line, rather than a set of straight lines.
Finally, generate a filename and a file for output:
Then write out the triangulation to this file. The last argument of the function is a pointer to a mapping object. This argument has a default value, and if no value is given a simple MappingQ1 object is taken, which we briefly described above. This would then result in a piecewise linear approximation of the true boundary in the output.
At the end of the loop, refine the mesh globally.
Now we proceed with the main part of the code, the approximation of \(\pi\). The area of a circle is of course given by \(\pi r^2\), so having a circle of radius 1, the area represents just the number that is searched for. The numerical computation of the area is performed by integrating the constant function of value 1 over the whole computational domain, i.e. by computing the areas \(\int_K 1 dx=\int_{\hat K} 1 \ \textrm{det}\ J(\hat x) d\hat x \approx \sum_i \textrm{det} \ J(\hat x_i)w(\hat x_i)\), where the sum extends over all quadrature points on all active cells in the triangulation, with \(w(x_i)\) being the weight of quadrature point \(x_i\). The integrals on each cell are approximated by numerical quadrature, hence the only additional ingredient we need is to set up a FEValues object that provides the corresponding ‘JxW’ values of each cell. (Note that ‘JxW’ is meant to abbreviate Jacobian determinant times weight; since in numerical quadrature the two factors always occur at the same places, we only offer the combined quantity, rather than two separate ones.) We note that here we won't use the FEValues object in its original purpose, i.e. for the computation of values of basis functions of a specific finite element at certain quadrature points. Rather, we use it only to gain the ‘JxW’ at the quadrature points, irrespective of the (dummy) finite element we will give to the constructor of the FEValues object. The actual finite element given to the FEValues object is not used at all, so we could give any.
For the numerical quadrature on all cells we employ a quadrature rule of sufficiently high degree. We choose QGauss that is of order 8 (4 points), to be sure that the errors due to numerical quadrature are of higher order than the order (maximal 6) that will occur due to the order of the approximation of the boundary, i.e. the order of the mappings employed. Note that the integrand, the Jacobian determinant, is not a polynomial function (rather, it is a rational one), so we do not use Gauss quadrature in order to get the exact value of the integral as done often in finite element computations, but could as well have used any quadrature formula of like order instead.
Now start by looping over polynomial mapping degrees=1..4:
First generate the triangulation, the boundary and the mapping object as already seen.
We now create a finite element. Unlike the rest of the example programs, we do not actually need to do any computations with shape functions; we only need the ‘JxW’ values from an FEValues object. Hence we use the special finite element class FE_Nothing which has exactly zero degrees of freedom per cell (as the name implies, the local basis on each cell is the empty set). A more typical usage of FE_Nothing is shown in step-46.
Likewise, we need to create a DoFHandler object. We do not actually use it, but it will provide us with ‘active_cell_iterators’ that are needed to reinitialize the FEValues object on each cell of the triangulation.
Now we set up the FEValues object, giving the Mapping, the dummy finite element and the quadrature object to the constructor, together with the update flags asking for the ‘JxW’ values at the quadrature points only. This tells the FEValues object that it needs not compute other quantities upon calling the reinit function, thus saving computation time.
The most important difference in the construction of the FEValues object compared to previous example programs is that we pass a mapping object as first argument, which is to be used in the computation of the mapping from unit to real cell. In previous examples, this argument was omitted, resulting in the implicit use of an object of type MappingQ1.
We employ an object of the ConvergenceTable class to store all important data like the approximated values for \(\pi\) and the error with respect to the true value of \(\pi\). We will also use functions provided by the ConvergenceTable class to compute convergence rates of the approximations to \(\pi\).
Now we loop over several refinement steps of the triangulation.
In this loop we first add the number of active cells of the current triangulation to the table. This function automatically creates a table column with superscription ‘cells’, in case this column was not created before.
Then we distribute the degrees of freedom for the dummy finite element. Strictly speaking we do not need this function call in our special case but we call it to make the DoFHandler happy – otherwise it would throw an assertion in the FEValues::reinit function below.
We define the variable area as ‘long double’ like we did for the pi variable before.
Now we loop over all cells, reinitialize the FEValues object for each cell, and add up all the ‘JxW’ values for this cell to ‘area’...
...and store the resulting area values and the errors in the table. We need a static cast to double as there is no add_value(string, long double) function implemented. Note that this also concerns the second call as the fabs function in the std namespace is overloaded on its argument types, so there exists a version taking and returning a long double, in contrast to the global namespace where only one such function is declared (which takes and returns a double).
We want to compute the convergence rates of the ‘error’ column. Therefore we need to omit the other columns from the convergence rate evaluation before calling ‘evaluate_all_convergence_rates’
Finally we set the precision and scientific mode for output of some of the quantities...
...and write the whole table to std::cout.
The following, second function also computes an approximation of \(\pi\) but this time via the perimeter \(2\pi r\) of the domain instead of the area. This function is only a variation of the previous function. So we will mainly give documentation for the differences.
We take the same order of quadrature but this time a ‘dim-1’ dimensional quadrature as we will integrate over (boundary) lines rather than over cells.
We loop over all degrees, create the triangulation, the boundary, the mapping, the dummy finite element and the DoFHandler object as seen before.
Then we create a FEFaceValues object instead of a FEValues object as in the previous function. Again, we pass a mapping as first argument.
Now we run over all cells and over all faces of each cell. Only the contributions of the ‘JxW’ values on boundary faces are added to the long double variable ‘perimeter’.
We reinit the FEFaceValues object with the cell iterator and the number of the face.
Then store the evaluated values in the table...
...and end this function as we did in the previous one:
The following main function just calls the above functions in the order of their appearance. Apart from this, it looks just like the main functions of previous tutorial programs.
The program performs two tasks, the first being to generate a visualization of the mapped domain, the second to compute pi by the two methods described. Let us first take a look at the generated graphics. They are generated in Gnuplot format, and can be viewed with the commands
or using one of the other filenames. The second line makes sure that the aspect ratio of the generated output is actually 1:1, i.e. a circle is drawn as a circle on your screen, rather than as an ellipse. The third line switches off the key in the graphic, as that will only print information (the filename) which is not that important right now. Similarly, the third and fourth disable tic marks. The plot is then generated with a specific line width ("`lw`", here set to 4) and line type ("`lt`", here chosen by saying that the line should be drawn using the RGB color "black").
The following table shows the triangulated computational domain for \(Q_1\), \(Q_2\), and \(Q_3\) mappings, for the original coarse grid (left), and a once uniformly refined grid (right).
These pictures show the obvious advantage of higher order mappings: they approximate the true boundary quite well also on rather coarse meshes. To demonstrate this a little further, here is part of the upper right quarter circle of the coarse meshes with \(Q_2\) and \(Q_3\) mappings, where the dashed red line marks the actual circle:
Obviously the quadratic mapping approximates the boundary quite well, while for the cubic mapping the difference between approximated domain and true one is hardly visible already for the coarse grid. You can also see that the mapping only changes something at the outer boundaries of the triangulation. In the interior, all lines are still represented by linear functions, resulting in additional computations only on cells at the boundary. Higher order mappings are therefore usually not noticeably slower than lower order ones, because the additional computations are only performed on a small subset of all cells.
The second purpose of the program was to compute the value of pi to good accuracy. This is the output of this part of the program:
One of the immediate observations from the output is that in all cases the values converge quickly to the true value of \(\pi=3.141592653589793238462643\). Note that for the \(Q_4\) mapping, we are already in the regime of roundoff errors and the convergence rate levels off, which is already quite a lot. However, also note that for the \(Q_1\) mapping, even on the finest grid the accuracy is significantly worse than on the coarse grid for a \(Q_3\) mapping!
The last column of the output shows the convergence order, in powers of the mesh width \(h\). In the introduction, we had stated that the convergence order for a \(Q_p\) mapping should be \(h^{p+1}\). However, in the example shown, the order is rather \(h^{2p}\)! This at first surprising fact is explained by the properties of the \(Q_p\) mapping. At order p, it uses support points that are based on the p+1 point Gauss-Lobatto quadrature rule that selects the support points in such a way that the quadrature rule converges at order 2p. Even though these points are here only used for interpolation of a pth order polynomial, we get a superconvergence effect when numerically evaluating the integral, resulting in the observed high order of convergence. (This effect is also discussed in detail in the following publication: A. Bonito, A. Demlow, and J. Owen: "A priori error estimates for finite element approximations to eigenvalues and eigenfunctions of the Laplace-Beltrami operator", submitted, 2018.)
1.8.17