A Discontinuous Galerkin Method by Patch Reconstruction for Biharmonic Problem

We propose a new discontinuous Galerkin method based on the least-squares patch reconstruction for the biharmonic problem. We prove the optimal error estimate of the proposed method. The two-dimensional and three-dimensional numerical examples are presented to confirm the accuracy and efficiency of the method with several boundary conditions and several types of polygon meshes and polyhedral meshes.


Introduction
The biharmonic boundary value problem is a fourth-order elliptic problem that models the thin plate bending problem in continuum mechanics, and describes slow flows of viscous incompressible fluids. Finite element methods have been employed to approximate this problem from its initial stage and by now there are many successful finite element methods for this problem [1].
Recently, the discontinuous Galerkin (DG) method [2,3,4,5,6,7] has been developed to solve the biharmonic problem. The DG method employs discontinuous basis functions that render great flexibility in the mesh partition and also provide a suitable framework for p-adaptivity. Higher order polynomials are easily implemented in DG method, which may efficiently capture the smooth solutions. To achieve higher accuracy, DG method requires a large number of degrees of freedom on a single element, which gives an extremely large linear system [8,9]. As a compromise of the standard FEM and the DG method, Brenner et al [10,11] developed the socalled C 0 interior penalty Galerkin method (C 0 IPG). This method applies standard continuous C 0 Lagrange finite elements to the interior penalty Galerkin variational formulation, which admits the optimal error estimate with less local degree of freedoms compared with standard DG method.
The aim of this work is to apply the patch reconstruction finite element method proposed in [12] to the biharmonic problem. The main idea of the proposed method is to construct a piecewise polynomial approximation space by patch reconstruction. The approximation space is discontinuous across the element face and has only one degree of freedom located inside each element, which is a sub-space of the commonly used discontinuous Galerkin finite element space. One advantage of the proposed method is the number of the unknown is maintained under a given mesh partition with the increasing of the order of accuracy. Moreover, the reconstruction procedure can be carried out over any mesh such as an arbitrary polygonal mesh. Given such reconstructed approximation space, we solve the biharmonic problem in the framework of interior penalty discontinuous Galerkin (IPDG). Based on the approximation properties of the approximation space established in [13] and [12], we analyze the proposed method in the framework of IPDG. The performance of the proposed method is verified by a series of numerical tests with different complexity, which is comparable with the C 0 IPG method while the basis functions in our method should be pre-computed, nevertheless, such basis functions may be reused.
The article is organized as follows. In Section 2, we demonstrate the reconstruction procedure of the approximation space and develop the corresponding approximation properties of such a space. Next, the IPDG method with such a reconstructed approximation space is proposed and analyzed in Section 3. In Section 4, we test the proposed method by several two-dimensional and three-dimensional biharmonic boundary value problems, which include smooth solution as well as nonsmooth solution for various boundary conditions and different types of meshes.
Throughout this paper, the constant C is a generic constant that may change from line to line, but is independent of the mesh size h.

Reconstruction Operator
Let Ω ⊂ R d with d = 2, 3 be a bounded convex domain. Let T h be a collection of N e polygonal elements that partition Ω. We denote all interior faces of T h as E i h and the set of the boundary faces as E b h , and then We assume that T h satisfies the shape-regular conditions used in [14,15], which read: there exist -an integer number N independent of h; -a real positive number σ independent of h; -a compatible sub-decomposition T h into shape-regular triangles or quadrilaterals, or mix of both triangles and quadrilaterals, such that (A1) any polygon K ∈ T h admits a decomposition T h formed by less than N triangles; (A2) any triangle T ∈ T h is shape-regular in the sense that there exists σ satisfying h K /ρ K ≤ σ, where ρ K is the radius of the largest ball inscribed in K. The above regularity assumptions lead to some useful consequences, which will be extensively used in the later analysis.
M1 For any triangle T ∈ T h , there exists ρ 1 ≥ 1 that depends on N and σ such that h K /h T ≤ ρ 1 . M2 [Agmon inequality] There exists C that depends on N and σ, but independent of h K such that M3 [Approximation property] There exists C that depends on N and σ, but independent of h K such that for any v ∈ H m+1 (K), there exists an approximating polynomial v m ∈ P m (K) such that M4 [Inverse inequality] For any v ∈ P m (K), there exists a constant C that depends only on N and σ such that Given the triangulation T h , we define the reconstruction operator as follows. First, for each element K, we prescribe a point x K as the collocation point. In particular, we may specify x K as the barycenter of K, although it could be more flexible. Second, we construct an element patch S(K) that consists of K itself and some elements around K. The element patch can be built in two ways. The first way is that we initialize S(K) with K, and we add all Von Neumann neighbours of the elements into S(K) in a recursive manner until we have collected  . Another way to construct S(K) has been reported in [13]. Denote by I K the set of the collocation points that belong to S(K). It is clear that #I K = #S(K). Let U h be the piecewise constant space associated with T h , i.e., where P r is the polynomial space of degree not greater than r. For any function g ∈ U h , we reconstruct a polynomial R K g of degree m on S(K) by solving the following least-squares: The uniqueness condition for Problem (2.4) relates to the location of the collocation points and #S(K). Following [13,12], we make the following assumption: Assumption 2.1. For all K ∈ T h and p ∈ P m (S(K)), (2.5) p| I K = 0 implies p| S(K) ≡ 0, The above assumption guarantees the uniqueness of the solution of Problem (2.4) if #S(K) is greater than dim P m . Hereafter, we assume that this assumption is always valid.
From (2.4) we define the global reconstruction operator R for g ∈ U h by restricting the polynomial R K g on K, Rg| K = (R K g)| K . Therefore, the operator R defines a map from U h to a discontinuous piecewise polynomial space with order m, which is denoted as V h = RU h .
Next we investigate the behaviour of the basis functions, which are generated by the reconstruction process. Let e K ∈ R N e , whose components are given by Then we denote {λ K |λ K = Re K for all K ∈ T h } as a group of basis functions. Given {λ K } we may write the reconstruction operator in an explicit way: The operator R may be defined for continuous functions posed on Ω as (2.6) . The support of λ K can be described as: A one-dimensional example is presented in Section 4 and Figure 4.2.
It is worthwhile to mention that firstly the support of the basis function may not be equal to the element patch, and vice versus; secondly, the basis function is discontinuous, which lends itself ideally to the DG framework.
To state the method, we introduce some notations. Define the broken Sobolev spaces: where H s (K) is the standard Sobolev space on element K associated with the norm · H s (K) and the seminorm |·| H s (K) . The associated broken norm and seminorm are defined respectively by For two neighbouring elements K + , K − that share a common face e = ∂K + ∩ ∂K − with e ∈ E i h , denote by n + and n − the outward unit normal vectors on e, corresponding to ∂K + and ∂K − , respectively. For any function q ∈ H 1 (Ω, T h ) and v ∈ [H 1 (Ω, T h )] d that may be discontinuous across interelement boundaries, we define the average operator {·} and the jump operator Here In the case e ∈ E b h , there exists an element K such that e = K ∩ ∂Ω, the definitions are changed to Following [4], for any w ∈ H 2 (Ω, T h ), we define the energy norm · h as This energy norm w h is equivalent to w H 2 (Ω,T h ) . Next, we turn to the stability estimate of the reconstruction operator R. By [13], we define Λ(m, I K ) for any element K ∈ T h as By Assumption 2.1, we may conclude that Λ(m, I K ) is finite. With some further conditions on S(K), we conclude that Λ(m, I K ) has a uniform upper bound, which is denoted by Λ m . One of such condition may be found in the following lemma. Lemma 2.1. [13, Lemma 3.5] If each element patch S(K) is convex and the triangulation is quasi-uniform, then Λ m is uniformly bounded.
The above condition is not necessary, we refer to [12] for the discussion on the uniform upper bound for non-convex element patch.
With the uniformly bounded Λ m , we have the following quasi-optimal approximation estimates for the reconstruction operator in the maximum norm.
Lemma 2.2. [13, Theorem 3.3] If Assumption 2.1 holds, the stability estimate holds true for any K ∈ T h and g ∈ C 0 (S(K)) as where C is independent of h but depends on #S(K).
The above lemma immediately implies the quasi-optimal approximation results in other norms.
Proof. By the standard interpolation estimate and (2.8), we obtain which gives (2.9) with q = 0. Let g m be the approximation polynomial in (2.2) for g, using (2.3), we obtain , which yields (2.9) with q = 1 by using (2.9) with q = 0 in the last step.
Using the above lemma, we obtain the following interpolation estimate for the reconstruction operator in the energy norm.

Error Estimate for the Biharmonic Problem
Let us consider the biharmonic problem: find u ∈ H 4 (Ω), such that where ∆ 2 u = ∆(∆u), n denotes the unit outward normal vector to ∂Ω, f ∈ L 2 (Ω), g D and g N are assumed to be suitably smooth such that under proper conditions on Ω, the boundary value problem (3.1) has a unique solution. We refer to [16] for precise description of such results.
The IPDG method [6] employs the following bilinear form B: for any v, w ∈ H 4 (Ω, T h ), The piecewise penalty parameters α and β are nonnegative and will be specified later on. The linear form l is defined for all v ∈ H 4 (Ω, T h ) as The discretized problem is to find u h ∈ U h such that for all v h ∈ U h , We begin the error estimate by defining the lifting operator L : H 4 (Ω, T h ) → P h [6]: where P h is defined by The following lemma [6] gives the stability of the lifting operator.
Lemma 3.1. For all w ∈ H 4 (Ω, T h ), the following estimate holds: , for piecewise constants γ, δ that are defined for all e ∈ E h by with sufficiently large positive constants C γ and C δ .
With the definition of the DG energy norm and the Lemma 3.1, we show the bilinear form B(·, ·) satisfies the boundedness and the stability condition.
We are ready to prove a priori error estimates for Problem (3.2).
Proof. We begin with the Galerkin orthogonality: for all v h ∈ U h , there holds

Using (3.7) and (3.4), we obtain
By the triangle inequality we immediately obtain which together with (2.11) gives the estimate (3.8).
The L 2 error estimate may be obtained by standard duality argument. Following [5], we make the regularity assumptions: the solution ψ of problem in Ω, belongs to H 4 (Ω), and there exists a positive constant C that only depends on Ω such that (3.11) ψ H 4 (Ω) ≤ C u − Ru h L 2 (Ω) .
Proof. An integration by parts gives where we have used the Galerkin orthogonality in the last step. Combining the interpolation estimate, the energy estimate (3.8) and the regularity estimate (3.11) gives (3.12) and (3.13).

Implementation and Numerical Results
In this section, we describe the implementation details of the reconstructed space and report some numerical examples to show the accuracy and performance of the proposed method. In all examples, we build element patches by the first method and use a sparse direct solver for the resulting sparse linear systems.
The local least-squares (2.4) on element K i is where A and q are given as follows. For i = 1, and for i = 2, 3, 4, and for i = 5, The basis functions {λ K } are plot in Figure 4.2. Thus we store (A T A) −1 A T for each element to represent the basis functions and it is the same when we deal with the high dimensional problem.

2D Smooth Solutions.
We firstly study the convergence rate for smooth solutions of two-dimensional problems.   u(x, y) = sin 2 (πx) sin 2 (πy), (x, y) ∈ Ω, and g D , g N and the source term f are chosen accordingly. The polynomial degree m is taken by m = 2, · · · , 6. And the domain Ω is partitioned into several regular disjoint elements for each mesh size h; see Figure 4.3, where h = 1/10, 1/20, 1/40 and h = 1/80. We choose #S(K) uniformly for all elements as in the second row of Table 4.1.
In Figure 4.4, we present the errors measured in both the DG norm and the L 2 norm. It is clear that the convergence rate in the DG norm is m − 1 for fixed m. The convergence rate in the L 2 norm is m + 1 for m ≥ 3, which converges quadratically when m = 2. Such convergence rate is consistent with the theoretical prediction in Theorem 3.1. Example 2. As we emphasize before, the local least-squares problem (2.4) is independent of the element geometry. In this example, we use the mesh generator in [17] to obtain a series of polygonal meshes from the Voronoi diagram of a given set of their reflections; see Figure 4.5. We also take (4.1) as the exact solution. The number #S(K) is chosen as in the third row of Table 4.1.
For each fixed m, we present the errors in both the DG norm and the L 2 norm against the number of elements; see Figure 4.6. It is clear that the numerical results still agree well with the theoretical results.
Example 3. In this test, we consider the biharmonic problem on Ω = (0, 1) 2 with the following boundary condition [18]: which is related to the bending of a simply supported plate. We take u(x, y) = sin(2πx) sin(2πy) as the exact solution. The mesh consists of a mixing of triangular and quadrilateral elements, which are generated by gmsh [19]; see Figure 4.7. The mesh size h varies equally from 1/10 to h = 1/80. We take #S(K) as in the last row of Table 4.1. The results in Figure 4.8 show the convergence rate for different m, which are also consistent with the theoretical results.
Example 4. In this test, we compare the C 0 IPG, IPDG and the proposed method by solving the biharmonic problem in the domain Ω = (0, 1) 2 with the same exact solution as for Example 1. The meshes used for this case are obtained by successively refining an initial mesh with h = 0.2. We measure the errors in both the broken H 2 norm and the L 2 norm. Figure 4.9 and Figure 4.10 show the performance of three methods by using spaces of polynomials of degree 2 and 3, respectively. We plot the errors in both norms against the number of degrees of freedom. It is clear that the proposed method behaves slightly better than other methods.        sin 2 (πx) sin 2 (πy) sin 2 (πz) and the function g D , g N and f are taken suitably. We build the reconstruction operator with different polynomial degrees m = 2, 3, 4. The number #S(K) is listed in Table 4.3.  The numerical results are presented in Table 4.4. The convergence rates also agree with the theoretical prediction.

Conclusion
We propose a new discontinuous Galerkin method to solve the biharmonic boundary value problem. A novelty of the method is a new discontinuous polynomial space that is reconstructed by solving local least-squares. The optimal error estimates in both the DG energy norm and the L 2 norm are proved, which are confirmed by a series of numerical examples with different complexity.