Basic concept

To define the shape functions, we need the following:

• A mesh on which to define shape functions. You have already seen how to generate and manipulate the objects that describe meshes in step-1 and step-2.
• A finite element that describes the shape functions we want to use on the reference cell (which in deal.II is always the unit interval \([0,1]\), the unit square \([0,1]^{2}\) or the unit cube \([0,1]^{3}\), depending on which space dimension you work in). In step-2, we had already used an object of type \(F E_{1} Q<2>\), which denotes the usual Lagrange elements that define shape functions by interpolation on support points The simplest one is FE_Q \(2>(1)\), which uses polynomial degree 1. In 2d, these are often referred to as bilinear, since they are linear in each of the two coordinates of the reference cell. (In \(1 d\), they would be linear and in \(3 \mathrm{~d}\) tri-linear, however, in the deal.II documentation, we will frequently not make this distinction and simply always call these functions "linear".)
• A DoFHandler object that enumerates all the degrees of freedom on the mesh, taking the reference cell description the finite element object provides as the basis. You've also already seen how to do this in step-2.
• A mapping that tells how the shape functions on the real cell are obtained from the shape functions defined by the finite element class on the reference cell. By default, unless you explicitly say otherwise, deal.II will use a (bi-, tri-)linear mapping for this, so in most cases you don't have to worry about this step.

Define the dimension specifically:

Triangulation<3> or Triangulation<2>

The basic class Function, which declares the common interface which all functions have to follow, would concrete classes that 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. For example:

1
2
3
4
5
6
7
8
9
10
11
12
13
14

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;
};

In essence, what is happening here is that Function 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.

Solve the wave equation

The wave equation in its prototypical form reads as follows: find \(u(x, t), x \in \Omega, t \in[0, T]\) that satisfies \[ \begin{aligned} \frac{\partial^{2} u}{\partial t^{2}}-\Delta u &=f & & \text { in } \Omega \times[0, T] \\ u(x, t) &=g & & \text { on } \partial \Omega \times[0, T] \\ u(x, 0) &=u_{0}(x) & & \text { in } \Omega \\ \frac{\partial u(x, 0)}{\partial t} &=u_{1}(x) & & \text { in } \Omega \end{aligned} \] Note that since this is an equation with second-order time derivatives, we need to pose two initial conditions, one for the value and one for the time derivative of the solution.

Rothe's method

Given these considerations, here is how we will proceed: let us first define a simple time stepping method for this second order problem, and then in a second step do the spatial discretization, i.e. we will follow Rothe's approach.

For the first step, let us take a little detour first: in order to discretize a second time derivative, we can either discretize it directly, or we can introduce an additional variable and transform the system into a first order system. In many cases, this turns out to be equivalent, but dealing with first order systems is often simpler. To this end, let us introduce \[ v=\frac{\partial u}{\partial t}, \] and call this variable the velocity for obvious reasons. We can then reformulate the original wave equation as follows: \[ \begin{aligned} \frac{\partial u}{\partial t}-v &=0 & & \text { in } \Omega \times[0, T], \\ \frac{\partial \nu}{\partial t}-\Delta u &=f & & \text { in } \Omega \times[0, T], \\ u(x, t) &=g & & \text { on } \partial \Omega \times[0, T], \\ u(x, 0) &=u_{0}(x) & & \text { in } \Omega, \\ v(x, 0) &=u_{1}(x) & & \text { in } \Omega . \end{aligned} \] The advantage of this formulation is that it now only contains first time derivatives for both variables, for which it is simple to write down time stepping schemes. Note that we do not have boundary conditions for \(v\) at first. However, we could enforce \(v=\frac{\partial g}{\partial t}\) on the boundary. It turns out in numerical examples that this is actually necessary: without doing so the solution doesn't look particularly wrong, but the CrankNicolson scheme does not conserve energy if one doesn't enforce these boundary conditions. With this formulation, let us introduce the following time discretization where a superscript \(n\) indicates the number of a time step and \(k=t_{n}-t_{n-1}\) is the length of the present time step: \[ \begin{aligned} &\frac{u^{n}-u^{n-1}}{k}-\left[\theta v^{n}+(1-\theta) v^{n-1}\right]=0, \\ &\frac{v^{n}-v^{n-1}}{k}-\Delta\left[\theta u^{n}+(1-\theta) u^{n-1}\right]=\theta f^{n}+(1-\theta) f^{n-1} . \end{aligned} \] Note how we introduced a parameter \(\theta\) here. If we chose \(\theta=0\), for example, the first equation would reduce to \(\frac{u^{n}-u^{n-1}}{k}-v^{n-1}=0\), which is well-known as the forward or explicit Euler method. On the other hand, if we set \(\theta=1\), then we would get \(\frac{u^{n}-u^{n-1}}{k}-v^{n}=0\), which corresponds to the backward or implicit Euler method. Both these methods are first order accurate methods. They are simple to implement, but they are not really very accurate.

The third case would be to choose \(\theta=\frac{1}{2}\). The first of the equations above would then read \(\frac{u^{n}-u^{n-1}}{k}-\frac{1}{2}\left[v^{n}+v^{n-1}\right]=0\). This method is known as the Crank-Nicolson method and has the advantage that it is second order accurate. In addition, it has the nice property that it preserves the energy in the solution (physically, the energy is the sum of the kinetic energy of the particles in the membrane plus the potential energy present due to the fact that it is locally stretched; this quantity is a conserved one in the continuous equation, but most time stepping schemes do not conserve it after time discretization). Since \(v^{n}\) also appears in the equation for \(u^{n}\), the Crank-Nicolson scheme is also implicit. In the program, we will leave \(\theta\) as a parameter, so that it will be easy to play with it. The results section will show some numerical evidence comparing the different schemes. The equations above (called the semidiscretized equations because we have only discretized the time, but not space), can be simplified a bit by eliminating \(v^{n}\) from the first equation and rearranging terms. We then get \[ \begin{aligned} \left[1-k^{2} \theta^{2} \Delta\right] u^{n} &=\left[1+k^{2} \theta(1-\theta) \Delta\right] u^{n-1}+k v^{n-1}+k^{2} \theta\left[\theta f^{n}+(1-\theta) f^{n-1}\right] \\ v^{n} &=v^{n-1}+k \Delta\left[\theta u^{n}+(1-\theta) u^{n-1}\right]+k\left[\theta f^{n}+(1-\theta) f^{n-1}\right] \end{aligned} \] In this form, we see that if we are given the solution \(u^{n-1}, v^{n-1}\) of the previous timestep, that we can then solve for the variables \(u^{n}, v^{n}\) separately, i.e. one at a time. This is convenient. In addition, we recognize that the operator in the first equation is positive definite, and the second equation looks particularly simple.

Space discretization

We have now derived equations that relate the approximate (semi-discrete) solution \(u^{n}(x)\) and its time derivative \(v^{n}(x)\) at time \(t_{n}\) with the solutions \(u^{n-1}(x), v^{n-1}(x)\) of the previous time step at \(t_{n-1}\). The next step is to also discretize the spatial variable using the usual finite element methodology. To this end, we multiply each equation with a test function, integrate over the entire domain, and integrate by parts where necessary. This leads to \[ \begin{aligned} \left(u^{n}, \varphi\right)+k^{2} \theta^{2}\left(\nabla u^{n}, \nabla \varphi\right) &=\left(u^{n-1}, \varphi\right)-k^{2} \theta(1-\theta)\left(\nabla u^{n-1}, \nabla \varphi\right)+k\left(v^{n-1}, \varphi\right)+k^{2} \theta\left[\theta\left(f^{n}, \varphi\right)+(1-\theta)\left(f^{n-1}, \varphi\right)\right], \\ \left(v^{n}, \varphi\right) &=\left(v^{n-1}, \varphi\right)-k\left[\theta\left(\nabla u^{n}, \nabla \varphi\right)+(1-\theta)\left(\nabla u^{n-1}, \nabla \varphi\right)\right]+k\left[\theta\left(f^{n}, \varphi\right)+(1-\theta)\left(f^{n-1}, \varphi\right)\right] . \end{aligned} \] It is then customary to approximate \(u^{n}(x) \approx u_{h}^{n}(x)=\sum_{i} U_{i}^{n} \phi_{i}^{n}(x)\), where \(\phi_{i}^{n}(x)\) are the shape functions used for the discretization of the \(n\)-th time step and \(U_{i}^{n}\) are the unknown nodal values of the solution. Similarly, \(v^{n}(x) \approx v_{h}^{n}(x)=\sum_{i} V_{i}^{n} \phi_{i}^{n}(x)\). Finally, we have the solutions of the previous time \(\operatorname{step}, u^{n-1}(x) \approx u_{h}^{n-1}(x)=\sum_{i} U_{i}^{n-1} \phi_{i}^{n-1}(x)\) and \(v^{n-1}(x) \approx v_{h}^{n-1}(x)=\sum_{i} V_{i}^{n-1} \phi_{i}^{n-1}(x)\). Note that since the solution of the previous time step has already been computed by the time we get to time step \(n, U^{n-1}, V^{n-1}\) are known. Furthermore, note that the solutions of the previous step may have been computed on a different mesh, so we have to use shape functions \(\phi_{i}^{n-1}(x)\). If we plug these expansions into above equations and test with the test functions from the present mesh, we get the following linear system: \[ \begin{aligned} \left(M^{n}+k^{2} \theta^{2} A^{n}\right) U^{n} &=M^{n, n-1} U^{n-1}-k^{2} \theta(1-\theta) A^{n, n-1} U^{n-1}+k M^{n, n-1} V^{n-1}+k^{2} \theta\left[\theta F^{n}+(1-\theta) F^{n-1}\right], \\ M^{n} V^{n} &=M^{n, n-1} V^{n-1}-k\left[\theta A^{n} U^{n}+(1-\theta) A^{n, n-1} U^{n-1}\right]+k\left[\theta F^{n}+(1-\theta) F^{n-1}\right], \end{aligned} \] where \[ \begin{aligned} M_{i j}^{n} &=\left(\phi_{i}^{n}, \phi_{j}^{n}\right), \\ A_{i j}^{n} &=\left(\nabla \phi_{i}^{n}, \nabla \phi_{j}^{n}\right), \\ M_{i j}^{n, n-1} &=\left(\phi_{i}^{n}, \phi_{j}^{n-1}\right), \\ A_{i j}^{n, n-1} &=\left(\nabla \phi_{i}^{n}, \nabla \phi_{j}^{n-1}\right), \\ F_{i}^{n} &=\left(f^{n}, \phi_{i}^{n}\right), \\ F_{i}^{n-1} &=\left(f^{n-1}, \phi_{i}^{n}\right) . \end{aligned} \] If we solve these two equations, we can move the solution one step forward and go on to the next time step.

It is worth noting that if we choose the same mesh on each time step (as we will in fact do in the program below), then we have the same shape functions on time step \(n\) and \(n-1\), i.e. \(\phi_{i}^{n}=\phi_{i}^{n-1}=\phi_{i}\). Consequently, we get \(M^{n}=M^{n, n-1}=M\) and \(A^{n}=A^{n, n-1}=A\). On the other hand, if we had used different shape functions, then we would have to compute integrals that contain shape functions defined on two meshes. This is a somewhat messy process that we omit here, but that is treated in some detail in step-28.

Files output and slicing

A senior graduate student (Chengjiang Yin) in our NJU FEM-NR group has establish a python script (almost a package) to automatically slice the data and save them in .csv files for efficiently ploting..