Reference documentation for deal.II version GIT relicensing-12-g89ba5bf9db 2024-02-24 00:00:03+00:00
\(\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-4 tutorial program

This tutorial depends on step-3.

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

Introduction

Note
The material presented here is also discussed in video lecture 12, video lecture 13. (All video lectures are also available here.)

deal.II has a unique feature which we call `‘dimension independent programming’'. You may have noticed in the previous examples that many classes had a number in angle brackets suffixed to them. This is to indicate that for example the triangulation in two and three space dimensions are different, but related data types. We could as well have called them Triangulation2d and Triangulation3d instead of Triangulation<2> and Triangulation<3> to name the two classes, but this has an important drawback: assume you have a function which does exactly the same functionality, but on 2d or 3d triangulations, depending on which dimension we would like to solve the equation in presently (if you don't believe that it is the common case that a function does something that is the same in all dimensions, just take a look at the code below - there are almost no distinctions between 2d and 3d!). We would have to write the same function twice, once working on Triangulation2d and once working with a Triangulation3d. This is an unnecessary obstacle in programming and leads to a nuisance to keep the two function in sync (at best) or difficult to find errors if the two versions get out of sync (at worst; this would probably the more common case).

Such obstacles can be circumvented by using some template magic as provided by the C++ language: templatized classes and functions are not really classes or functions but only a pattern depending on an as-yet undefined data type parameter or on a numerical value which is also unknown at the point of definition. However, the compiler can build proper classes or functions from these templates if you provide it with the information that is needed for that. Of course, parts of the template can depend on the template parameters, and they will be resolved at the time of compilation for a specific template parameter. For example, consider the following piece of code:

template <int dim>
{
};
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

At the point where the compiler sees this function, it does not know anything about the actual value of dim. The only thing the compiler has is a template, i.e. a blueprint, to generate functions make_grid if given a particular value of dim. Since dim has an unknown value, there is no code the compiler can generate for the moment.

However, if later down the compiler would encounter code that looks, for example, like this,

then the compiler will deduce that the function make_grid for dim==2 was requested and will compile the template above into a function with dim replaced by 2 everywhere, i.e. it will compile the function as if it were defined as

However, it is worth to note that the function GridGenerator::hyper_cube depends on the dimension as well, so in this case, the compiler will call the function GridGenerator::hyper_cube<2> while if dim were 3, it would call GridGenerator::hyper_cube<3> which might be (and actually is) a totally unrelated function.

The same can be done with member variables. Consider the following function, which might in turn call the above one:

template <int dim>
void make_grid_and_dofs (Triangulation<dim> &triangulation)
{
make_grid (triangulation);
...
};

This function has a member variable of type DoFHandler<dim>. Again, the compiler can't compile this function until it knows for which dimension. If you call this function for a specific dimension as above, the compiler will take the template, replace all occurrences of dim by the dimension for which it was called, and compile it. If you call the function several times for different dimensions, it will compile it several times, each time calling the right make_grid function and reserving the right amount of memory for the member variable; note that the size of a DoFHandler might, and indeed does, depend on the space dimension.

The deal.II library is built around this concept of dimension-independent programming, and therefore allows you to program in a way that will not need to distinguish between the space dimensions. It should be noted that in only a very few places is it necessary to actually compare the dimension using ifs or switches. However, since the compiler has to compile each function for each dimension separately, even there it knows the value of dim at the time of compilation and will therefore be able to optimize away the if statement along with the unused branch.

In this example program, we will show how to program dimension independently (which in fact is even simpler than if you had to take care about the dimension) and we will extend the Laplace problem of the last example to a program that runs in two and three space dimensions at the same time. Other extensions are the use of a non-constant right hand side function and of non-zero boundary values.

Note
When using templates, C++ imposes all sorts of syntax constraints that make it sometimes a bit difficult to understand why exactly something has to be written this way. A typical example is the need to use the keyword typename in so many places. If you are not entirely familiar with this already, then several of these difficulties are explained in the deal.II Frequently Asked Questions (FAQ) linked to from the deal.II homepage.

The commented program

Include files

The first few (many?) include files have already been used in the previous example, so we will not explain their meaning here again.

  #include <deal.II/grid/tria.h>
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/grid/grid_generator.h>
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/base/quadrature_lib.h>
  #include <deal.II/base/function.h>
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/numerics/matrix_tools.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/lac/full_matrix.h>
  #include <deal.II/lac/sparse_matrix.h>
  #include <deal.II/lac/dynamic_sparsity_pattern.h>
  #include <deal.II/lac/solver_cg.h>
  #include <deal.II/lac/precondition.h>
 
  #include <deal.II/numerics/data_out.h>
  #include <fstream>
  #include <iostream>
 

This is new, however: in the previous example we got some unwanted output from the linear solvers. If we want to suppress it, we have to include this file and add a single line somewhere to the program (see the main() function below for that):

  #include <deal.II/base/logstream.h>
 

The final step, as in previous programs, is to import all the deal.II class and function names into the global namespace:

  using namespace dealii;
 

The Step4 class template

This is again the same Step4 class as in the previous example. The only difference is that we have now declared it as a class with a template parameter, and the template parameter is of course the spatial dimension in which we would like to solve the Laplace equation. Of course, several of the member variables depend on this dimension as well, in particular the Triangulation class, which has to represent quadrilaterals or hexahedra, respectively. Apart from this, everything is as before.

  template <int dim>
  class Step4
  {
  public:
  Step4();
  void run();
 
  private:
  void make_grid();
  void setup_system();
  void assemble_system();
  void solve();
  void output_results() const;
 
  FE_Q<dim> fe;
  DoFHandler<dim> dof_handler;
 
  SparsityPattern sparsity_pattern;
  SparseMatrix<double> system_matrix;
 
  Vector<double> solution;
  Vector<double> system_rhs;
  };
 
 
Definition fe_q.h:550

Right hand side and boundary values

In the following, we declare two more classes denoting the right hand side and the non-homogeneous Dirichlet boundary values. Both are functions of a dim-dimensional space variable, so we declare them as templates as well.

Each of these classes is derived from a common, abstract base class Function, which declares the common interface which all functions have to follow. In particular, concrete classes have to overload the value function, which takes a point in dim-dimensional space as parameters and returns the value at that point as a double variable.

The value function takes a second argument, which we have here named component: This is only meant for vector-valued functions, where you may want to access a certain component of the vector at the point p. However, our functions are scalar, so we need not worry about this parameter and we will not use it in the implementation of the functions. Inside the library's header files, the Function base class's declaration of the value function has a default value of zero for the component, so we will access the value function of the right hand side with only one parameter, namely the point where we want to evaluate the function. A value for the component can then simply be omitted for scalar functions.

Function objects are used in lots of places in the library (for example, in step-3 we used a Functions::ZeroFunction instance as an argument to VectorTools::interpolate_boundary_values) and this is the first tutorial where we define a new class that inherits from Function. Since we only ever call Function::value(), we could get away with just a plain function (and this is what is done in step-5), but since this is a tutorial we inherit from Function for the sake of example.

  template <int dim>
  class RightHandSide : public Function<dim>
  {
  public:
  virtual double value(const Point<dim> &p,
  const unsigned int component = 0) const override;
  };
 
 
 
  template <int dim>
  class BoundaryValues : public Function<dim>
  {
  public:
  virtual double value(const Point<dim> &p,
  const unsigned int component = 0) const override;
  };
 
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Definition point.h:111

If you are not familiar with what the keywords virtual and override in the function declarations above mean, you will probably want to take a look at your favorite C++ book or an online tutorial such as http://www.cplusplus.com/doc/tutorial/polymorphism/ . In essence, what is happening here is that Function<dim> is an "abstract" base class that declares a certain "interface" – a set of functions one can call on objects of this kind. But it does not actually implement these functions: it just says "this is how Function objects look like", but what kind of function it actually is, is left to derived classes that implement the value() function.

Deriving one class from another is often called an "is-a" relationship function. Here, the RightHandSide class "is a" Function class because it implements the interface described by the Function base class. (The actual implementation of the value() function is in the code block below.) The virtual keyword then means "Yes, the function here is one that can be overridden by derived classes", and the override keyword means "Yes, this is in fact a function we know has been declared as part of the base class". The override keyword is not strictly necessary, but is an insurance against typos: If we get the name of the function or the type of one argument wrong, the compiler will warn us by stating "You say that this function overrides one in a base class, but I don't actually know any such function with this name and these arguments."

But back to the concrete case here: For this tutorial, we choose as right hand side the function \(4(x^4+y^4)\) in 2d, or \(4(x^4+y^4+z^4)\) in 3d. We could write this distinction using an if-statement on the space dimension, but here is a simple way that also allows us to use the same function in 1d (or in 4D, if you should desire to do so), by using a short loop. Fortunately, the compiler knows the size of the loop at compile time (remember that at the time when you define the template, the compiler doesn't know the value of dim, but when it later encounters a statement or declaration RightHandSide<2>, it will take the template, replace all occurrences of dim by 2 and compile the resulting function). In other words, at the time of compiling this function, the number of times the body will be executed is known, and the compiler can minimize the overhead needed for the loop; the result will be as fast as if we had used the formulas above right away.

The last thing to note is that a Point<dim> denotes a point in dim-dimensional space, and its individual components (i.e. \(x\), \(y\), ... coordinates) can be accessed using the () operator (in fact, the [] operator will work just as well) with indices starting at zero as usual in C and C++.

  template <int dim>
  double RightHandSide<dim>::value(const Point<dim> &p,
  const unsigned int /*component*/) const
  {
  double return_value = 0.0;
  for (unsigned int i = 0; i < dim; ++i)
  return_value += 4.0 * std::pow(p[i], 4.0);
 
  return return_value;
  }
 
 
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)

As boundary values, we choose \(x^2+y^2\) in 2d, and \(x^2+y^2+z^2\) in 3d. This happens to be equal to the square of the vector from the origin to the point at which we would like to evaluate the function, irrespective of the dimension. So that is what we return:

  template <int dim>
  double BoundaryValues<dim>::value(const Point<dim> &p,
  const unsigned int /*component*/) const
  {
  return p.square();
  }
 
 
 
constexpr numbers::NumberTraits< Number >::real_type square() const

Implementation of the Step4 class

Next for the implementation of the class template that makes use of the functions above. As before, we will write everything as templates that have a formal parameter dim that we assume unknown at the time we define the template functions. Only later, the compiler will find a declaration of Step4<2> (in the main function, actually) and compile the entire class with dim replaced by 2, a process referred to as ‘instantiation of a template’. When doing so, it will also replace instances of RightHandSide<dim> by RightHandSide<2> and instantiate the latter class from the class template.

In fact, the compiler will also find a declaration Step4<3> in main(). This will cause it to again go back to the general Step4<dim> template, replace all occurrences of dim, this time by 3, and compile the class a second time. Note that the two instantiations Step4<2> and Step4<3> are completely independent classes; their only common feature is that they are both instantiated from the same general template, but they are not convertible into each other, for example, and share no code (both instantiations are compiled completely independently).

Step4::Step4

After this introduction, here is the constructor of the Step4 class. It specifies the desired polynomial degree of the finite elements and associates the DoFHandler to the triangulation just as in the previous example program, step-3:

  template <int dim>
  Step4<dim>::Step4()
  : fe(1)
  , dof_handler(triangulation)
  {}
 
 

Step4::make_grid

Grid creation is something inherently dimension dependent. However, as long as the domains are sufficiently similar in 2d or 3d, the library can abstract for you. In our case, we would like to again solve on the square \([-1,1]\times [-1,1]\) in 2d, or on the cube \([-1,1] \times [-1,1] \times [-1,1]\) in 3d; both can be termed GridGenerator::hyper_cube(), so we may use the same function in whatever dimension we are. Of course, the functions that create a hypercube in two and three dimensions are very much different, but that is something you need not care about. Let the library handle the difficult things.

  template <int dim>
  void Step4<dim>::make_grid()
  {
  triangulation.refine_global(4);
 
  std::cout << " Number of active cells: " << triangulation.n_active_cells()
  << std::endl
  << " Total number of cells: " << triangulation.n_cells()
  << std::endl;
  }
 

Step4::setup_system

This function looks exactly like in the previous example, although it performs actions that in their details are quite different if dim happens to be 3. The only significant difference from a user's perspective is the number of cells resulting, which is much higher in three than in two space dimensions!

  template <int dim>
  void Step4<dim>::setup_system()
  {
  dof_handler.distribute_dofs(fe);
 
  std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
  << std::endl;
 
  DynamicSparsityPattern dsp(dof_handler.n_dofs());
  DoFTools::make_sparsity_pattern(dof_handler, dsp);
  sparsity_pattern.copy_from(dsp);
 
  system_matrix.reinit(sparsity_pattern);
 
  solution.reinit(dof_handler.n_dofs());
  system_rhs.reinit(dof_handler.n_dofs());
  }
 
 
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)

Step4::assemble_system

Unlike in the previous example, we would now like to use a non-constant right hand side function and non-zero boundary values. Both are tasks that are readily achieved with only a few new lines of code in the assemblage of the matrix and right hand side.

More interesting, though, is the way we assemble matrix and right hand side vector dimension independently: there is simply no difference to the two-dimensional case. Since the important objects used in this function (quadrature formula, FEValues) depend on the dimension by way of a template parameter as well, they can take care of setting up properly everything for the dimension for which this function is compiled. By declaring all classes which might depend on the dimension using a template parameter, the library can make nearly all work for you and you don't have to care about most things.

  template <int dim>
  void Step4<dim>::assemble_system()
  {
  QGauss<dim> quadrature_formula(fe.degree + 1);
 

We wanted to have a non-constant right hand side, so we use an object of the class declared above to generate the necessary data. Since this right hand side object is only used locally in the present function, we declare it here as a local variable:

  RightHandSide<dim> right_hand_side;
 

Compared to the previous example, in order to evaluate the non-constant right hand side function we now also need the quadrature points on the cell we are presently on (previously, we only required values and gradients of the shape function from the FEValues object, as well as the quadrature weights, FEValues::JxW() ). We can tell the FEValues object to do for us by also giving it the update_quadrature_points flag:

  FEValues<dim> fe_values(fe,
  quadrature_formula,
 
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.

We then again define the same abbreviation as in the previous program. The value of this variable of course depends on the dimension which we are presently using, but the FiniteElement class does all the necessary work for you and you don't have to care about the dimension dependent parts:

  const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 
  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
  Vector<double> cell_rhs(dofs_per_cell);
 
  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 

Next, we again have to loop over all cells and assemble local contributions. Note, that a cell is a quadrilateral in two space dimensions, but a hexahedron in 3d. In fact, the active_cell_iterator data type is something different, depending on the dimension we are in, but to the outside world they look alike and you will probably never see a difference. In any case, the real type is hidden by using auto:

  for (const auto &cell : dof_handler.active_cell_iterators())
  {
  fe_values.reinit(cell);
  cell_matrix = 0;
  cell_rhs = 0;
 

Now we have to assemble the local matrix and right hand side. This is done exactly like in the previous example, but now we revert the order of the loops (which we can safely do since they are independent of each other) and merge the loops for the local matrix and the local vector as far as possible to make things a bit faster.

Assembling the right hand side presents the only significant difference to how we did things in step-3: Instead of using a constant right hand side with value 1, we use the object representing the right hand side and evaluate it at the quadrature points:

  for (const unsigned int q_index : fe_values.quadrature_point_indices())
  for (const unsigned int i : fe_values.dof_indices())
  {
  for (const unsigned int j : fe_values.dof_indices())
  cell_matrix(i, j) +=
  (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
  fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
  fe_values.JxW(q_index)); // dx
 
  const auto &x_q = fe_values.quadrature_point(q_index);
  cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
  right_hand_side.value(x_q) * // f(x_q)
  fe_values.JxW(q_index)); // dx
  }

As a final remark to these loops: when we assemble the local contributions into cell_matrix(i,j), we have to multiply the gradients of shape functions \(i\) and \(j\) at point number q_index and multiply it with the scalar weights JxW. This is what actually happens: fe_values.shape_grad(i,q_index) returns a dim dimensional vector, represented by a Tensor<1,dim> object, and the operator* that multiplies it with the result of fe_values.shape_grad(j,q_index) makes sure that the dim components of the two vectors are properly contracted, and the result is a scalar floating point number that then is multiplied with the weights. Internally, this operator* makes sure that this happens correctly for all dim components of the vectors, whether dim be 2, 3, or any other space dimension; from a user's perspective, this is not something worth bothering with, however, making things a lot simpler if one wants to write code dimension independently.

With the local systems assembled, the transfer into the global matrix and right hand side is done exactly as before, but here we have again merged some loops for efficiency:

  cell->get_dof_indices(local_dof_indices);
  for (const unsigned int i : fe_values.dof_indices())
  {
  for (const unsigned int j : fe_values.dof_indices())
  system_matrix.add(local_dof_indices[i],
  local_dof_indices[j],
  cell_matrix(i, j));
 
  system_rhs(local_dof_indices[i]) += cell_rhs(i);
  }
  }
 

As the final step in this function, we wanted to have non-homogeneous boundary values in this example, unlike the one before. This is a simple task, we only have to replace the Functions::ZeroFunction used there by an object of the class which describes the boundary values we would like to use (i.e. the BoundaryValues class declared above):

The function VectorTools::interpolate_boundary_values() will only work on faces that have been marked with boundary indicator 0 (because that's what we say the function should work on with the second argument below). If there are faces with boundary id other than 0, then the function interpolate_boundary_values will do nothing on these faces. For the Laplace equation doing nothing is equivalent to assuming that on those parts of the boundary a zero Neumann boundary condition holds.

  std::map<types::global_dof_index, double> boundary_values;
  0,
  BoundaryValues<dim>(),
  boundary_values);
  system_matrix,
  solution,
  system_rhs);
  }
 
 
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})

Step4::solve

Solving the linear system of equations is something that looks almost identical in most programs. In particular, it is dimension independent, so this function is copied verbatim from the previous example.

  template <int dim>
  void Step4<dim>::solve()
  {
  SolverControl solver_control(1000, 1e-12);
  SolverCG<Vector<double>> solver(solver_control);
  solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
 

We have made one addition, though: since we suppress output from the linear solvers, we have to print the number of iterations by hand.

  std::cout << " " << solver_control.last_step()
  << " CG iterations needed to obtain convergence." << std::endl;
  }
 
 

Step4::output_results

This function also does what the respective one did in step-3. No changes here for dimension independence either.

Since the program will run both 2d and 3d versions of the Laplace solver, we use the dimension in the filename to generate distinct filenames for each run (in a better program, one would check whether dim can have other values than 2 or 3, but we neglect this here for the sake of brevity).

  template <int dim>
  void Step4<dim>::output_results() const
  {
  DataOut<dim> data_out;
 
  data_out.attach_dof_handler(dof_handler);
  data_out.add_data_vector(solution, "solution");
 
  data_out.build_patches();
 
  std::ofstream output(dim == 2 ? "solution-2d.vtk" : "solution-3d.vtk");
  data_out.write_vtk(output);
  }
 
 
 
void attach_dof_handler(const DoFHandler< dim, spacedim > &)

Step4::run

This is the function which has the top-level control over everything. Apart from one line of additional output, it is the same as for the previous example.

  template <int dim>
  void Step4<dim>::run()
  {
  std::cout << "Solving problem in " << dim << " space dimensions."
  << std::endl;
 
  make_grid();
  setup_system();
  assemble_system();
  solve();
  output_results();
  }
 
 

The main function

And this is the main function. It also looks mostly like in step-3, but if you look at the code below, note how we first create a variable of type Step4<2> (forcing the compiler to compile the class template with dim replaced by 2) and run a 2d simulation, and then we do the whole thing over in 3d.

In practice, this is probably not what you would do very frequently (you probably either want to solve a 2d problem, or one in 3d, but not both at the same time). However, it demonstrates the mechanism by which we can simply change which dimension we want in a single place, and thereby force the compiler to recompile the dimension independent class templates for the dimension we request. The emphasis here lies on the fact that we only need to change a single place. This makes it rather trivial to debug the program in 2d where computations are fast, and then switch a single place to a 3 to run the much more computing intensive program in 3d for ‘real’ computations.

Each of the two blocks is enclosed in braces to make sure that the laplace_problem_2d variable goes out of scope (and releases the memory it holds) before we move on to allocate memory for the 3d case. Without the additional braces, the laplace_problem_2d variable would only be destroyed at the end of the function, i.e. after running the 3d problem, and would needlessly hog memory while the 3d run could actually use it.

  int main()
  {
  {
  Step4<2> laplace_problem_2d;
  laplace_problem_2d.run();
  }
 
  {
  Step4<3> laplace_problem_3d;
  laplace_problem_3d.run();
  }
 
  return 0;
  }

Results

The output of the program looks as follows (the number of iterations may vary by one or two, depending on your computer, since this is often dependent on the round-off accuracy of floating point operations, which differs between processors):

Solving problem in 2 space dimensions.
Number of active cells: 256
Total number of cells: 341
Number of degrees of freedom: 289
26 CG iterations needed to obtain convergence.
Solving problem in 3 space dimensions.
Number of active cells: 4096
Total number of cells: 4681
Number of degrees of freedom: 4913
30 CG iterations needed to obtain convergence.

It is obvious that in three spatial dimensions the number of cells and therefore also the number of degrees of freedom is much higher. What cannot be seen here, is that besides this higher number of rows and columns in the matrix, there are also significantly more entries per row of the matrix in three space dimensions. Together, this leads to a much higher numerical effort for solving the system of equation, which you can feel in the run time of the two solution steps when you actually run the program.

The program produces two files: solution-2d.vtk and solution-3d.vtk, which can be viewed using the programs VisIt or Paraview (in case you do not have these programs, you can easily change the output format in the program to something which you can view more easily). Visualizing solutions is a bit of an art, but it can also be fun, so you should play around with your favorite visualization tool to get familiar with its functionality. Here's what I have come up with for the 2d solution:

(See also video lecture 11, video lecture 32.) The picture shows the solution of the problem under consideration as a 3D plot. As can be seen, the solution is almost flat in the interior of the domain and has a higher curvature near the boundary. This, of course, is due to the fact that for Laplace's equation the curvature of the solution is equal to the right hand side and that was chosen as a quartic polynomial which is nearly zero in the interior and is only rising sharply when approaching the boundaries of the domain; the maximal values of the right hand side function are at the corners of the domain, where also the solution is moving most rapidly. It is also nice to see that the solution follows the desired quadratic boundary values along the boundaries of the domain. It can also be useful to verify a computed solution against an analytical solution. For an explanation of this technique, see step-7.

On the other hand, even though the picture does not show the mesh lines explicitly, you can see them as little kinks in the solution. This clearly indicates that the solution hasn't been computed to very high accuracy and that to get a better solution, we may have to compute on a finer mesh.

In three spatial dimensions, visualization is a bit more difficult. The left picture shows the solution and the mesh it was computed on on the surface of the domain. This is nice, but it has the drawback that it completely hides what is happening on the inside. The picture on the right is an attempt at visualizing the interior as well, by showing surfaces where the solution has constant values (as indicated by the legend at the top left). Isosurface pictures look best if one makes the individual surfaces slightly transparent so that it is possible to see through them and see what's behind.

Note
A final remark on visualization: the idea of visualization is to give insight, which is not the same as displaying information. In particular, it is easy to overload a picture with information, but while it shows more information it makes it also more difficult to glean insight. As an example, the program I used to generate these pictures, VisIt, by default puts tick marks on every axis, puts a big fat label "X Axis" on the \(x\) axis and similar for the other axes, shows the file name from which the data was taken in the top left and the name of the user doing so and the time and date on the bottom right. None of this is important here: the axes are equally easy to make out because the tripod at the bottom left is still visible, and we know from the program that the domain is \([-1,1]^3\), so there is no need for tick marks. As a consequence, I have switched off all the extraneous stuff in the picture: the art of visualization is to reduce the picture to those parts that are important to see what one wants to see, but no more.

Postprocessing: What to do with the solution?

This tutorial – like most of the other programs – principally only shows how to numerically approximate the solution of a partial differential equation, and then how to visualize this solution graphically. But solving a PDE is of course not the goal in most practical applications (unless you are a numerical methods developer and the method is the goal): We generally want to solve a PDE because we want to extract information from it. Examples for what people are interested in from solutions include the following:

  • Let's say you solve the equations of elasticity (which we will do in step-8), then that's presumably because you want to know about the deformation of an elastic object under a given load. From an engineering perspective, what you then presumably want to learn is the degree of deformation of the object, say at a specific point; or you may want to know the maximum stress in order to determine whether the applied load exceeds the safe maximal stress the material can withstand.
  • If you are solving fluid flow problems (such as in step-22, step-57, step-67, and step-69), then you might be interested in the fluid velocity at specific points, and oftentimes the forces the fluid exerts on the boundary of the fluid domain. The latter is important in many applications: If the fluid in question is the air flowing around an airplane, then we are typically interested in the drag and lift forces on the fuselage and wings. If the fluid is water flowing around a ship, then we typically care about the drag force on the ship.
  • If you are solving the Maxwell equations of electromagnetics, you are typically interested in how much energy is radiated in certain directions (say, in order to know the range of information transmission via an antenna, or to determine the radar cross section of planes or ships).

The point here is that from an engineering perspective, solving the PDE is only the first step. The second step is to evaluate the computed solution in order to extract relevant numbers that allow us to either optimize a design, or to make decisions. This second step is often called "postprocessing the solution".

This program does not solve a solid or fluid mechanics problem, so we should try to illustrate postprocessing with something that makes sense in the context of the equation we solve here. The Poisson equation in two space dimensions is a model for the vertical deformation of a membrane that is clamped at the boundary and is subject to a vertical force. For this kind of situation, it makes sense to evaluate the average vertical displacement,

\[ \bar u_h = \frac{\int_\Omega u_h(\mathbf x) \, dx}{|\Omega|}, \]

where \(|\Omega| = \int_\Omega 1 \, dx\) is the area of the domain. To compute \(\bar u_h\), as usual we replace integrals over the domain by a sum of integrals over cells,

\[ \int_\Omega u_h(\mathbf x) \, dx = \sum_K \int_K u_h(\mathbf x) \, dx, \]

and then integrals over cells are approximated by quadrature:

\begin{align*} \int_\Omega u_h(\mathbf x) \, dx &= \sum_K \int_K u_h(\mathbf x) \, dx, \\ &= \sum_K \sum_q u_h(\mathbf x_q^K) w_q^K, \end{align*}

where \(w_q^K\) is the weight of the \(q\)th quadrature point evaluated on cell \(K\). All of this is as always provided by the FEValues class – the entry point for all integrals in deal.II.

The actual implementation of this is straightforward once you know how to get the values of the solution \(u\) at the quadrature points of a cell. This functionality is provided by FEValues::get_function_values(), a function that takes a global vector of nodal values as input and returns a vector of function values at the quadrature points of the current cell. Using this function, to see how it all works together you can place the following code snippet anywhere in the program after the solution has been computed (the output_results() function seems like a good place to also do postprocessing, for example):

QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
std::vector<double> solution_values(quadrature_formula.size());
double integral_of_u = 0;
double volume_of_omega = 0;
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
fe_values.get_function_values(solution, solution_values);
for (const unsigned int q_point : fe_values.quadrature_point_indices())
{
integral_of_u += solution_values[q_point] * fe_values.JxW(q_point);
volume_of_omega += 1 * fe_values.JxW(q_point);
}
}
std::cout << " Mean value of u=" << integral_of_u / volume_of_omega
<< std::endl;

In this code snippet, we also compute the volume (or, since we are currently thinking about a two-dimensional situation: the area) \(|\Omega|\) by computing the integral \(|\Omega| = \int_\Omega 1 \, dx\) in exactly the same way, via quadrature. (We could avoid having to compute \(|\Omega|\) by hand here, using the fact that deal.II has a function for this: GridTools::volume(). That said, it is efficient to compute the two integrals concurrently in the same loop, and so that's what we do.)

This program of course also solves the same Poisson equation in three space dimensions. In this situation, the Poisson equation is often used as a model for diffusion of either a physical species (say, of ink in a tank of water, or a pollutant in the air) or of energy (specifically, of thermal energy in a solid body). In that context, the quantity

\[ \Phi_h = \int_{\partial\Omega} \nabla u_h(\mathbf x) \cdot \mathbf n(\mathbf x) \; dx \]

is the flux of this species or energy across the boundary. (In actual physical models, one would also have to multiply the right hand side by a diffusivity or conductivity constant, but let us ignore this here.) In much the same way as before, we compute such integrals by splitting it over integrals of faces of cells, and then applying quadrature:

\begin{align*} \Phi_h &= \int_{\partial\Omega} \nabla u_h(\mathbf x) \cdot \mathbf n(\mathbf x) \; dx \\ &= \sum_K \sum_{f \in \text{faces of @f$K@f$}, f\subset\partial\Omega} \int_f \nabla u_h(\mathbf x) \cdot \mathbf n(\mathbf x) \; dx \\ &= \sum_K \sum_{f \in \text{faces of @f$K@f$}, f\subset\partial\Omega} \sum_q \nabla u_h(\mathbf x_q^f) \cdot \mathbf n(\mathbf x_q^f) w_q^f, \end{align*}

where now \(\mathbf x_q^f\) are the quadrature points located on face \(f\), and \(w_q^f\) are the weights associated with these faces. The second of the sum symbols loops over all faces of cell \(K\), but restricted to those that are actually at the boundary.

This all is easily implemented by the following code that replaces the use of the FEValues class (which is used for integrating over cells – i.e., domain integrals) by the FEFaceValues class (which is used for integrating over faces – i.e., boundary integrals):

QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
FEFaceValues<dim> fe_face_values(fe,
face_quadrature_formula,
std::vector<Tensor<1, dim>> solution_gradients(face_quadrature_formula.size());
double flux = 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);
fe_face_values.get_function_gradients(solution, solution_gradients);
for (const unsigned int q_point :
fe_face_values.quadrature_point_indices())
{
flux += solution_gradients[q_point] *
fe_face_values.normal_vector(q_point) *
fe_face_values.JxW(q_point);
}
}
std::cout << " Flux=" << flux << std::endl;
@ update_normal_vectors
Normal vectors.

If you add these two code snippets to the code, you will get output like the following when you run the program:

Solving problem in 2 space dimensions.
Number of active cells: 256
Total number of cells: 341
Number of degrees of freedom: 289
26 CG iterations needed to obtain convergence.
Mean value of u=1.33303
Flux=-3.68956
Solving problem in 3 space dimensions.
Number of active cells: 4096
Total number of cells: 4681
Number of degrees of freedom: 4913
30 CG iterations needed to obtain convergence.
Mean value of u=1.58058
Flux=-8.29435

This makes some sense: If you look, for example, at the 2d output above, the solution varies between values of 1 and 2, but with a larger part of the solution closer to one than two; so an average value of 1.33 for the mean value is reasonable. For the flux, recall that \(\nabla u \cdot \mathbf n\) is the directional derivative in the normal direction – in other words, how the solution changes as we move from the interior of the domain towards the boundary. If you look at the 2d solution, you will realize that for most parts of the boundary, the solution decreases as we approach the boundary, so the normal derivative is negative – so if we integrate along the boundary, we should expect (and obtain!) a negative value.

Possibilities for extensions

There are many ways with which one can play with this program. The simpler ones include essentially all the possibilities already discussed in the Possibilities for extensions in the documentation of step 3, except that you will have to think about whether something now also applies to the 3d case discussed in the current program.

It is also worthwhile considering the postprocessing options discussed above. The documentation states two numbers (the mean value and the normal flux) for both the 2d and 3d cases. Can we trust these numbers? We have convinced ourselves that at least the mean value is reasonable, and that the sign of the flux is probably correct. But are these numbers accurate?

A general rule is that we should never trust a number unless we have verified it in some way. From the theory of finite element methods, we know that as we make the mesh finer and finer, the numerical solution \(u_h\) we compute here must converge to the exact solution \(u\). As a consequence, we also expect that \(\bar u_h \rightarrow \bar u\) and \(\Phi_h \rightarrow \Phi\), but that does not mean that for any given mesh \(\bar u_h\) or \(\Phi_h\) are particularly accurate approximations.

To test this kind of thing, we have already considered the convergence of a point value in step-3. We can do the same here by selecting how many times the mesh is globally refined in the make_grid() function of this program. For the mean value of the solution, we then get the following numbers:

# of refinements \(\bar u_h\) in 2d \(\bar u_h\) in 3d
4 1.33303 1.58058
5 1.33276 1.57947
6 1.3327 1.5792
7 1.33269 1.57914
8 1.33268
9 1.33268

I did not have the patience to run the last two values for the 3d case – one needs quite a fine mesh for this, with correspondingly long run times. But we can be reasonably assured that values around 1.33 (for the 2d case) and 1.58 (for the 3d case) are about right – and at least for engineering applications, three digits of accuracy are good enough.

The situation looks very different for the flux. Here, we get results such as the following:

# of refinements \(\Phi_h\) in 2d \(\Phi_h\) in 3d
4 -3.68956 -8.29435
5 -4.90147 -13.0691
6 -5.60745 -15.9171
7 -5.99111 -17.4918
8 -6.19196
9 -6.29497
10 -6.34721
11 -6.37353

So this is not great. For the 2d case, we might infer that perhaps a value around -6.4 might be right if we just refine the mesh enough – though 11 refinements already leads to some 4,194,304 cells. In any case, the first number (the one shown in the beginning where we discussed postprocessing) was off by almost a factor of 2!

For the 3d case, the last number shown was on a mesh with 2,097,152 cells; the next one would have had 8 times as many cells. In any case, the numbers mean that we can't even be sure that the first digit of that last number is correct! In other words, it was worth checking, or we would have just believed all of these numbers. In fact, that last column isn't even doing a particularly good job convincing us that the code might be correctly implemented.

If you keep reading through the other tutorial programs, you will find many ways to make these sorts of computations more accurate and to come to believe that the flux actually does converge to its correct value. For example, we can dramatically increase the accuracy of the computation by using adaptive mesh refinement (step-6) near the boundary, and in particular by using higher polynomial degree finite elements (also step-6, but also step-7). Using the latter, using cubic elements (polynomial degree 3), we can actually compute the flux pretty accurately even in 3d: \(\Phi_h=-19.0148\) with 4 global refinement steps, and \(\Phi_h=-19.1533\) with 5 refinement steps. These numbers are already pretty close together and give us a reasonable idea of the first two correct digits of the "true" answer.

Note
We would be remiss to not also comment on the fact that there are good theoretical reasons why computing the flux accurately appears to be so much more difficult than the average value. This has to do with the fact that finite element theory provides us with the estimate \(\|u-u_h\|_{L_2(\Omega)} \le C h^2 \|\nabla^2u\|_{L_2(\Omega)}\) when using the linear elements this program uses – that is, for every global mesh refinement, \(h\) is reduced by a factor of two and the error goes down by a factor of 4. Now, the \(L_2\) error is not equivalent to the error in the mean value, but the two are related: They are both integrals over the domain, using the value of the solution. We expect the mean value to converge no worse than the \(L_2\) norm of the error. At the same time, theory also provides us with this estimate: \(\|\nabla (u-u_h)\|_{L_2(\partial\Omega)} \le C h^{1/2} \|\nabla^2u\|_{L_2(\Omega)}\). The move from values to gradients reduces the convergence rates by one order, and the move from domain to boundary by another half order. Here, then, each refinement step reduces the error not by a factor of 4 any more, by only by a factor of \(\sqrt{2} \approx 1.4\). It takes a lot of global refinement steps to reduce the error by, say, a factor ten or hundred, and this is reflected in the very slow convergence evidenced by the table. On the other hand, for cubic elements (i.e., polynomial degree 3), we would get \(\|u-u_h\|_{L_2(\Omega)} \le C h^4 \|\nabla^4u\|_{L_2(\Omega)}\) and after reduction by 1.5 orders, we would still have \(\|\nabla (u-u_h)\|_{L_2(\partial\Omega)} \le C h^{2+1/2} \|\nabla^4u\|_{L_2(\Omega)}\). This rate, \({\cal O}(h^{2.5})\) is still quite rapid, and it is perhaps not surprising that we get much better answers with these higher order elements. This also illustrates that when trying to approximate anything that relates to a gradient of the solution, using linear elements (polynomial degree one) is really not a good choice at all.
In this very specific case, it turns out that we can actually compute the exact value of \(\Phi\). This is because for the Poisson equation we compute the solution of here, \(-\Delta u = f\), we can integrate over the domain, \(-\int_\Omega \Delta u = \int_\Omega f\), and then use that \(\Delta = \text{div}\;\text{grad}\); this allows us to use the divergence theorem followed by multiplying by minus one to find \(\int_{\partial\Omega} \nabla u \cdot n = -\int_\Omega f\). The left hand side happens to be \(\Phi\). For the specific right hand side \(f(x_1,x_2)=4(x_1^4+x_2^4)\) we use in 2d, we then get \(-\int_\Omega f = -\int_{-1}^{1} \int_{-1}^{1} 4(x_1^4+x_2^4) \; dx_2\; dx_1 = -16 \left[\int_{-1}^{1} x^4 \; dx\right] = -16\times\frac 25\), which has a numerical value of exactly -6.4 – right on with our guess above. In 3d, we can do the same and get that the exact value is \(-\int_\Omega f = -\int_{-1}^{1} \int_{-1}^{1} \int_{-1}^{1} 4(x_1^4+x_2^4+x_3^4) \; dx_3 \; dx_2\; dx_1 = -48\times\frac 25=-19.2\). What we found with cubic elements is then quite close to this exact value. Of course, in practice we almost never have exact values to compare with: If we could compute something on a piece of paper, we wouldn't have to solve the PDE numerically. But these sorts of situations make for excellent test cases that help us verify that our numerical solver works correctly. In many other cases, the literature contains numbers where others have already computed an answer accurately using their own software, and these are also often useful to compare against in verifying the correctness of our codes.

The plain program

/* ------------------------------------------------------------------------
*
* SPDX-License-Identifier: LGPL-2.1-or-later
* Copyright (C) 1999 - 2024 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.
*
* ------------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, University of Heidelberg, 1999
*/
#include <fstream>
#include <iostream>
using namespace dealii;
template <int dim>
class Step4
{
public:
Step4();
void run();
private:
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results() const;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
template <int dim>
class RightHandSide : public Function<dim>
{
public:
virtual double value(const Point<dim> &p,
const unsigned int component = 0) const override;
};
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
virtual double value(const Point<dim> &p,
const unsigned int component = 0) const override;
};
template <int dim>
double RightHandSide<dim>::value(const Point<dim> &p,
const unsigned int /*component*/) const
{
double return_value = 0.0;
for (unsigned int i = 0; i < dim; ++i)
return_value += 4.0 * std::pow(p[i], 4.0);
return return_value;
}
template <int dim>
double BoundaryValues<dim>::value(const Point<dim> &p,
const unsigned int /*component*/) const
{
return p.square();
}
template <int dim>
Step4<dim>::Step4()
: fe(1)
, dof_handler(triangulation)
{}
template <int dim>
void Step4<dim>::make_grid()
{
triangulation.refine_global(4);
std::cout << " Number of active cells: " << triangulation.n_active_cells()
<< std::endl
<< " Total number of cells: " << triangulation.n_cells()
<< std::endl;
}
template <int dim>
void Step4<dim>::setup_system()
{
dof_handler.distribute_dofs(fe);
std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
DynamicSparsityPattern dsp(dof_handler.n_dofs());
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
template <int dim>
void Step4<dim>::assemble_system()
{
QGauss<dim> quadrature_formula(fe.degree + 1);
RightHandSide<dim> right_hand_side;
FEValues<dim> fe_values(fe,
quadrature_formula,
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell_rhs = 0;
for (const unsigned int q_index : fe_values.quadrature_point_indices())
for (const unsigned int i : fe_values.dof_indices())
{
for (const unsigned int j : fe_values.dof_indices())
cell_matrix(i, j) +=
(fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
fe_values.JxW(q_index)); // dx
const auto &x_q = fe_values.quadrature_point(q_index);
cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
right_hand_side.value(x_q) * // f(x_q)
fe_values.JxW(q_index)); // dx
}
cell->get_dof_indices(local_dof_indices);
for (const unsigned int i : fe_values.dof_indices())
{
for (const unsigned int j : fe_values.dof_indices())
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
std::map<types::global_dof_index, double> boundary_values;
0,
BoundaryValues<dim>(),
boundary_values);
system_matrix,
solution,
system_rhs);
}
template <int dim>
void Step4<dim>::solve()
{
SolverControl solver_control(1000, 1e-12);
SolverCG<Vector<double>> solver(solver_control);
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
std::cout << " " << solver_control.last_step()
<< " CG iterations needed to obtain convergence." << std::endl;
}
template <int dim>
void Step4<dim>::output_results() const
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");
data_out.build_patches();
std::ofstream output(dim == 2 ? "solution-2d.vtk" : "solution-3d.vtk");
data_out.write_vtk(output);
}
template <int dim>
void Step4<dim>::run()
{
std::cout << "Solving problem in " << dim << " space dimensions."
<< std::endl;
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
int main()
{
{
Step4<2> laplace_problem_2d;
laplace_problem_2d.run();
}
{
Step4<3> laplace_problem_3d;
laplace_problem_3d.run();
}
return 0;
}
void write_vtk(std::ostream &out) const
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition data_out.cc:1062
spacedim const Point< spacedim > & p
Definition grid_tools.h:981
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74