Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
The step-10 tutorial program

This tutorial depends on step-7.

Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program

Introduction

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 a unit magnitude constant function ( \(f = 1\)) 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.

Note
This tutorial shows in essence how to choose a particular mapping for integrals, by attaching a particular geometry to the triangulation (as had already been done in step-1, for example) and then passing a mapping argument to the FEValues class that is used for all integrals in deal.II. The geometry we choose is a circle, for which deal.II already has a class (SphericalManifold) that can be used. If you want to define your own geometry, for example because it is complicated and cannot be described by the classes already available in deal.II, you will want to read through step-53.

The commented program

The first of the following include files are probably well-known by now and need no further explanation.

  #include <deal.II/base/quadrature_lib.h>
  #include <deal.II/base/convergence_table.h>
  #include <deal.II/grid/grid_generator.h>
  #include <deal.II/grid/manifold_lib.h>
  #include <deal.II/grid/tria.h>
  #include <deal.II/grid/grid_out.h>
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/fe/fe_values.h>
 

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.

  #include <deal.II/fe/fe_nothing.h>
 

The following header file is also new: in it, we declare the MappingQ class which we will use for polynomial mappings of arbitrary order:

  #include <deal.II/fe/mapping_q.h>
 

And this again is C++:

  #include <iostream>
  #include <fstream>
  #include <cmath>
 

The last step is as in previous programs:

  namespace Step10
  {
  using namespace dealii;
 

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 and throw an exception when attempted to use for any other spatial dimension.

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.

  template <int dim>
  void gnuplot_output()
  {
  std::cout << "Output of grids into gnuplot files:" << std::endl
  << "===================================" << std::endl;
 

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:

 
void hyper_ball(Triangulation< dim > &tria, const Point< dim > &center=Point< dim >(), const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

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 (unsigned int refinement = 0; refinement < 2; ++refinement)
  {
  std::cout << "Refinement level: " << refinement << std::endl;
 
  std::string filename_base = "ball_" + std::to_string(refinement);
 
  for (unsigned int degree = 1; degree < 4; ++degree)
  {
  std::cout << "Degree = " << degree << std::endl;
 

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.

  const MappingQ<dim> mapping(degree);

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.

  GridOut grid_out;
  GridOutFlags::Gnuplot gnuplot_flags(false, 60);
  grid_out.set_flags(gnuplot_flags);
 

Finally, generate a filename and a file for output:

  std::string filename =
  filename_base + "_mapping_q_" + std::to_string(degree) + ".dat";
  std::ofstream gnuplot_file(filename);
 

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.

  grid_out.write_gnuplot(triangulation, gnuplot_file, &mapping);
  }
  std::cout << std::endl;
 

At the end of the loop, refine the mesh globally.

  }
  }
 
void refine_global(const unsigned int times=1)

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.

  template <int dim>
  void compute_pi_by_area()
  {
  std::cout << "Computation of Pi by the area:" << std::endl
  << "==============================" << std::endl;
 

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.

  const QGauss<dim> quadrature(4);
 

Now start by looping over polynomial mapping degrees=1..4:

  for (unsigned int degree = 1; degree < 5; ++degree)
  {
  std::cout << "Degree = " << degree << std::endl;
 

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.

  const FE_Nothing<dim> fe;
 

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.

  FEValues<dim> fe_values(mapping, fe, quadrature, update_JxW_values);
 
@ update_JxW_values
Transformed quadrature weights.

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.

  for (unsigned int refinement = 0; refinement < 6;
  ++refinement, triangulation.refine_global(1))
  {

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.

  table.add_value("cells", triangulation.n_active_cells());
 
unsigned int n_active_cells() const

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.

  dof_handler.distribute_dofs(fe);
 

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

  double area = 0;
  for (const auto &cell : dof_handler.active_cell_iterators())
  {
  fe_values.reinit(cell);
  for (unsigned int i = 0; i < fe_values.n_quadrature_points; ++i)
  area += fe_values.JxW(i);
  }
 

...and store the resulting area values and the errors in the table:

  table.add_value("eval.pi", area);
  table.add_value("error", std::fabs(area - numbers::PI));
  }
 
static constexpr double PI
Definition numbers.h:259

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

  table.omit_column_from_convergence_rate_evaluation("cells");
  table.omit_column_from_convergence_rate_evaluation("eval.pi");
  table.evaluate_all_convergence_rates(
 

Finally we set the precision and scientific mode for output of some of the quantities...

  table.set_precision("eval.pi", 16);
  table.set_scientific("error", true);
 

...and write the whole table to std::cout.

  table.write_text(std::cout);
 
  std::cout << std::endl;
  }
  }
 
 

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.

  template <int dim>
  void compute_pi_by_perimeter()
  {
  std::cout << "Computation of Pi by the perimeter:" << std::endl
  << "===================================" << std::endl;
 

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.

  const QGauss<dim - 1> quadrature(4);
 

We loop over all degrees, create the triangulation, the boundary, the mapping, the dummy finite element and the DoFHandler object as seen before.

  for (unsigned int degree = 1; degree < 5; ++degree)
  {
  std::cout << "Degree = " << degree << std::endl;
 
  const MappingQ<dim> mapping(degree);
  const FE_Nothing<dim> fe;
 
 

Then we create a FEFaceValues object instead of a FEValues object as in the previous function. Again, we pass a mapping as first argument.

  FEFaceValues<dim> fe_face_values(mapping,
  fe,
  quadrature,
 
  for (unsigned int refinement = 0; refinement < 6;
  ++refinement, triangulation.refine_global(1))
  {
  table.add_value("cells", triangulation.n_active_cells());
 
  dof_handler.distribute_dofs(fe);
 
void add_value(const std::string &key, const T value)

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 variable perimeter.

  double perimeter = 0;
  for (const auto &cell : dof_handler.active_cell_iterators())
  for (const auto &face : cell->face_iterators())
  if (face->at_boundary())
  {

We reinit the FEFaceValues object with the cell iterator and the number of the face.

  fe_face_values.reinit(cell, face);
  for (unsigned int i = 0;
  i < fe_face_values.n_quadrature_points;
  ++i)
  perimeter += fe_face_values.JxW(i);
  }

Then store the evaluated values in the table...

  table.add_value("eval.pi", static_cast<double>(perimeter / 2.0));
  table.add_value("error", std::fabs(perimeter / 2.0 - numbers::PI));
  }
 

...and end this function as we did in the previous one:

  table.omit_column_from_convergence_rate_evaluation("cells");
  table.omit_column_from_convergence_rate_evaluation("eval.pi");
  table.evaluate_all_convergence_rates(
 
  table.set_precision("eval.pi", 16);
  table.set_scientific("error", true);
 
  table.write_text(std::cout);
 
  std::cout << std::endl;
  }
  }
  } // namespace Step10
 
 

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.

  int main()
  {
  try
  {
  std::cout.precision(16);
 
  const unsigned int dim = 2;
 
  Step10::gnuplot_output<dim>();
 
  Step10::compute_pi_by_area<dim>();
  Step10::compute_pi_by_perimeter<dim>();
  }
  catch (std::exception &exc)
  {
  std::cerr << std::endl
  << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Exception on processing: " << std::endl
  << exc.what() << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
 
  return 1;
  }
  catch (...)
  {
  std::cerr << std::endl
  << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Unknown exception!" << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
 
  return 0;
  }

Results

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

set style data lines
set size ratio -1
unset key
unset xtics
unset ytics
plot [-1:1][-1:1] "ball_0_mapping_q_1.dat" lw 4 lt rgb "black"

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 fourth and fifth disable tick 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).

Five-cell discretization of the disk.
20-cell discretization of the disk (i.e., five cells
              refined once).
Five-cell discretization of the disk with quadratic edges. The
              boundary is nearly indistinguishable from the actual circle.
20-cell discretization with quadratic edges.
Five-cell discretization of the disk with cubic edges. The
              boundary is nearly indistinguishable from the actual circle.
20-cell discretization with cubic edges.

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:

Close-up of quadratic discretization. The distance between the
         quadratic interpolant and the actual circle is small.
Close-up of cubic discretization. The distance between the
         cubic interpolant and the actual circle is very small.

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 the last time it was updated:

Output of grids into gnuplot files:
===================================
Refinement level: 0
Degree = 1
Degree = 2
Degree = 3
Refinement level: 1
Degree = 1
Degree = 2
Degree = 3
Computation of Pi by the area:
==============================
Degree = 1
cells eval.pi error
5 1.9999999999999984 1.1416e+00 -
20 2.8284271247461947 3.1317e-01 1.87
80 3.0614674589207289 8.0125e-02 1.97
320 3.1214451522579965 2.0148e-02 1.99
1280 3.1365484905461294 5.0442e-03 2.00
5120 3.1403311569547410 1.2615e-03 2.00
Degree = 2
cells eval.pi error
5 3.1045694996615865 3.7023e-02 -
20 3.1391475703122351 2.4451e-03 3.92
80 3.1414377167038401 1.5494e-04 3.98
320 3.1415829366418513 9.7169e-06 4.00
1280 3.1415920457578785 6.0783e-07 4.00
5120 3.1415926155920988 3.7998e-08 4.00
Degree = 3
cells eval.pi error
5 3.1410033851499315 5.8927e-04 -
20 3.1415830393583946 9.6142e-06 5.94
80 3.1415925017363939 1.5185e-07 5.98
320 3.1415926512106185 2.3792e-09 6.00
1280 3.1415926535527783 3.7015e-11 6.01
5120 3.1415926535891936 5.9952e-13 5.95
Degree = 4
cells eval.pi error
5 3.1415871927401144 5.4608e-06 -
20 3.1415926314742491 2.2116e-08 7.95
80 3.1415926535026268 8.7166e-11 7.99
320 3.1415926535894005 3.9257e-13 7.79
1280 3.1415926535899774 1.8430e-13 1.09
5120 3.1415926535897669 2.6201e-14 2.81
Computation of Pi by the perimeter:
===================================
Degree = 1
cells eval.pi error
5 2.8284271247461903 3.1317e-01 -
20 3.0614674589207187 8.0125e-02 1.97
80 3.1214451522580493 2.0148e-02 1.99
320 3.1365484905459371 5.0442e-03 2.00
1280 3.1403311569547530 1.2615e-03 2.00
5120 3.1412772509327755 3.1540e-04 2.00
Degree = 2
cells eval.pi error
5 3.1248930668550590 1.6700e-02 -
20 3.1404050605605454 1.1876e-03 3.81
80 3.1415157631807005 7.6890e-05 3.95
320 3.1415878042798600 4.8493e-06 3.99
1280 3.1415923498174543 3.0377e-07 4.00
5120 3.1415926345931982 1.8997e-08 4.00
Degree = 3
cells eval.pi error
5 3.1414940401456053 9.8613e-05 -
20 3.1415913432549156 1.3103e-06 6.23
80 3.1415926341726905 1.9417e-08 6.08
320 3.1415926532906924 2.9910e-10 6.02
1280 3.1415926535851346 4.6585e-12 6.00
5120 3.1415926535897158 7.7272e-14 5.91
Degree = 4
cells eval.pi error
5 3.1415921029432572 5.5065e-07 -
20 3.1415926513737582 2.2160e-09 7.96
80 3.1415926535810699 8.7232e-12 7.99
320 3.1415926535897576 3.5527e-14 7.94
1280 3.1415926535897896 3.5527e-15 3.32
5120 3.1415926535897940 8.8818e-16 2.00
unsigned int level
Definition grid_out.cc:4626
Note
Once the error reaches a level on the order of \(10^{-13}\) to \(10^{-15}\), it is essentially dominated by round-off and consequently dominated by what precisely the library is doing in internal computations. Since these things change, the precise values and errors change from release to release at these round-off levels, though the overall order of errors should of course remain the same. See also the comment below in the section on Possibilities for extensions about how to compute these results more accurately.

One of the immediate observations from the output above 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.)

Possibilities for extensions

As the table of numbers copied from the output of the program shows above, it is not very difficult to compute the value of \(\pi\) to 13 or 15 digits. But, the output also shows that once we approach the level of accuracy with which double precision numbers store information (namely, with roughly 16 digits of accuracy), we no longer see the expected convergence order and the error no longer decreases with mesh refinement as anticipated. This is because both within this code and within the many computations that happen within deal.II itself, each operation incurs an error on the order of \(10^{-16}\); adding such errors many times over then results in an error that may be on the order of \(10^{-14}\), which will dominate the discretization error after a number of refinement steps and consequently destroy the convergence rate.

The question is whether one can do anything about this. One thought is to use a higher-precision data type. For example, one could think of declaring both the area and perimeter variables in compute_pi_by_area() and compute_pi_by_perimeter() with data type long double. long double is a data type that is not well specified in the C++ standard but at least on Intel processors has around 19, instead of around 16, digits of accuracy. If we were to do that, we would get results that differ from the ones shown above. However, maybe counter-intuitively, they are not uniformly better. For example, when computing \(\pi\) by the area, at the time of writing these sentences we get these values with double precision for degree 4:

5 3.1415871927401144 5.4608e-06 -
20 3.1415926314742491 2.2116e-08 7.95
80 3.1415926535026268 8.7166e-11 7.99
320 3.1415926535894005 3.9257e-13 7.79
1280 3.1415926535899774 1.8430e-13 1.09
5120 3.1415926535897669 2.6201e-14 2.81

On the other hand, the results are as follows when using long double:

cells eval.pi error
5 3.1415871927401136 5.4608e-06 -
20 3.1415926314742446 2.2116e-08 7.95
80 3.1415926535026215 8.7172e-11 7.99
320 3.1415926535894516 3.4157e-13 8.00
1280 3.1415926535897918 1.5339e-15 7.80
5120 3.1415926535897927 5.2649e-16 1.54

Indeed, here we get results that are approximately 50 times as accurate. On the other hand, when computing \(\pi\) by the perimeter, we get this with double precision:

5 3.1415921029432572 5.5065e-07 -
20 3.1415926513737582 2.2160e-09 7.96
80 3.1415926535810699 8.7232e-12 7.99
320 3.1415926535897576 3.5527e-14 7.94
1280 3.1415926535897896 3.5527e-15 3.32
5120 3.1415926535897940 8.8818e-16 2.00

Whereas we get the following with long double:

5 3.1415921029432572 5.5065e-07 -
20 3.1415926513737595 2.2160e-09 7.96
80 3.1415926535810703 8.7230e-12 7.99
320 3.1415926535897576 3.5705e-14 7.93
1280 3.1415926535897918 1.3785e-15 4.70
5120 3.1415926535897944 1.3798e-15 -0.00

Here, using double precision is more accurate by about a factor of two. (Of course, in all cases, we have computed \(\pi\) with more accuracy than any engineer would ever want to know.)

What explains this unpredictability? In general, round-off errors can be thought of as random, and add up in ways that are not worth thinking too much about; we should therefore always treat any accuracy beyond, say, thirteen digits as suspect. Thus, it is probably not worth spending too much time on wondering why we get different winners and losers in the data type exchange from double and long double. The accuracy of the results is also largely not determined by the precision of the data type in which we accumulate each cell's (or face's) contributions, but the accuracy of what deal.II gives us via FEValues::JxW() and FEFaceValues::JxW(), which always uses double precision and which we cannot directly affect.

But there are cases where one can do something about the precision, and it is worth at least mentioning the name of the most well-known algorithm in this area. Specifically, what we are doing when we add contributions into the area and perimeter values is that we are adding together positive numbers as we do here. In general, the round-off errors associated with each of these numbers is random, and if we add up contributions of substantially different sizes, then we will likely be dominated by the error in the largest contributions. One can avoid this by adding up numbers sorted by their size, and this may then result in marginally more accurate end results. The algorithm that implements this is typically called Kahan's summation algorithm. While one could play with it in the current context, it is likely not going to improve the accuracy in ways that will truly matter.

The plain program

/* ------------------------------------------------------------------------
*
* SPDX-License-Identifier: LGPL-2.1-or-later
* Copyright (C) 2001 - 2023 by the deal.II authors
*
* This file is part of the deal.II library.
*
* Part of the source code is dual licensed under Apache-2.0 WITH
* LLVM-exception OR LGPL-2.1-or-later. Detailed license information
* governing the source code and code contributions can be found in
* LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
*
* ------------------------------------------------------------------------
*
* Authors: Wolfgang Bangerth, Ralf Hartmann, University of Heidelberg, 2001
*/
#include <iostream>
#include <fstream>
#include <cmath>
namespace Step10
{
using namespace dealii;
template <int dim>
void gnuplot_output()
{
std::cout << "Output of grids into gnuplot files:" << std::endl
<< "===================================" << std::endl;
for (unsigned int refinement = 0; refinement < 2; ++refinement)
{
std::cout << "Refinement level: " << refinement << std::endl;
std::string filename_base = "ball_" + std::to_string(refinement);
for (unsigned int degree = 1; degree < 4; ++degree)
{
std::cout << "Degree = " << degree << std::endl;
const MappingQ<dim> mapping(degree);
GridOut grid_out;
GridOutFlags::Gnuplot gnuplot_flags(false, 60);
grid_out.set_flags(gnuplot_flags);
std::string filename =
filename_base + "_mapping_q_" + std::to_string(degree) + ".dat";
std::ofstream gnuplot_file(filename);
grid_out.write_gnuplot(triangulation, gnuplot_file, &mapping);
}
std::cout << std::endl;
}
}
template <int dim>
void compute_pi_by_area()
{
std::cout << "Computation of Pi by the area:" << std::endl
<< "==============================" << std::endl;
const QGauss<dim> quadrature(4);
for (unsigned int degree = 1; degree < 5; ++degree)
{
std::cout << "Degree = " << degree << std::endl;
const MappingQ<dim> mapping(degree);
const FE_Nothing<dim> fe;
FEValues<dim> fe_values(mapping, fe, quadrature, update_JxW_values);
for (unsigned int refinement = 0; refinement < 6;
++refinement, triangulation.refine_global(1))
{
dof_handler.distribute_dofs(fe);
double area = 0;
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
for (unsigned int i = 0; i < fe_values.n_quadrature_points; ++i)
area += fe_values.JxW(i);
}
table.add_value("eval.pi", area);
table.add_value("error", std::fabs(area - numbers::PI));
}
table.set_precision("eval.pi", 16);
table.set_scientific("error", true);
table.write_text(std::cout);
std::cout << std::endl;
}
}
template <int dim>
void compute_pi_by_perimeter()
{
std::cout << "Computation of Pi by the perimeter:" << std::endl
<< "===================================" << std::endl;
const QGauss<dim - 1> quadrature(4);
for (unsigned int degree = 1; degree < 5; ++degree)
{
std::cout << "Degree = " << degree << std::endl;
const MappingQ<dim> mapping(degree);
const FE_Nothing<dim> fe;
FEFaceValues<dim> fe_face_values(mapping,
fe,
quadrature,
for (unsigned int refinement = 0; refinement < 6;
++refinement, triangulation.refine_global(1))
{
dof_handler.distribute_dofs(fe);
double perimeter = 0;
for (const auto &cell : dof_handler.active_cell_iterators())
for (const auto &face : cell->face_iterators())
if (face->at_boundary())
{
fe_face_values.reinit(cell, face);
for (unsigned int i = 0;
i < fe_face_values.n_quadrature_points;
++i)
perimeter += fe_face_values.JxW(i);
}
table.add_value("eval.pi", static_cast<double>(perimeter / 2.0));
table.add_value("error", std::fabs(perimeter / 2.0 - numbers::PI));
}
table.set_precision("eval.pi", 16);
table.set_scientific("error", true);
table.write_text(std::cout);
std::cout << std::endl;
}
}
} // namespace Step10
int main()
{
try
{
std::cout.precision(16);
const unsigned int dim = 2;
Step10::gnuplot_output<dim>();
Step10::compute_pi_by_area<dim>();
Step10::compute_pi_by_perimeter<dim>();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
void evaluate_all_convergence_rates(const std::string &reference_column_key, const RateMode rate_mode)
void omit_column_from_convergence_rate_evaluation(const std::string &key)
void set_flags(const GridOutFlags::DX &flags)
Definition grid_out.cc:471
void write_gnuplot(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition grid_out.cc:4608
void write_text(std::ostream &out, const TextOutputFormat format=table_with_headers) const
void set_scientific(const std::string &key, const bool scientific)
void set_precision(const std::string &key, const unsigned int precision)