<a name="Intro"></a>
<h1>Introduction</h1>
-<p>
+
This is a rather short example which only shows some aspects of using
-higher order mappings. By <it>mapping</it> we mean the transformation
+higher order mappings. By <em>mapping</em> 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
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 <em>Q<sub>2</sub></em> approximation. If we
-use piecewise graphs of cubic polynomials, then this is a <em>Q<sub>3</sub></em>
+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.
-</p>
-<p>
+
+
For some differential equations, it is known that piecewise linear
-approximations of the boundary, i.e. Q<sub>1</sub> mappings, are not
+approximations of the boundary, i.e. $Q_1$ mappings, are not
sufficient if the boundary of the domain is curved. Examples are the
-biharmonic equation using C<sup>1</sup> elements, or the Euler
+biharmonic equation using $C^1$ elements, or the Euler
equation on domains with curved reflective boundaries. In these cases,
it is necessary to compute the integrals using a higher order
mapping. The reason, of course, is that if we do not use a higher
the order of convergence of the entire numerical scheme, irrespective
of the order of convergence of the discretization in the interior of
the domain.
-</p>
-<p>
+
+
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... by two
+calculating the value of $\pi=3.141592653589793238462643\ldots$ by two
different methods.
-</p>
-<p>
+
+
The first method uses a triangulated approximation of the circle with
unit radius and integrates the unit function 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 is not exactly pi. However, it is known that as
-we refine the triangulation, a Q<sub>p</sub> mapping approximates the boundary
-with an order <it>h<sup>p+1</sup></it>, where <it>h</it> is the mesh
+we refine the triangulation, a $Q_p$ mapping approximates the boundary
+with an order $h^{p+1}$, where $h$ is the mesh
width. 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.
-</p>
-<p>
+
+
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 the radius is equal to one.
-</p>
+
<a name="Results"></a>
<h1>Results</h1>
-<p>
+
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
-<pre>
+@code
set data style lines
set size 0.721, 1
set nokey
plot [-1:1][-1:1] "ball0_mapping_q1.dat"
-</pre>
+@endcode
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.
-</p>
-<p>
+
+
The following table shows the triangulated computational domain for
Q1, Q2, and Q3 mappings, for the original coarse grid (left), and a
once uniformly refined grid (right). If your browser does not display
<table "width=80%" align="center">
<tr>
<td>
- <img src="step-10.data/ball_mapping_q1_ref0.jpg" alt="Q1 mapping, coarse grid" height="300">
+ @image html step-10.ball_mapping_q1_ref0.png
</td>
<td>
- <img src="step-10.data/ball_mapping_q1_ref1.jpg" alt="Q1 mapping, once refined grid" height="300">
+ @image html step-10.ball_mapping_q1_ref1.png
</td>
</tr>
<tr>
<td>
- <img src="step-10.data/ball_mapping_q2_ref0.jpg" alt="Q2 mapping, coarse grid" height="300">
+ @image html step-10.ball_mapping_q2_ref0.png
</td>
<td>
- <img src="step-10.data/ball_mapping_q2_ref1.jpg" alt="Q2 mapping, once refined grid" height="300">
+ @image html step-10.ball_mapping_q2_ref1.png
</td>
</tr>
<tr>
<td>
- <img src="step-10.data/ball_mapping_q3_ref0.jpg" alt="Q3 mapping, coarse grid" height="300">
+ @image html step-10.ball_mapping_q3_ref0.png
</td>
<td>
- <img src="step-10.data/ball_mapping_q3_ref1.jpg" alt="Q3 mapping, once refined grid" height="300">
+ @image html step-10.ball_mapping_q3_ref1.png
</td>
</tr>
</table>
<table "width=80%" align="center">
<tr>
<td>
- <img src="step-10.data/quarter-q1.jpg" alt="Q1 mapping, coarse grid" height="300">
+ @image html step-10.quarter-q1.png
</td>
<td>
- <img src="step-10.data/quarter-q2.jpg" alt="Q2 mapping, coarse grid" height="300">
+ @image html step-10.quarter-q2.png
</td>
<td>
- <img src="step-10.data/quarter-q3.jpg" alt="Q3 mapping, coarse grid" height="300">
+ @image html step-10.quarter-q3.png
</td>
</tr>
</table>
usually not noticably slower than lower order ones, because the
additional computations are only performed on a small subset of all
cells.
-</p>
-<p>
+
+
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:
-<pre>
+@code
Computation of Pi by the area:
==============================
Degree = 1
320 3.1415926540544636 4.6467e-10 6.00
1280 3.1415926535970535 7.2602e-12 6.00
5120 3.1415926535899010 1.0805e-13 6.07
-</pre>
-</p>
+@endcode
+
+
-<p>
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<sub>4</sub> mapping, the last
+$\pi=3.141592653589793238462643$. Note that for the $Q_4$ mapping, the last
number is correct to 13 digits in both computations, which is already
-quite a lot. However, also note that for the Q<sub>1</sub> mapping, even on the
+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<sub>4</sub> mapping!
-</p>
+grid for a $Q_4$ mapping!
+
+
-<p>
The last column of the output shows the convergence order, in powers
-of the mesh width <it>h</it>. In the introduction, we had stated that
-the convergence order for a <it>Q<sub>p</sub></it> mapping should be
-<it>h<sup>p+1</sup></it>. However, in the example shown, the Q<sub>2</sub> and Q<sub>4</sub>
-mappings show a convergence order of <it>h<sup>p+2</sup></it>! This at
+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 $Q_2$ and $Q_4$
+mappings show a convergence order of $h^{p+2}$! This at
first surprising fact is readily explained by the particular boundary
we have chosen in this example. In fact, the circle is described by the function
-<it>sqrt(1-x<sup>2</sup>)</it>, which has the series expansion
-<it>1-x<sup>2</sup>/2-x<sup>4</sup>/8-x<sup>6</sup>/16+...</it>
-around <it>x=0</it>. Thus, for the quadratic mapping where the
+$\sqrt{1-x^2}$, which has the series expansion
+$1-x^2/2-x^4/8-x^6/16+\ldots$
+around $x=0$. Thus, for the quadratic mapping where the
truncation error of the quadratic approximation should be cubic, there
is no such term but only a quartic one, which raises the convergence
-order to 4, instead of 3. The same happens for the Q<sub>4</sub> mapping.
-</p>
+order to 4, instead of 3. The same happens for the $Q_4$ mapping.
+
<a name="Results"></a>
<h1>Results</h1>
-<p>
+
The output of this program consist of the console output, the eps
files including the grids, and the solutions given in gnuplot format.
-<code>
-<pre>
+@code
Cycle 0:
Number of active cells: 64
Number of degrees of freedom: 256
Time of assemble_system1: 0.05
Time of assemble_system2: 0.04
solution1 and solution2 coincide.
-Writing grid to <grid-0.eps>...
-Writing solution to <sol-0.gnuplot>...
+Writing grid to <grid-0.eps>...
+Writing solution to <sol-0.gnuplot>...
Cycle 1:
Number of active cells: 112
Time of assemble_system1: 0.09
Time of assemble_system2: 0.07
solution1 and solution2 coincide.
-Writing grid to <grid-1.eps>...
-Writing solution to <sol-1.gnuplot>...
+Writing grid to <grid-1.eps>...
+Writing solution to <sol-1.gnuplot>...
Cycle 2:
Number of active cells: 214
Time of assemble_system1: 0.17
Time of assemble_system2: 0.14
solution1 and solution2 coincide.
-Writing grid to <grid-2.eps>...
-Writing solution to <sol-2.gnuplot>...
+Writing grid to <grid-2.eps>...
+Writing solution to <sol-2.gnuplot>...
Cycle 3:
Number of active cells: 415
Time of assemble_system1: 0.32
Time of assemble_system2: 0.28
solution1 and solution2 coincide.
-Writing grid to <grid-3.eps>...
-Writing solution to <sol-3.gnuplot>...
+Writing grid to <grid-3.eps>...
+Writing solution to <sol-3.gnuplot>...
Cycle 4:
Number of active cells: 796
Time of assemble_system1: 0.62
Time of assemble_system2: 0.52
solution1 and solution2 coincide.
-Writing grid to <grid-4.eps>...
-Writing solution to <sol-4.gnuplot>...
+Writing grid to <grid-4.eps>...
+Writing solution to <sol-4.gnuplot>...
Cycle 5:
Number of active cells: 1561
Time of assemble_system1: 1.23
Time of assemble_system2: 1.03
solution1 and solution2 coincide.
-Writing grid to <grid-5.eps>...
-Writing solution to <sol-5.gnuplot>...
+Writing grid to <grid-5.eps>...
+Writing solution to <sol-5.gnuplot>...
+@endcode
-</pre>
-</code>
We see that, as expected, on each refinement step the two solutions
coincide. The difference measured in time of treating each face only
once (second version of the DG method) in comparison with treating
each face twice within a nested loop over all cells and all faces of
each cell (first version), is much less than one might have
expected. The gain is less than 20% on the last few refinement steps.
-</p>
-<p> First we show the solutions on the initial mesh, the mesh after two
+
+ First we show the solutions on the initial mesh, the mesh after two
and after five adaptive refinement steps.
-<p align="center">
- <img src="step-12.data/sol-0.gif" alt="sol-0.gnuplot" width="75%">
- <img src="step-12.data/sol-2.gif" alt="sol-2.gnuplot" width="75%">
- <img src="step-12.data/sol-5.gif" alt="sol-5.gnuplot" width="75%">
-</p>
-<p> Then we show the final grid (after 5 refinement steps). The
+@image step-12.sol-0.png
+@image step-12.sol-2.png
+@image step-12.sol-5.png
+
+
+Then we show the final grid (after 5 refinement steps). The
grid is largely concentrated in the vicinity of the jump of the
solution.
-<p align="center">
- <img src="step-12.data/grid-5.gif" alt="grid-5.eps" width="50%">
-</p>
-
-<p> And finally we show a plot of a 3d computation.
-<p align="center">
- <img src="step-12.data/sol-5-3d.gif" alt="sol-5-3d.eps" width="50%">
-</p>
+
+@image step-12.grid-5.png
+
+
+And finally we show a plot of a 3d computation.
+
+@image html step-12.sol-5-3d.png
+
<a name="Intro"></a>
<h1>Introduction</h1>
-<p>
+
This example shows the basic usage of the multilevel functions in
-<tt>deal.II</tt>. It solves the Helmholtz equation with Neumann boundary conditions
+deal.II. It solves the Helmholtz equation with Neumann boundary conditions
to avoid additional complications due to Dirichlet boundary conditions (there,
some library functions are missing). Therefore, the solution is the constant
function with value unity. In all other respects, it is similar to step 5.
-</p>
-<p> In order to allow sufficient flexibility in conjunction with systems of
+
+ In order to allow sufficient flexibility in conjunction with systems of
differential equations and block preconditioners, quite a few different objects
have to be created before starting the multilevel method. These are
<ul>
<li>MGMatrix, the matrix object having a special level multiplication, i.e. we
basically store one matrix per grid level and allow multiplication with it.
</ul>
-</p>
-<p>
+
+
These objects are combined in an object of type Multigrid, containing the
implementation of the V-cycle, which is in turn used by the preconditioner
PreconditionMG, ready for plug-in into a linear solver of the LAC library.
-</p>
-<p>
-The multilevel method in <tt>deal.II</tt> follows in many respects the outlines
+
+
+The multilevel method in deal.II follows in many respects the outlines
of the various publications by James Bramble, Joseph Pasciak and Jinchao Xu. In
order to understand many of the options, a rough familiarity with their work is
quite helpful.
-</p>
+
<a name="Results"></a>
<h1>Results</h1>
-<p>
+
There is not much to be said about the results of this program, apart
from that they look nice. All images were made using GMV from the
output files that the program wrote to disk. The first picture shows
the displacement as a vector field, where one vector is shown at each
vertex of the grid:
-</p>
-<p ALIGN=CENTER>
-<IMG SRC="step-8.data/displacement-vectors.jpg" ALT="displacement-vectors" WIDTH="300">
-</p>
-<p>
+@image html step-8.vectors.png
+
+
You can clearly see the sources of x-displacement around x=0.5 and
x=-0.5, and of y-displacement at the origin. The next image shows the
final grid after eight steps of refinement:
-</p>
-<p ALIGN=CENTER>
-<IMG SRC="step-8.data/final-grid.jpg" ALT="final-grid" WIDTH="300">
-</p>
-<p>
+@image html step-8.grid.png
+
+
+
Finally, the x-displacement and y-displacement are displayed separately:
-</p>
-<p>
+
+
<TABLE WIDTH="100%">
<tr>
<td>
-<IMG SRC="step-8.data/displacement-x.jpg" ALT="displacement-x" WIDTH="300">
+@image html step-8.x.png
</td>
<td>
-<IMG SRC="step-8.data/displacement-y.jpg" ALT="displacement-y" WIDTH="300">
+@image html step-8.y.png
</td>
</tr>
</table>
-</p>
-<p>
+
+
It should be noted that intuitively one would have expected the
solution to be symmetric about the x- and y-axes since the x- and
y-forces are symmetric with respect to these axes. However, the force
considered as a vector is not symmetric and consequently neither is
the solution.
-</p>
+
<a name="Results"></a>
<h1>Results</h1>
-<p>
+
The results of this program are not particularly spectacular. They
consist of the console output, some grid files, and the solution on
the finest grid. First for the console output:
-<code>
-<pre>
+@code
Cycle 0:
Number of active cells: 256
Number of degrees of freedom: 289
Cycle 5:
Number of active cells: 26002
Number of degrees of freedom: 30035
-</pre>
-</code>
+@endcode
As can be seen, quite a number of cells is used on the finest level to
resolve the features of the solution. The final grid showing this is
displayed in the following picture:
-</p>
-<p align="center">
- <img src="step-9.data/final-grid.gif" alt="final grid" width="50%">
-</p>
-<p>
+@image html step-9.grid.png
+
+
+
The structure of the grid will be understandable by looking at the
solution itself:
-</p>
-<p align="center">
- <img src="step-9.data/final-solution.jpeg" alt="final solution" width="50%">
-</p>
-<p>
+@image html step-9.solution.png
+
+
+
Note that the solution is created by that part that is transported
along the wiggled advection field from the left and lower boundaries
to the top right, and the part that is created by the source in the
lower left corner, and the results of which are also transported
along. The grid shown above is well-adapted to resolve these
features.
-</p>
\ No newline at end of file