 |
Reference documentation for deal.II version 9.2.0
|
|
Reference documentation for deal.II version 9.2.0
|
This tutorial depends on step-12.
| Table of contents | |
|---|---|
This program was contributed by Natasha Sharma, Guido Kanschat, Timo Heister, Wolfgang Bangerth, and Zhuoran Wang.
The first author would like to acknowledge the support of NSF Grant No. DMS-1520862. Timo Heister and Wolfgang Bangerth acknowledge support through NSF awards DMS-1821210, EAR-1550901, and OAC-1835673.
This program deals with the biharmonic equation,
\begin{align*} \Delta^2 u(\mathbf x) &= f(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \Omega. \end{align*}
This equation appears in the modeling of thin structures such as roofs of stadiums. These objects are of course in reality three-dimensional with a large aspect ratio of lateral extent to perpendicular thickness, but one can often very accurately model these structures as two dimensional by making assumptions about how internal forces vary in the perpendicular direction. These assumptions lead to the equation above.
The model typically comes in two different kinds, depending on what kinds of boundary conditions are imposed. The first case,
\begin{align*} u(\mathbf x) &= g(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \\ \Delta u(\mathbf x) &= h(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \end{align*}
corresponds to the edges of the thin structure attached to the top of a wall of height \(g(\mathbf x)\) in such a way that the bending forces that act on the structure are \(h(\mathbf x)\); in most physical situations, one will have \(h=0\), corresponding to the structure simply sitting atop the wall.
In the second possible case of boundary values, one would have
\begin{align*} u(\mathbf x) &= g(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \\ \frac{\partial u(\mathbf x)}{\partial \mathbf n} &= j(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega. \end{align*}
This corresponds to a "clamped" structure for which a nonzero \(j(\mathbf x)\) implies a certain angle against the horizontal.
As with Dirichlet and Neumann boundary conditions for the Laplace equation, it is of course possible to have one kind of boundary conditions on one part of the boundary, and the other on the remainder.
The fundamental issue with the equation is that it takes four derivatives of the solution. In the case of the Laplace equation we treated in step-3, step-4, and several other tutorial programs, one multiplies by a test function, integrates, integrates by parts, and ends up with only one derivative on both the test function and trial function – something one can do with functions that are continuous globally, but may have kinks at the interfaces between cells: The derivative may not be defined at the interfaces, but that is on a lower-dimensional manifold (and so doesn't show up in the integrated value).
But for the biharmonic equation, if one followed the same procedure using integrals over the entire domain (i.e., the union of all cells), one would end up with two derivatives on the test functions and trial functions each. If one were to use the usual piecewise polynomial functions with their kinks on cell interfaces, the first derivative would yield a discontinuous gradient, and the second derivative with delta functions on the interfaces – but because both the second derivatives of the test functions and of the trial functions yield a delta function, we would try to integrate the product of two delta functions. For example, in 1d, where \(\varphi_i\) are the usual piecewise linear "hat functions", we would get integrals of the sort
\begin{align*} \int_0^L (\Delta \varphi_i) (\Delta \varphi_j) = \int_0^L \frac 1h \left[\delta(x-x_{i-1}) - 2\delta(x-x_i) + \delta(x-x_{i+1})\right] \frac 1h \left[\delta(x-x_{j-1}) - 2\delta(x-x_j) + \delta(x-x_{j+1})\right] \end{align*}
where \(x_i\) is the node location at which the shape function \(\varphi_i\) is defined, and \(h\) is the mesh size (assumed uniform). The problem is that delta functions in integrals are defined using the relationship
\begin{align*} \int_0^L \delta(x-\hat x) f(x) \; dx = f(\hat x). \end{align*}
But that only works if (i) \(f(\cdot)\) is actually well defined at \(\hat x\), and (ii) if it is finite. On the other hand, an integral of the form
\begin{align*} \int_0^L \delta(x-x_i) \delta (x-x_i) \end{align*}
does not make sense. Similar reasoning can be applied for 2d and 3d situations.
In other words: This approach of trying to integrate over the entire domain and then integrating by parts can't work.
Historically, numerical analysts have tried to address this by inventing finite elements that are "C<sup>1</sup> continuous", i.e., that use shape functions that are not just continuous but also have continuous first derivatives. This is the realm of elements such as the Argyris element, the Clough-Tocher element and others, all developed in the late 1960s. From a twenty-first century perspective, they can only be described as bizarre in their construction. They are also exceedingly cumbersome to implement if one wants to use general meshes. As a consequence, they have largely fallen out of favor and deal.II currently does not contain implementations of these shape functions.
So how does one approach solving such problems then? That depends a bit on the boundary conditions. If one has the first set of boundary conditions, i.e., if the equation is
\begin{align*} \Delta^2 u(\mathbf x) &= f(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \Omega, \\ u(\mathbf x) &= g(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \\ \Delta u(\mathbf x) &= h(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \end{align*}
then the following trick works (at least if the domain is convex, see below): In the same way as we obtained the mixed Laplace equation of step-20 from the regular Laplace equation by introducing a second variable, we can here introduce a variable \(v=\Delta u\) and can then replace the equations above by the following, "mixed" system:
\begin{align*} -\Delta u(\mathbf x) +v(\mathbf x) &= 0 \qquad \qquad &&\forall \mathbf x \in \Omega, \\ -\Delta v(\mathbf x) &= -f(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \Omega, \\ u(\mathbf x) &= g(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \\ v(\mathbf x) &= h(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega. \end{align*}
In other words, we end up with what is in essence a system of two coupled Laplace equations for \(u,v\), each with Dirichlet-type boundary conditions. We know how to solve such problems, and it should not be very difficult to construct good solvers and preconditioners for this system either using the techniques of step-20 and step-22. So this case is pretty simple to deal with.
The more complicated situation is if we have the "clamped" boundary conditions, i.e., if the equation looks like this:
\begin{align*} \Delta^2 u(\mathbf x) &= f(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \Omega, \\ u(\mathbf x) &= g(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \\ \frac{\partial u(\mathbf x)}{\partial \mathbf n} &= j(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega. \end{align*}
The same trick with the mixed system does not work here, because we would end up with both Dirichlet and Neumann boundary conditions for \(u\), but none for \(v\).
The solution to this conundrum arrived with the Discontinuous Galerkin method wave in the 1990s and early 2000s: In much the same way as one can use discontinuous shape functions for the Laplace equation by penalizing the size of the discontinuity to obtain a scheme for an equation that has one derivative on each shape function, we can use a scheme that uses continuous (but not \(C^1\) continuous) shape functions and penalize the jump in the derivative to obtain a scheme for an equation that has two derivatives on each shape function. In analogy to the Interior Penalty (IP) method for the Laplace equation, this scheme for the biharmonic equation is typically called the \(C^0\) IP (or C0IP) method, since it uses \(C^0\) (continuous but not continuously differentiable) shape functions with an interior penalty formulation.
We base this program on the \(C^0\) IP method presented by Susanne Brenner and Li-Yeng Sung in the paper "C \(^0\) Interior Penalty Method for Linear Fourth Order Boundary Value Problems on polygonal domains'' [12] , where the method is derived for the biharmonic equation with "clamped" boundary conditions.
As mentioned, this method relies on the use of \(C^0\) Lagrange finite elements where the \(C^1\) continuity requirement is relaxed and has been replaced with interior penalty techniques. To derive this method, we consider a \(C^0\) shape function \(v_h\) which vanishes on \(\partial\Omega\). We introduce notation \( \mathbb{F} \) as the set of all faces of \(\mathbb{T}\), \( \mathbb{F}^b \) as the set of boundary faces, and \( \mathbb{F}^i \) as the set of interior faces for use further down below. Since the higher order derivatives of \(v_h\) have two values on each interface \(e\in \mathbb{F}\) (shared by the two cells \(K_{+},K_{-} \in \mathbb{T}\)), we cope with this discontinuity by defining the following single-valued functions on \(e\):
\begin{align*} \jump{\frac{\partial^k v_h}{\partial \mathbf n^k}} &= \frac{\partial^k v_h|_{K_+}}{\partial \mathbf n^k} \bigg |_e - \frac{\partial^k v_h|_{K_-}}{\partial \mathbf n^k} \bigg |_e, \\ \average{\frac{\partial^k v_h}{\partial \mathbf n^k}} &= \frac{1}{2} \bigg( \frac{\partial^k v_h|_{K_+}}{\partial \mathbf n^k} \bigg |_e + \frac{\partial^k v_h|_{K_-}}{\partial \mathbf n^k} \bigg |_e \bigg ) \end{align*}
for \(k =1,2\) (i.e., for the gradient and the matrix of second derivatives), and where \(\mathbf n\) denotes a unit vector normal to \(e\) pointing from \(K_+\) to \(K_-\). In the literature, these functions are referred to as the "jump" and "average" operations, respectively.
To obtain the \(C^0\) IP approximation \(u_h\), we left multiply the biharmonic equation by \(v_h\), and then integrate over \(\Omega\). As explained above, we can't do the integration by parts on all of \(\Omega\) with these shape functions, but we can do it on each cell individually since the shape functions are just polynomials on each cell. Consequently, we start by using the following integration-by-parts formula on each mesh cell \(K \in {\mathbb{T}}\):
\begin{align*} \int_K v_h (\Delta^2 w_h) &= \int_K v_h (\nabla\cdot\nabla) (\Delta w_h) \\ &= -\int_K \nabla v_h \cdot (\nabla \Delta w_h) +\int_{\partial K} v_h (\nabla \Delta w_h \cdot \mathbf n). \end{align*}
At this point, we have two options: We can integrate the domain term's \(\nabla\Delta w_h\) one more time to obtain
\begin{align*} \int_K v_h (\Delta^2 w_h) &= \int_K (\Delta v_h) (\Delta w_h) +\int_{\partial K} v_h (\nabla \Delta w_h \cdot \mathbf n) -\int_{\partial K} (\nabla v_h \cdot \mathbf n) \Delta w_h. \end{align*}
For a variety of reasons, this turns out to be a variation that is not useful for our purposes.
Instead, what we do is recognize that \(\nabla\Delta w_h = \text{grad}\,(\text{div}\,\text{grad}\, w_h)\), and we can re-sort these operations as \(\nabla\Delta w_h = \text{div}\,(\text{grad}\,\text{grad}\, w_h)\) where we typically write \(\text{grad}\,\text{grad}\, w_h = D^2 w_h\) to indicate that this is the "Hessian" matrix of second derivatives. With this re-ordering, we can now integrate the divergence, rather than the gradient operator, and we get the following instead:
\begin{align*} \int_K v_h (\Delta^2 w_h) &= \int_K (\nabla \nabla v_h) : (\nabla \nabla w_h) +\int_{\partial K} v_h (\nabla \Delta w_h \cdot \mathbf n) -\int_{\partial K} (\nabla v_h \otimes \mathbf n) : (\nabla\nabla w_h) \\ &= \int_K (D^2 v_h) : (D^2 w_h) +\int_{\partial K} v_h (\nabla \Delta w_h \cdot \mathbf n) -\int_{\partial K} (\nabla v_h) \cdot (D^2 w_h \mathbf n). \end{align*}
Here, the colon indicates a double-contraction over the indices of the matrices to its left and right, i.e., the scalar product between two tensors. The outer product of two vectors \(a \otimes b\) yields the matrix \((a \otimes b)_{ij} = a_i b_j\).
Then, we sum over all cells \(K \in \mathbb{T}\), and take into account that this means that every interior face appears twice in the sum. If we therefore split everything into a sum of integrals over cell interiors and a separate sum over cell interfaces, we can use the jump and average operators defined above. There are two steps left: First, because our shape functions are continuous, the gradients of the shape functions may be discontinuous, but the continuity guarantees that really only the normal component of the gradient is discontinuous across faces whereas the tangential component(s) are continuous. Second, the discrete formulation that results is not stable as the mesh size goes to zero, and to obtain a stable formulation that converges to the correct solution, we need to add the following terms:
\begin{align*} -\sum_{e \in \mathbb{F}} \int_{e} \average{\frac{\partial^2 v_h}{\partial \mathbf n^2}} \jump{\frac{\partial u_h}{\partial \mathbf n}} + \sum_{e \in \mathbb{F}} \frac{\gamma}{h_e}\int_e \jump{\frac{\partial v_h}{\partial \mathbf n}} \jump{\frac{\partial u_h}{\partial \mathbf n}}. \end{align*}
Then, after making cancellations that arise, we arrive at the following C0IP formulation of the biharmonic equation: find \(u_h\) such that \(u_h = g\) on \(\partial \Omega\) and
\begin{align*} \mathcal{A}(v_h,u_h)&=\mathcal{F}(v_h) \quad \text{holds for all test functions } v_h, \end{align*}
where
\begin{align*} \mathcal{A}(v_h,u_h):=&\sum_{K \in \mathbb{T}}\int_K D^2v_h:D^2u_h \ dx \\ & -\sum_{e \in \mathbb{F}} \int_{e} \jump{\frac{\partial v_h}{\partial \mathbf n}} \average{\frac{\partial^2 u_h}{\partial \mathbf n^2}} \ ds -\sum_{e \in \mathbb{F}} \int_{e} \average{\frac{\partial^2 v_h}{\partial \mathbf n^2}} \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds \\ &+ \sum_{e \in \mathbb{F}} \frac{\gamma}{h_e} \int_e \jump{\frac{\partial v_h}{\partial \mathbf n}} \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds, \end{align*}
and
\begin{align*} \mathcal{F}(v_h)&:=\sum_{K \in \mathbb{T}}\int_{K} v_h f \ dx - \sum_{e \in \mathbb{F}, e\subset\partial\Omega} \int_e \average{\frac{\partial^2 v_h}{\partial \mathbf n^2}} j \ ds + \sum_{e \in \mathbb{F}, e\subset\partial\Omega} \frac{\gamma}{h_e} \int_e \jump{\frac{\partial v_h}{\partial \mathbf n}} j \ ds. \end{align*}
Here, \(\gamma\) is the penalty parameter which both weakly enforces the boundary condition
\begin{align*} \frac{\partial u(\mathbf x)}{\partial \mathbf n} = j(\mathbf x) \end{align*}
on the boundary interfaces \(e \in \mathbb{F}^b\), and also ensures that in the limit \(h\rightarrow 0\), \(u_h\) converges to a \(C^1\) continuous function. \(\gamma\) is chosen to be large enough to guarantee the stability of the method. We will discuss our choice in the program below.
On polygonal domains, the weak solution \(u\) to the biharmonic equation lives in \(H^{2 +\alpha}(\Omega)\) where \(\alpha \in(1/2, 2]\) is determined by the interior angles at the corners of \(\Omega\). For instance, whenever \(\Omega\) is convex, \(\alpha=1\); \(\alpha\) may be less than one if the domain has re-entrant corners but \(\alpha\) is close to \(1\) if one of all interior angles is close to \(\pi\).
Now suppose that the \(C^0\) IP solution \(u_h\) is approximated by \(C^0\) shape functions with polynomial degree \(p \ge 2\). Then the discretization outlined above yields the convergence rates as discussed below.
Convergence in the \(C^0\) IP-norm
Ideally, we would like to measure convergence in the "energy norm" \(\|D^2(u-u_h)\|\). However, this does not work because, again, the discrete solution \(u_h\) does not have two (weak) derivatives. Instead, one can define a discrete ( \(C^0\) IP) seminorm that is "equivalent" to the energy norm, as follows:
\begin{align*} |u_h|_{h}^2 := \sum\limits_{K \in \mathbb{T}} \big|u_h\big|_{H^2(K)}^2 + \sum\limits_{e \in \mathbb{F} } \frac{\gamma }{h_e} \left\| \jump{\frac{\partial u_h}{\partial \mathbf n}} \right\|_{L^2(e)}^2. \end{align*}
In this seminorm, the theory in the paper mentioned above yields that we can expect
\begin{align*} |u-u_h|_{h}^2 = {\cal O}(h^{p-1}), \end{align*}
much as one would expect given the convergence rates we know are true for the usual discretizations of the Laplace equation.
Of course, this is true only if the exact solution is sufficiently smooth. Indeed, if \(f \in H^m(\Omega)\) with \(m \ge 0\), \(u \in H^{2+\alpha}(\Omega)\) where \( 2 < 2+\alpha \le m+4\), then the convergence rate of the \(C^0\) IP method is \(\mathcal{O}(h^{\min\{p-1, \alpha\}})\). In other words, the optimal convergence rate can only be expected if the solution is so smooth that \(\alpha\ge p-1\); this can only happen if (i) the domain is convex with a sufficiently smooth boundary, and (ii) \(m\ge p-3\). In practice, of course, the solution is what it is (independent of the polynomial degree we choose), and the last condition can then equivalently be read as saying that there is definitely no point in choosing \(p\) large if \(m\) is not also large. In other words, the only reasonably choices for \(p\) are \(p\le m+3\) because larger polynomial degrees do not result in higher convergence orders.
For the purposes of this program, we're a bit too lazy to actually implement this equivalent seminorm – though it's not very difficult and would make for a good exercise. Instead, we'll simply check in the program what the "broken" \(H^2\) seminorm
\begin{align*} \left(|u_h|^\circ_{H^2}\right)^2 := \sum\limits_{K \in \mathbb{T}} \big|u_h\big|_{H^2(K)}^2 = \sum\limits_{K \in \mathbb{T}} \big|D^2 u_h\big|_{L_2}^2 \end{align*}
yields. The convergence rate in this norm can, from a theoretical perspective, of course not be worse than the one for \(|\cdot|_h\) because it contains only a subset of the necessary terms, but it could at least conceivably be better. It could also be the case that we get the optimal convergence rate even though there is a bug in the program, and that that bug would only show up in sub-optimal rates for the additional terms present in \(|\cdot|_h\). But, one might hope that if we get the optimal rate in the broken norm and the norms discussed below, then the program is indeed correct. The results section will demonstrate that we obtain optimal rates in all norms shown.
Convergence in the \(L_2\)-norm
The optimal convergence rate in the \(L_2\)-norm is \(\mathcal{O}(h^{p+1})\) provided \(p \ge 3\). More details can be found in Theorem 4.6 of [23] .
The default in the program below is to choose \(p=2\). In that case, the theorem does not apply, and indeed one only gets \(\mathcal{O}(h^2)\) instead of \(\mathcal{O}(h^3)\) as we will show in the results section.
Convergence in the \(H^1\)-seminorm
Given that we expect \(\mathcal{O}(h^{p-1})\) in the best of cases for a norm equivalent to the \(H^2\) seminorm, and \(\mathcal{O}(h^{p+1})\) for the \(L_2\) norm, one may ask about what happens in the \(H^1\) seminorm that is intermediate to the two others. A reasonable guess is that one should expect \(\mathcal{O}(h^{p})\). There is probably a paper somewhere that proves this, but we also verify that this conjecture is experimentally true below.
We remark that the derivation of the \(C^0\) IP method for the biharmonic equation with other boundary conditions – for instance, for the first set of boundary conditions namely \(u(\mathbf x) = g(\mathbf x)\) and \(\Delta u(\mathbf x)= h(\mathbf x)\) on \(\partial\Omega\) – can be obtained with suitable modifications to \(\mathcal{A}(\cdot,\cdot)\) and \(\mathcal{F}(\cdot)\) described in the book chapter [14] .
The last step that remains to describe is what this program solves for. As always, a trigonometric function is both a good and a bad choice because it does not lie in any polynomial space in which we may seek the solution while at the same time being smoother than real solutions typically are (here, it is in \(C^\infty\) while real solutions are typically only in \(H^3\) or so on convex polygonal domains, or somewhere between \(H^2\) and \(H^3\) if the domain is not convex). But, since we don't have the means to describe solutions of realistic problems in terms of relatively simple formulas, we just go with the following, on the unit square for the domain \(\Omega\):
\begin{align*} u = \sin(\pi x) \sin(\pi y). \end{align*}
As a consequence, we then need choose as boundary conditions the following:
\begin{align*} g &= u|_{\partial\Omega} = \sin(\pi x) \sin(\pi y)|_{\partial\Omega}, \\ j &= \frac{\partial u}{\partial\mathbf n}|_{\partial\Omega} \\ &= \left.\begin{pmatrix} \pi\cos(\pi x) \sin(\pi y) \\ \pi\sin(\pi x) \cos(\pi y) \end{pmatrix}\right|_{\partial\Omega} \cdot \mathbf n. \end{align*}
The right hand side is easily computed as
\begin{align*} f = \Delta^2 u = 4 \pi^4 \sin(\pi x) \sin(\pi y). \end{align*}
The program has classes ExactSolution::Solution and ExactSolution::RightHandSide that encode this information.
The first few include files have already been used in the previous example, so we will not explain their meaning here again. The principal structure of the program is very similar to that of, for example, step-4 and so we include many of the same header files.
The two most interesting header files will be these two:
The first of these is responsible for providing the class FEInterfaceValues that can be used to evaluate quantities such as the jump or average of shape functions (or their gradients) across interfaces between cells. This class will be quite useful in evaluating the penalty terms that appear in the C0IP formulation.
In the following namespace, let us define the exact solution against which we will compare the numerically computed one. It has the form \(u(x,y) = \sin(\pi x) \sin(\pi y)\) (only the 2d case is implemented), and the namespace also contains a class that corresponds to the right hand side that produces this solution.
The following is the principal class of this tutorial program. It has the structure of many of the other tutorial programs and there should really be nothing particularly surprising about its contents or the constructor that follows it.
Next up are the functions that create the initial mesh (a once refined unit square) and set up the constraints, vectors, and matrices on each mesh. Again, both of these are essentially unchanged from many previous tutorial programs.
The following pieces of code are more interesting. They all relate to the assembly of the linear system. While assembling the cell-interior terms is not of great difficulty – that works in essence like the assembly of the corresponding terms of the Laplace equation, and you have seen how this works in step-4 or step-6, for example – the difficulty is with the penalty terms in the formulation. These require the evaluation of gradients of shape functions at interfaces of cells. At the least, one would therefore need to use two FEFaceValues objects, but if one of the two sides is adaptively refined, then one actually needs an FEFaceValues and one FESubfaceValues objects; one also needs to keep track which shape functions live where, and finally we need to ensure that every face is visited only once. All of this is a substantial overhead to the logic we really want to implement (namely the penalty terms in the bilinear form). As a consequence, we will make use of the FEInterfaceValues class – a helper class in deal.II that allows us to abstract away the two FEFaceValues or FESubfaceValues objects and directly access what we really care about: jumps, averages, etc.
But this doesn't yet solve our problem of having to keep track of which faces we have already visited when we loop over all cells and all of their faces. To make this process simpler, we use the MeshWorker::mesh_loop() function that provides a simple interface for this task: Based on the ideas outlined in the WorkStream namespace documentation, MeshWorker::mesh_loop() requires three functions that do work on cells, interior faces, and boundary faces. These functions work on scratch objects for intermediate results, and then copy the result of their computations into copy data objects from where a copier function copies them into the global matrix and right hand side objects.
The following structures then provide the scratch and copy objects that are necessary for this approach. You may look up the WorkStream namespace as well as the Parallel computing with multiple processors module for more information on how they typically work.
The more interesting part is where we actually assemble the linear system. Fundamentally, this function has five parts:
cell_worker lambda function, a small function that is defined within the assemble_system() function and that will be responsible for computing the local integrals on an individual cell. It will work on a copy of the ScratchData class and put its results into the corresponding CopyData object.face_worker lambda function that does the integration of all terms that live on the interfaces between cells.boundary_worker function that does the same but for cell faces located on the boundary of the domain.copier function that is responsible for copying all of the data the previous three functions have put into copy objects for a single cell, into the global matrix and right hand side.The fifth part is the one where we bring all of this together.
Let us go through each of these pieces necessary for the assembly in turns.
The first piece is the cell_worker that does the assembly on the cell interiors. It is a (lambda) function that takes a cell (input), a scratch object, and a copy object (output) as arguments. It looks like the assembly functions of many other of the tutorial programs, or at least the body of the loop over all cells.
The terms we integrate here are the cell contribution
\begin{align*} A^K_{ij} = \int_K \nabla^2\varphi_i(x) : \nabla^2\varphi_j(x) dx \end{align*}
to the global matrix, and
\begin{align*} f^K_i = \int_K \varphi_i(x) f(x) dx \end{align*}
to the right hand side vector.
We use the same technique as used in the assembly of step-22 to accelerate the function: Instead of calling fe_values.shape_hessian(i, qpoint) in the innermost loop, we instead create a variable hessian_i that evaluates this value once in the loop over i and re-use the so-evaluated value in the loop over j. For symmetry, we do the same with a variable hessian_j, although it is indeed only used once and we could have left the call to fe_values.shape_hessian(j,qpoint) in the instruction that computes the scalar product between the two terms.
The next building block is the one that assembles penalty terms on each of the interior faces of the mesh. As described in the documentation of MeshWorker::mesh_loop(), this function receives arguments that denote a cell and its neighboring cell, as well as (for each of the two cells) the face (and potentially sub-face) we have to integrate over. Again, we also get a scratch object, and a copy object for putting the results in.
The function has three parts itself. At the top, we initialize the FEInterfaceValues object and create a new CopyData::FaceData object to store our input in. This gets pushed to the end of the copy_data.face_data variable. We need to do this because the number of faces (or subfaces) over which we integrate for a given cell differs from cell to cell, and the sizes of these matrices also differ, depending on what degrees of freedom are adjacent to the face or subface. As discussed in the documentation of MeshWorker::mesh_loop(), the copy object is reset every time a new cell is visited, so that what we push to the end of copy_data.face_data() is really all that the later copier function gets to see when it copies the contributions of each cell to the global matrix and right hand side objects.
The second part deals with determining what the penalty parameter should be. By looking at the units of the various terms in the bilinear form, it is clear that the penalty has to have the form \(\frac{\gamma}{h_K}\) (i.e., one over length scale), but it is not a priori obvious how one should choose the dimension-less number \(\gamma\). From the discontinuous Galerkin theory for the Laplace equation, one might conjecture that the right choice is \(\gamma=p(p+1)\) is the right choice, where \(p\) is the polynomial degree of the finite element used. We will discuss this choice in a bit more detail in the results section of this program.
In the formula above, \(h_K\) is the size of cell \(K\). But this is not quite so straightforward either: If one uses highly stretched cells, then a more involved theory says that \(h\) should be replaced by the diameter of cell \(K\) normal to the direction of the edge in question. It turns out that there is a function in deal.II for that. Secondly, \(h_K\) may be different when viewed from the two different sides of a face.
To stay on the safe side, we take the maximum of the two values. We will note that it is possible that this computation has to be further adjusted if one were to use hanging nodes resulting from adaptive mesh refinement.
Finally, and as usual, we loop over the quadrature points and indices i and j to add up the contributions of this face or sub-face. These are then stored in the copy_data.face_data object created above. As for the cell worker, we pull the evaluation of averages and jumps out of the loops if possible, introducing local variables that store these results. The assembly then only needs to use these local variables in the innermost loop. Regarding the concrete formula this code implements, recall that the interface terms of the bilinear form were as follows:
\begin{align*} -\sum_{e \in \mathbb{F}} \int_{e} \jump{ \frac{\partial v_h}{\partial \mathbf n}} \average{\frac{\partial^2 u_h}{\partial \mathbf n^2}} \ ds -\sum_{e \in \mathbb{F}} \int_{e} \average{\frac{\partial^2 v_h}{\partial \mathbf n^2}} \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds + \sum_{e \in \mathbb{F}} \frac{\gamma}{h_e} \int_e \jump{\frac{\partial v_h}{\partial \mathbf n}} \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds. \end{align*}
The third piece is to do the same kind of assembly for faces that are at the boundary. The idea is the same as above, of course, with only the difference that there are now penalty terms that also go into the right hand side.
As before, the first part of the function simply sets up some helper objects:
Positively, because we now only deal with one cell adjacent to the face (as we are on the boundary), the computation of the penalty factor \(\gamma\) is substantially simpler:
The third piece is the assembly of terms. This is now slightly more involved since these contains both terms for the matrix and for the right hand side. The former is exactly the same as for the interior faces stated above if one just defines the jump and average appropriately (which is what the FEInterfaceValues class does). The latter requires us to evaluate the boundary conditions \(j(\mathbf x)\), which in the current case (where we know the exact solution) we compute from \(j(\mathbf x) = \frac{\partial u(\mathbf x)}{\partial {\mathbf n}}\). The term to be added to the right hand side vector is then \(\frac{\gamma}{h_e}\int_e \jump{\frac{\partial v_h}{\partial \mathbf n}} j \ ds\).
Part 4 is a small function that copies the data produced by the cell, interior, and boundary face assemblers above into the global matrix and right hand side vector. There really is not very much to do here: We distribute the cell matrix and right hand side contributions as we have done in almost all of the other tutorial programs using the constraints objects. We then also have to do the same for the face matrix contributions that have gained content for the faces (interior and boundary) and that the face_worker and boundary_worker have added to the copy_data.face_data array.
Having set all of this up, what remains is to just create a scratch and copy data object and call the MeshWorker::mesh_loop() function that then goes over all cells and faces, calls the respective workers on them, and then the copier function that puts things into the global matrix and right hand side. As an additional benefit, MeshWorker::mesh_loop() does all of this in parallel, using as many processor cores as your machine happens to have.
The show is essentially over at this point: The remaining functions are not overly interesting or novel. The first one simply uses a direct solver to solve the linear system (see also step-29):
The next function evaluates the error between the computed solution and the exact solution (which is known here because we have chosen the right hand side and boundary values in a way so that we know the corresponding solution). In the first two code blocks below, we compute the error in the \(L_2\) norm and the \(H^1\) semi-norm.
Now also compute an approximation to the \(H^2\) seminorm error. The actual \(H^2\) seminorm would require us to integrate second derivatives of the solution \(u_h\), but given the Lagrange shape functions we use, \(u_h\) of course has kinks at the interfaces between cells, and consequently second derivatives are singular at interfaces. As a consequence, we really only integrate over the interior of cells and ignore the interface contributions. This is not an equivalent norm to the energy norm for the problem, but still gives us an idea of how fast the error converges.
We note that one could address this issue by defining a norm that is equivalent to the energy norm. This would involve adding up not only the integrals over cell interiors as we do below, but also adding penalty terms for the jump of the derivative of \(u_h\) across interfaces, with an appropriate scaling of the two kinds of terms. We will leave this for later work.
Equally uninteresting is the function that generates graphical output. It looks exactly like the one in step-6, for example.
The same is true for the run() function: Just like in previous programs.
Finally for the main() function. There is, again, not very much to see here: It looks like the ones in previous tutorial programs. There is a variable that allows selecting the polynomial degree of the element we want to use for solving the equation. Because the C0IP formulation we use requires the element degree to be at least two, we check with an assertion that whatever one sets for the polynomial degree actually makes sense.
We run the program with right hand side and boundary values as discussed in the introduction. These will produce the solution \(u = \sin(\pi x) \sin(\pi y)\) on the domain \(\Omega = (0,1)^2\). We test this setup using \(Q_2\), \(Q_3\), and \(Q_4\) elements, which one can change via the fe_degree variable in the main() function. With mesh refinement, the \(L_2\) convergence rates, \(H^1\)-seminorm rate, and \(H^2\)-seminorm convergence of \(u\) should then be around 2, 2, 1 for \(Q_2\) (with the \(L_2\) norm sub-optimal as discussed in the introduction); 4, 3, 2 for \(Q_3\); and 5, 4, 3 for \(Q_4\), respectively.
From the literature, it is not immediately clear what the penalty parameter \(\gamma\) should be. For example, [13] state that it needs to be larger than one, and choose \(\gamma=5\). The FEniCS/Dolphin tutorial chooses it as \(\gamma=8\), see https://fenicsproject.org/docs/dolfin/1.6.0/python/demo/documented/biharmonic/python/documentation.html . [67] uses a value for \(\gamma\) larger than the number of edges belonging to an element for Kirchhoff plates (see their Section 4.2). This suggests that maybe \(\gamma = 1\), \(2\), are too small; on the other hand, a value \(p(p+1)\) would be reasonable, where \(p\) is the degree of polynomials. The last of these choices is the one one would expect to work by comparing to the discontinuous Galerkin formulations for the Laplace equation, and it will turn out to also work here. But we should check what value of \(\gamma\) is right, and we will do so below; changing \(\gamma\) is easy in the two face_worker and boundary_worker functions defined in assemble_system().
We run the code with differently refined meshes and get the following convergence rates.
| Number of refinements | \(\|u-u_h^\circ\|_{L_2}\) | Conv. rates | \(|u-u_h|_{H^1}\) | Conv. rates | \(|u-u_h|_{H^2}\) | Conv. rates |
|---|---|---|---|---|---|---|
| 2 | 8.780e-03 | 7.095e-02 | 1.645 | |||
| 3 | 3.515e-03 | 1.32 | 2.174e-02 | 1.70 | 8.121e-01 | 1.018 |
| 4 | 1.103e-03 | 1.67 | 6.106e-03 | 1.83 | 4.015e-01 | 1.016 |
| 5 | 3.084e-04 | 1.83 | 1.622e-03 | 1.91 | 1.993e-01 | 1.010 |
We can see that the \(L_2\) convergence rates are around 2, \(H^1\)-seminorm convergence rates are around 2, and \(H^2\)-seminorm convergence rates are around 1. The latter two match the theoretically expected rates; for the former, we have no theorem but are not surprised that it is sub-optimal given the remark in the introduction.
| Number of refinements | \(\|u-u_h^\circ\|_{L_2}\) | Conv. rates | \(|u-u_h|_{H^1}\) | Conv. rates | \(|u-u_h|_{H^2}\) | Conv. rates |
|---|---|---|---|---|---|---|
| 2 | 2.045e-04 | 4.402e-03 | 1.641e-01 | |||
| 3 | 1.312e-05 | 3.96 | 5.537e-04 | 2.99 | 4.096e-02 | 2.00 |
| 4 | 8.239e-07 | 3.99 | 6.904e-05 | 3.00 | 1.023e-02 | 2.00 |
| 5 | 5.158e-08 | 3.99 | 8.621e-06 | 3.00 | 2.558e-03 | 2.00 |
We can see that the \(L_2\) convergence rates are around 4, \(H^1\)-seminorm convergence rates are around 3, and \(H^2\)-seminorm convergence rates are around 2. This, of course, matches our theoretical expectations.
| Number of refinements | \(\|u-u_h^\circ\|_{L_2}\) | Conv. rates | \(|u-u_h|_{H^1}\) | Conv. rates | \(|u-u_h|_{H^2}\) | Conv. rates |
|---|---|---|---|---|---|---|
| 2 | 6.510e-06 | 2.215e-04 | 1.275e-02 | |||
| 3 | 2.679e-07 | 4.60 | 1.569e-05 | 3.81 | 1.496e-03 | 3.09 |
| 4 | 9.404e-09 | 4.83 | 1.040e-06 | 3.91 | 1.774e-04 | 3.07 |
| 5 | 7.943e-10 | 3.56 | 6.693e-08 | 3.95 | 2.150e-05 | 3.04 |
We can see that the \(L_2\) norm convergence rates are around 5, \(H^1\)-seminorm convergence rates are around 4, and \(H^2\)-seminorm convergence rates are around 3. On the finest mesh, the \(L_2\) norm convergence rate is much smaller than our theoretical expectations because the linear solver becomes the limiting factor due to round-off. Of course the \(L_2\) error is also very small already in that case.
For comparison with the results above, let us now also consider the case where we simply choose \(\gamma=1\):
| Number of refinements | \(\|u-u_h^\circ\|_{L_2}\) | Conv. rates | \(|u-u_h|_{H^1}\) | Conv. rates | \(|u-u_h|_{H^2}\) | Conv. rates |
|---|---|---|---|---|---|---|
| 2 | 7.350e-02 | 7.323e-01 | 10.343 | |||
| 3 | 6.798e-03 | 3.43 | 1.716e-01 | 2.09 | 4.836 | 1.09 |
| 4 | 9.669e-04 | 2.81 | 6.436e-02 | 1.41 | 3.590 | 0.430 |
| 5 | 1.755e-04 | 2.46 | 2.831e-02 | 1.18 | 3.144 | 0.19 |
Although \(L_2\) norm convergence rates of \(u\) more or less follows the theoretical expectations, the \(H^1\)-seminorm and \(H^2\)-seminorm do not seem to converge as expected. Comparing results from \(\gamma = 1\) and \(\gamma = p(p+1)\), it is clear that \(\gamma = p(p+1)\) is a better penalty. Given that \(\gamma=1\) is already too small for \(Q_2\) elements, it may not be surprising that if one repeated the experiment with a \(Q_3\) element, the results are even more disappointing: One again only obtains convergence rates of 2, 1, zero – i.e., no better than for the \(Q_2\) element (although the errors are smaller in magnitude). Maybe surprisingly, however, one obtains more or less the expected convergence orders when using \(Q_4\) elements. Regardless, this uncertainty suggests that \(\gamma=1\) is at best a risky choice, and at worst an unreliable one and that we should choose \(\gamma\) larger.
Since \(\gamma=1\) is clearly too small, one might conjecture that \(\gamma=2\) might actually work better. Here is what one obtains in that case:
| Number of refinements | \(\|u-u_h^\circ\|_{L_2}\) | Conv. rates | \(|u-u_h|_{H^1}\) | Conv. rates | \(|u-u_h|_{H^2}\) | Conv. rates |
|---|---|---|---|---|---|---|
| 2 | 4.133e-02 | 2.517e-01 | 3.056 | |||
| 3 | 6.500e-03 | 2.66 | 5.916e-02 | 2.08 | 1.444 | 1.08 |
| 4 | 6.780e-04 | 3.26 | 1.203e-02 | 2.296 | 6.151e-01 | 1.231 |
| 5 | 1.622e-04 | 2.06 | 2.448e-03 | 2.297 | 2.618e-01 | 1.232 |
In this case, the convergence rates more or less follow the theoretical expectations, but, compared to the results from \(\gamma = p(p+1)\), are more variable. Again, we could repeat this kind of experiment for \(Q_3\) and \(Q_4\) elements. In both cases, we will find that we obtain roughly the expected convergence rates. Of more interest may then be to compare the absolute size of the errors. While in the table above, for the \(Q_2\) case, the errors on the finest grid are comparable between the \(\gamma=p(p+1)\) and \(\gamma=2\) case, for \(Q_3\) the errors are substantially larger for \(\gamma=2\) than for \(\gamma=p(p+1)\). The same is true for the \(Q_4\) case.
The conclusions for which of the "reasonable" choices one should use for the penalty parameter is that \(\gamma=p(p+1)\) yields the expected results. It is, consequently, what the code uses as currently written.
There are a number of obvious extensions to this program that would make sense:
Similar to the "clamped" boundary condition addressed in the implementation, we will derive the \(C^0\) IP finite element scheme for simply supported plates:
\begin{align*} \Delta^2 u(\mathbf x) &= f(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \Omega, u(\mathbf x) &= g(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \\ \Delta u(\mathbf x) &= h(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega. \end{align*}
We multiply the biharmonic equation by the test function \(v_h\) and integrate over \( K \) and get:
\begin{align*} \int_K v_h (\Delta^2 u_h) &= \int_K (D^2 v_h) : (D^2 u_h) + \int_{\partial K} v_h \frac{\partial (\Delta u_h)}{\partial \mathbf{n}} -\int_{\partial K} (\nabla v_h) \cdot (\frac{\partial \nabla u_h}{\partial \mathbf{n}}). \end{align*}
Summing up over all cells \(K \in \mathbb{T}\),since normal directions of \(\Delta u_h\) are pointing at opposite directions on each interior edge shared by two cells and \(v_h = 0\) on \(\partial \Omega\),
\begin{align*} \sum_{K \in \mathbb{T}} \int_{\partial K} v_h \frac{\partial (\Delta u_h)}{\partial \mathbf{n}} = 0, \end{align*}
and by the definition of jump over cell interfaces,
\begin{align*} -\sum_{K \in \mathbb{T}} \int_{\partial K} (\nabla v_h) \cdot (\frac{\partial \nabla u_h}{\partial \mathbf{n}}) = -\sum_{e \in \mathbb{F}} \int_{e} \jump{\frac{\partial v_h}{\partial \mathbf{n}}} (\frac{\partial^2 u_h}{\partial \mathbf{n^2}}). \end{align*}
We separate interior faces and boundary faces of the domain,
\begin{align*} -\sum_{K \in \mathbb{T}} \int_{\partial K} (\nabla v_h) \cdot (\frac{\partial \nabla u_h}{\partial \mathbf{n}}) = -\sum_{e \in \mathbb{F}^i} \int_{e} \jump{\frac{\partial v_h}{\partial \mathbf{n}}} (\frac{\partial^2 u_h}{\partial \mathbf{n^2}}) - \sum_{e \in \partial \Omega} \int_{e} \jump{\frac{\partial v_h}{\partial \mathbf{n}}} h, \end{align*}
where \(\mathbb{F}^i\) is the set of interior faces. This leads us to
\begin{align*} \sum_{K \in \mathbb{T}} \int_K (D^2 v_h) : (D^2 u_h) \ dx - \sum_{e \in \mathbb{F}^i} \int_{e} \jump{\frac{\partial v_h}{\partial \mathbf{n}}} (\frac{\partial^2 u_h}{\partial \mathbf{n^2}}) \ ds = \sum_{K \in \mathbb{T}}\int_{K} v_h f \ dx + \sum_{e\subset\partial\Omega} \int_{e} \jump{\frac{\partial v_h}{\partial \mathbf{n}}} h \ ds. \end{align*}
In order to symmetrize and stabilize the discrete problem, we add symmetrization and stabilization term. We finally get the \(C^0\) IP finite element scheme for the biharmonic equation: find \(u_h\) such that \(u_h =g\) on \(\partial \Omega\) and
\begin{align*} \mathcal{A}(v_h,u_h)&=\mathcal{F}(v_h) \quad \text{holds for all test functions } v_h, \end{align*}
where
\begin{align*} \mathcal{A}(v_h,u_h)&=\mathcal{F}(v_h) \quad \text{holds for all test functions } v_h, \end{align*}
where
\begin{align*} \mathcal{A}(v_h,u_h):=&\sum_{K \in \mathbb{T}}\int_K D^2v_h:D^2u_h \ dx \\ & -\sum_{e \in \mathbb{F}^i} \int_{e} \jump{\frac{\partial v_h}{\partial \mathbf n}} \average{\frac{\partial^2 u_h}{\partial \mathbf n^2}} \ ds -\sum_{e \in \mathbb{F}^i} \int_{e} \average{\frac{\partial^2 v_h}{\partial \mathbf n^2}} \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds \\ &+ \sum_{e \in \mathbb{F}^i} \frac{\gamma}{h_e} \int_e \jump{\frac{\partial v_h}{\partial \mathbf n}} \jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds, \end{align*}
and
\begin{align*} \mathcal{F}(v_h)&:=\sum_{K \in \mathbb{T}}\int_{K} v_h f \ dx + \sum_{e\subset\partial\Omega} \int_e \jump{\frac{\partial v_h}{\partial \mathbf n}} h \ ds. \end{align*}
The implementation of this boundary case is similar to "clamped" version except for boundary_worker is no longer needed for system assembling and the right hand side is changed according to the formulation.
1.8.17