A Discontinuous Galerkin Method by Patch Reconstruction for Biharmonic Problem
Abstract.
We propose a new discontinuous Galerkin method based on the leastsquares patch reconstruction for the biharmonic problem. We prove the optimal error estimate of the proposed method. The twodimensional and threedimensional 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.
keyword: Leastsquares problem Reconstructed basis function discontinuous Galerkin method Biharmonic problem
MSC2010: 49N45; 65N21
1. Introduction
The biharmonic boundary value problem is a fourthorder 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 elements for this problem [18].
Recently, the discontinuous Galerkin (DG) method [10, 14, 16, 15, 8, 6] have been developed to solve the biharmonic problem. The DG method employs discontinuous basis functions that renders great flexibility in the mesh partition and also provides a suitable framework for 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 [11, 19]. As a compromise of the standard FEM and the DG method, Brenner et al [7, 3] developed the socalled interior penalty Galerkin method ( IPG). This method applies standard continuous 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 locates inside each element, which is a subspace 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 triangulation 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 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 space. Next, the interior penalty discontinuous Galerkin method with such reconstructed approximation space is proposed and analyzed in Section 3. In Section 4, we test the proposed method by several twodimensional and threedimensional biharmonic boundary value problems, which includes smooth solution as well as nonsmooth solution for various boundary conditions and different types of meshes.
2. Reconstruction operator
Let with be a bounded convex domain. Let be a collection of polygonal elements that triangulate . We denote all interior faces of as and the set of the boundary faces as , and then is then a collection of all dimensional faces of all elements in . Let and . We assume that satisfies the shaperegular conditions used in [1, 5], which reads: there exist

an integer number independent of ;

a real positive number independent of ;

a compatible subdecomposition into shaperegular triangles or quadrilaterals, or mix of both triangles and quadrilaterals,
such that

any polygon admits a decomposition formed by less than triangles;

any triangle is shaperegular in the sense that there exists satisfying , where is the radius of the largest ball inscribed in .
The above regularity assumptions lead to some useful consequences, which will be extensively used in the later analysis.

For any triangle , there exists that depends on and such that .

[Agmon inequality] There exists that depends on and , but independent of such that
(2.1) 
[Approximation property] There exists that depends on and , but independent of such that for any , there exists an approximating polynomial such that
(2.2) 
[Inverse inequality] For any , there exists a constant that depends only on and such that
(2.3)
Given the triangulation , we define the reconstruction operator as follows. First, for each element , we prescribe a point as the collocation point. In particular, we may specify as the barycenter of for all collocation points, although it could be more flexible. Second, we construct an element patch that consists of itself and some elements around . Denote by the number of elements inside . The element patch can be built in two ways. The first one is that we define the element patch in a recursive manner as in [13]. The other way is that we initialize with , and we add all Von Neumann neighbours of the elements in in a recursive manner until we collect enough elements in ; See Figure 2.1 for . Denote by the set of the collocation points belong to . It is clear that .
Let be the piecewise constant space associated with .
where is the polynomial space of degree not greater than .
For any function , we reconstruct a polynomial of degree on by solving the following local leastsquares:
(2.4) 
The uniqueness condition for the problem (2.4) relates to the location of the collocation points and . Following [13, 12], we make the following assumption:
Assumption 2.1.
For and ,
(2.5) 
The above Assumption guarantees the uniqueness of solution to problem (2.4) if is greater than . Hereafter, we assume that the assumption is always valid.
From (2.4) we define the global reconstruction operator for by restricting the polynomial on , .
We introduce some notations that are used in the traditional DG framework. First, we define the following broken Sobolev spaces:
where is the standard Sobolev space on element . The associated broken norm and seminorm are defined respectively by
where and are the standard Sobolev norm and seminorm on element , respectively.
For two neighbouring elements that share an face with , denote by and the outward unit normal vectors on , corresponding to and , respectively. For any function and that may be discontinuous across interelement boundaries, we define the average operator and the jump operator as
and
Here and . In the case , there exists an element such that , the definitions are modified as
and
We summarize the properties of the reconstruction operator proved in [13] and [12]. For any element , we define as
(2.6) 
By Assumption 2.5, we conclude that is finite. With some conditions on , we conclude that has a uniform upper bound, which is denoted by . We refer to [12, Theorem 2.2] for details. Given , we have the following interpolate estimate for the reconstruction operator.
Lemma 2.1.
Let and , there exists such that for ,
(2.7) 
and
(2.8) 
where .
Using the above lemma and Agmon inequality (2.1), we obtain the following interpolation estimate for the reconstruction operator in the energy norm.
Theorem 2.1.
Let . There exists a constant such that
(2.9) 
where .
3. Error Estimate for the Biharmonic Problem
Let us consider the biharmonic problem: find , such that
(3.1) 
where , denotes the unit outward normal vector to , , and 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 [2] for precise description of such results.
The interior penalty discontinuous Galerkin method [14, 15, 8] emplyes the following bilinear form : for any ,
The piecewise penalty parameters and are nonnegative and will be specified later on. The linear form is defined for all as
The discretized problem is to find such that
(3.2) 
The bilinear form is bounded and coercive if the penalty parameters are taken properly. We state this result in the following lemma and refer to [8] for the proof.
Lemma 3.1.
For any , let
(3.3) 
with positive constants and . For sufficiently big and , there exists constants and such that
(3.4)  
(3.5) 
We are ready to prove a priori error estimates for Problem (3.2).
Theorem 3.1.
Let be the solution of Problem (3.2), and assume that the solution to the Problem belongs to the broken Sobolev space . Furthermore, let the penalty functions and satisfy the condition (3.3) in Lemma 3.1. Then
(3.6) 
where and .
If the duality hypothesis holds [15], then
(3.7)  
(3.8) 
4. Numerical Results
In this section, we report some numerical examples to show the accuracy and performance of the proposed method. In all examples, we build element patches by the second method and use a sparse direct solver for the resulting sparse linear systems.
4.1. 2D Smooth Solutions
We firstly study the convergence rate for smooth solutions of twodimensional problems.
polynomial degree  2  3  4  5  6  
Example 1  9  15  22  29  38  
Example 2  9  16  23  32  45  
Example 3  9  20  28  38  49 
Example 1. Consider the biharmonic problem on the domain with Dirichlet boundary condition. The exact solution is taken as
(4.1) 
and , and the source term are chosen accordingly. The polynomial degree is taken by . And the domain is partitioned into several regular disjoint elements for each mesh size ; See Figure 4.1, where and . We select uniformly for all elements as in the second row of Table 4.1.
In Figure 4.2, we present the error measured in the DG norm and the norm. It is clear that the convergence rate in the DG norm is for fixed . The convergence rate in the L norm is for , which converges quadratically when . 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.3. We also take (4.1) as the reference solution. The number is chosen as in the third row of Table 4.1.
For each fixed , we present the error in the DG norm and the norm against the number of elements; See Figure 4.4. 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 with the following boundary condition [4]:
which is related to the bending of the simply supported plate. We take as the reference solution. The mesh consists of a mixing of triangular and quadrilateral elements, which are generated by gmsh[9]; See Figure 4.5. The mesh size varies equally from to . We take as in the last row of Table 4.1. The results in Figure 4.6 show the convergence rate different , which are also consistent with the theoretical results.
4.2. Lshaped domain with known exact solution
In this example, we study the performance of the proposed method with the problem with a corner singularity. Let be the Lshaped domain and we use triangular meshes; See Figure 4.7. Following [8], we let
in polar coordinate and impose Dirichlet boundary condition. At the corner the exact solution contains a singularity which indicates only belongs to for . The number is chosen as in the second row of Table 4.1.
In Table 4.2, we list the error measured in the DG norm and the norm against the mesh size for . Here we observe that the error in norm decreases at the rate while the error in DG norm decreases at the rate . It seems the convergence rate agree with that in [8].
Norm  Dofs  Dofs  Order  Dofs  Order  Dofs  Order  Dofs  Order  
Error  Error  Error  Error  Error  
2  3  4  1.17  4  1.23  4  1.20  5  1.19  
1  1  0.72  1  0.72  2  0.69  2  0.68  
3  4  4  1.47  4  1.35  5  1.19  5  1.21  
1  1  0.88  1  0.64  2  0.66  2  0.66  
4  3  4  1.76  4  1.52  5  1.28  5  1.22  
1  1  0.96  1  0.71  2  0.67  2  0.67 
4.3. 3D Smooth Solution
In this example, we solve a three dimensional biharmonic problem on a unit cube . The domain is partitioned into tetrahedral meshes with mesh size and by gmsh. The exact solution is chosen as and the function and are taken suitably. We build the reconstruction operator with different polynomial degrees . The number is listed in Table 4.3.
polynomial degree  2  3  4 
all  21  40  62 
The numerical results are presented in Table 4.4. The convergence rates are also agree with the theoretical prediction.
m  Norm  Dofs  Dofs  Order  Dofs  Order  Dofs  Order 
384  3072  24576  196608  
Error  Error  Error  Error  
2  2  2  2.36  3  2.10  3  2.05  
0  0  0.95  0  0.93  0  0.98  
3  2  3  3.15  4  3.91  5  3.97  
0  0  1.66  1  1.90  1  1.97  
4  2  4  4.98  5  5.19  7  5.09  
0  0  2.31  1  2.86  2  3.08 
5. Conclusion
We proposed a new discontinuous method based on combination of a leastsquares reconstruction procedure and the variational formulation of IPDG. The method is employed to solve the biharmonic boundary value problem. The optimal error estimate in the DG energy norm and the norm are proved. A series of numerical examples are reported to show the efficiency and wide application of the proposed method, which also confirm the theoretical results.
References
 P. F. Antonietti, L. Beirão da Veiga, and M. Verani. A mimetic discretization of elliptic obstacle problems. Math. Comp., 82(283):1379–1400, 2013.
 H. Blum and R. Rannacher. On the boundary value problem of the biharmonic operator on domains with angular corners. Math. Methods Appl. Sci., 2(4):556–581, 1980.
 S. C. Brenner. interior penalty methods. In Frontiers in Numerical Analysis—Durham 2010, volume 85 of Lect. Notes Comput. Sci. Eng., pages 79–147. Springer, Heidelberg, 2012.
 S. C. Brenner, P. Monk, and J. Sun. interior penalty Galerkin method for biharmonic eigenvalue problems. In Spectral and High Order Methods for Partial Differential Equations—ICOSAHOM 2014, volume 106 of Lect. Notes Comput. Sci. Eng., pages 3–15. Springer, Cham, 2015.
 F. Brezzi, K. Lipnikov, and V. Simoncini. A family of mimetic finite difference methods on polygonal and polyhedral meshes. Math. Models Methods Appl. Sci., 15(10):1533–1551, 2005.
 B. Cockburn, B. Dong, and J. Guzmán. A hybridizable and superconvergent discontinuous Galerkin method for biharmonic problems. J. Sci. Comput., 40(13):141–187, 2009.
 G. Engel, K. Garikipati, T. J. R. Hughes, M. G. Larson, L. Mazzei, and R. L. Taylor. Continuous/discontinuous finite element approximations of fourthorder elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity. Comput. Methods Appl. Mech. Engrg., 191(34):3669–3750, 2002.
 E. H. Georgoulis and P. Houston. Discontinuous Galerkin methods for the biharmonic problem. IMA J. Numer. Anal., 29(3):573–594, 2009.
 C. Geuzaine and J.F. Remacle. Gmsh: A 3D finite element mesh generator with builtin pre and postprocessing facilities. Internat. J. Numer. Methods Engrg., 79(11):1309–1331, 2009.
 P. Hansbo and M. G. Larson. A discontinuous Galerkin method for the plate equation. Calcolo, 39(1):41–59, 2002.
 T. J. R. Hughes, G. Engel, L. Mazzei, and M. G. Larson. A comparison of discontinuous and continuous Galerkin methods based on error estimates, conservation, robustness and efficiency. In Discontinuous Galerkin methods (Newport, RI, 1999), volume 11 of Lect. Notes Comput. Sci. Eng., pages 135–146. Springer, Berlin, 2000.
 R. Li, P. Ming, Z. Sun, and Z. Yang. A discontinuous finite element space by patch reconstruction. In Preparation, 2017.
 R. Li, P. Ming, and F. Tang. An efficient high order heterogeneous multiscale method for elliptic problems. Multiscale Model. Simul., 10(1):259–283, 2012.
 I. Mozolevski and E. Süli. A priori error analysis for the version of the discontinuous Galerkin finite element method for the biharmonic equation. Comput. Methods Appl. Math., 3(4):596–607, 2003.
 I. Mozolevski, E. Süli, and P. R. Bösing. version a priori error analysis of interior penalty discontinuous Galerkin finite element approximations to the biharmonic equation. J. Sci. Comput., 30(3):465–491, 2007.
 E. Süli and I. Mozolevski. version interior penalty DGFEMs for the biharmonic equation. Comput. Methods Appl. Mech. Engrg., 196(1316):1851–1863, 2007.
 C. Talischi, G. H. Paulino, A. Pereira, and I. F. M. Menezes. PolyMesher: a generalpurpose mesh generator for polygonal elements written in Matlab. Struct. Multidiscip. Optim., 45(3):309–328, 2012.
 O. C. Zienkiewicz, R. L. Taylor, and D. D. Fox. The Finite Element Method for Solid and Structural Mechanics. Elsevier/Butterworth Heinemann, Amsterdam, seventh edition, 2014.
 O. C. Zienkiewicz, R. L. Taylor, S. J. Sherwin, and J. Peiró. On discontinuous Galerkin methods. Internat. J. Numer. Methods Engrg., 58(8):1119–1148, 2003.