A $C^0$-Weak Galerkin Finite Element Method for the Biharmonic Equation

Abstract

A $C^0$-weak Galerkin (WG) method is introduced and analyzed in this article for solving the biharmonic equation in 2D and 3D. A discrete weak Laplacian is defined for $C^0$ functions, which is then used to design the weak Galerkin finite element scheme. This WG finite element formulation is symmetric, positive definite and parameter free. Optimal order error estimates are established for the weak Galerkin finite element solution in both a discrete $H^2$ norm and the standard $H^1$ and $L^2$ norms with appropriate regularity assumptions. Numerical results are presented to confirm the theory. As a technical tool, a refined Scott-Zhang interpolation operator is constructed to assist the corresponding error estimates. This refined interpolation preserves the volume mass of order $(k+1-d)$ and the surface mass of order $(k+2-d)$ for the $P_{k+2}$ finite element functions in $d$-dimensional space.

Appendix: A Mass-Preserving Scott-Zhang Operator

We will prove the existence of an interpolation $Q_0$ used in (16) and in the previous section, which is a special Scott-Zhang operator [26]. The new Scott-Zhang operator preserves the mass on each element and on each face, of four orders and three orders less, respectively, when interpolating $H^1(\Omega )$ functions to the finite element $V_h$ functions. We shall derive the optimal-order approximation properties for the interpolation in the section, which leads to a quasi-optimal convergence of the weak Galerkin finite element method (13).

The original Scott-Zhang operator maps $u\in H^1(\Omega )$ functions to $C^0$-Lagrange finite element functions, preserving the zero boundary condition if $u\in H^1(\Omega )$. It is an Lagrange interpolation. All the Lagrange nodes ([6]) on one element are classified into three types:

$$\begin{aligned}&\hbox {corner nodes} \, c_j: \quad \quad \, \, \hbox { 3 vertex nodes in 2D, or all edge nodes in 3D}, \\&\hbox {middle nodes} \ m_j: \quad \quad \hbox {all mid-edge nodes in 2D, or mid-triangle nodes in 3D}, \\&\hbox {internal nodes} \ i_j: \quad \quad \hbox {all internal nodes in the triangle/tetrahedra. } \end{aligned}$$

The three types of nodes are illustrated in Figs. 1 and 2. In simple words, $\{c_j\}$ are nodes shared by possibly more than two elements, $\{m_j\}$ are nodes shared by no more than two elements, and $\{i_j\}$ are nodes internal to one element.

Fig. 1

Averaging patches ($C_j$) in 2D (an edge or a triangle), for corner nodes ($c_j$), middle nodes ($m_j$) and internal nodes ($i_j$), in 2D

Fig. 2

Averaging patches $(C_j)$ in 3D (a triangle or a tetrahedron), for corner nodes $(c_j)$, middle nodes $(m_j)$ and internal nodes $(i_j)$, in 3D

A Lagrange nodal basis function $\phi _j$ is a $P_{k+2}$ polynomial which assumes value 1 at one node $c_j$, but vanishes at all other $\dim P_{k+2}^d-1$ nodes. For example, a $P_4^2$ nodal basis function on the reference triangle $\{ 0 \le x, y, 1-x-y \le 1 \}$, at node $(1/4,0)$, c.f. Fig. 3, is

$$\begin{aligned} \phi _2(x,y) = \frac{x (1-x-y)(3/4-x-y)(2/4-x-y)}{(1/4) (1-1/4-0)(3/4-1/4-0)(2/4-1/4-0)}. \end{aligned}$$

(46)

The restriction of a nodal basis $\phi _j$ on a lower dimensional simplex, a triangle or an edge or a vertex, is also a nodal basis function on that lower dimensional finite element. For example, this node basis function (46) is the restriction of the following 3D nodal basis function (at node $(1/4,0,0)$ on tetrahedron $\{ 0 \le x, y, z, 1-x-y-z \le 1 \}$) on the reference triangle,

$$\begin{aligned} \phi _j(x,y,z) = \frac{x (1-x-y-z)(3/4-x-y-z)(2/4-x-y-z)}{(1/4) (1-1/4-0-0)(3/4-1/4-0-0)(2/4-1/4-0-0)}. \end{aligned}$$

(47)

The restriction of 2D basis function $\phi _2$ in (46) in 1D is, c.f. Fig. 3,

$$\begin{aligned} \phi _{j'}(x) = \frac{x (1-x )(3/4-x )(2/4-x )}{(1/4) (1-1/4 )(3/4-1/4 )(2/4-1/4 )}. \end{aligned}$$

(48)

Fig. 3

A 2D nodal basis $\phi _2$ (46) and its restriction in 1D, $\phi _{j'}$ (48)

On each element $T$ (an edge, a triangle, or a tetrahedron), the $P_k$ Lagrange basis $\{\phi _j\}$ has a dual basis $\{\psi _j\in P_k^d \}$, satisfying

$$\begin{aligned} \int \limits _T \phi _j \psi _{j'} dx = \delta _{jj'}=\left\{ \begin{array}{ll} 1 &{} \hbox { if } j=j', \\ 0 &{} \hbox { if } j\ne j'. \end{array}\right. \end{aligned}$$

(49)

In other words, if writing $\{ \psi _j\}$ as linear combinations of Lagrange basis $\{\phi _j\}$, the coefficients are simply the inverse matrix of the mass matrix, the $L^2$-matrix of $\{\phi _j\}$. For example, the dual basis function $\psi _2$ for the nodal basis function $\phi _2$ in (46) (2D) is

$$\begin{aligned} \psi _2^{[2D]}(x,y)&= \frac{2835}{4} x - \frac{12285}{4}x^2 + \frac{8505}{2} x^3 + 8505 x^2 y +\frac{8505}{2} xy^2 \\&- 1890 x^4 - 5670 x^3 y - 5670 x^2 y^2 - 1890 x y^3. \end{aligned}$$

We can compute the dual of $\psi _{j'}$ in (48) in 1D to get

$$\begin{aligned} \psi _{j'}^{[1D]}(x)&= -\frac{485}{128} +\frac{2865}{32} x -\frac{21105}{64}x^2 +\frac{13615}{32}x^3 -\frac{11655}{64}x^4\\&= -\frac{485}{128}\phi _0+ \frac{64225}{16384}\phi _1 +\frac{345}{1024}\phi _2 -\frac{4255}{16384}\phi _3 -\frac{85}{128}\phi _4,\nonumber \end{aligned}$$

(50)

where $\phi _i$ is the nodal basis on $[0,1]$ at $x_i=i/4, i=0,1,2,3,4.$ The dual function $\psi _{j'}^{[1D]}(x)$ in (50) is plotted in Fig. 4. Similarly we can compute the dual basis function for (47) in 3D. We note that both Lagrange nodal basis and its dual basis are affine invariant. That is, the Lagrange basis on the reference triangle is also the Lagrange basis on a general triangle after an affine mapping. For simplicity, we use the same notations $\phi _j$ and $\psi _j^{[2D]}$ for the nodal basis and the dual basis functions on the reference triangle and on a general triangle.

Fig. 4

A 1D nodal basis $\phi _{j'}$ (48) and its dual basis in 1D, $\psi _{j'}^{[1D]}$ (50), $\int _0^1 \phi _{j'} \psi _{j'}^{[1D]}=1$

We now define the Scott-Zhang interpolation operator:

$$\begin{aligned} Q_0 : H^1 (\Omega ) \rightarrow V_h, \end{aligned}$$

where $V_h$ is the $C^0$-$P_{k+2}$ finite element space defined in (9). $Q_0v$ is defined by the nodal values at three types nodes.

• 1. For each corner node $c_j$ (shared by possibly more than two elements), we select one boundary $(d-1)$-dimensional simplex $C_j$ if $c_j$ is a boundary node, or any one $(d-1)$-dimensional face simplex $C_j$ on which the node is, as $c_j$’s averaging patch. C.f. Fig. 1, the boundary node $c_j$ has a boundary edge $C_j$, while the corner node $c_{j''}$ can choose any one of four edges passing it, as its averaging patch. In Fig. 2, a corner node $c_j$ has a triangle $C_j$ as its averaging patch. In both 2D and 3D, we use a same definition

  $$\begin{aligned} Q_0 v(c_j) = \int \limits _{C_j} \psi _j^{[(d-1)D]} v(x) d x. \end{aligned}$$

  (51)

• 2. For each middle node $m_{j'}$, the averaging patch $C_j$ is the unique $(d-1)$-dimensional simplex $C_j$ containing $m_j$, see Figs. 1 and 2. The interpolated nodal value is then determined by the unique solution of linear equations:

  $$\begin{aligned}&\int \limits _{C_j}\!\!\Bigg ( \sum _{m_{j'} \in C_j} Q_0 v(m_{j'})\phi _{j'}(x) \Bigg ) p_{i}(x) dx \\&\quad = \int \limits _{C_j}\!\!\Bigg ( v(x) - \sum _{c_{j } \in C_j} Q_0 v(c_{j })\phi _{j }(x) \Bigg ) p_{i}(x) dx \nonumber \end{aligned}$$

  (52)

  for all degree $(k+2-d)$ polynomials $p_{i}$ on ($d-1$)-simplex $C_j$, where $Q_0v(c_j)$ is defined in (51). In 2D, c.f. Fig. 1, after we determine the nodal values at the two end points ($c_j$), we determine the middle-edge points’ nodal value by (52).

• 3. After determine all nodal values on the surface of each element, we define the interpolation inside the element: The nodal values at internal nodes ($Q_0(i_{j''})$) are determined by the unique solution of the following linear equations

  $$\begin{aligned}&\int \limits _{C_j}\!\!\Bigg ( \sum _{i_{j''} \in C_j}Q_0 v(i_{j''})\phi _{j''}\Bigg ) p_{i} dx\\&\quad = \int \limits _{C_j} \!\!\Bigg (v - \sum _{m_{j'} \in C_j} Q_0 v(m_{j'})\phi _{j'} - \sum _{c_{j } \in C_j} Q_0 v(c_{j })\phi _{j } \Bigg ) p_{i} dx \nonumber \end{aligned}$$

  (53)

  for all degree $(k+1-d)$ polynomials $p_{i}$ on $d$-simplex $C_j$.

By (51)–(53), the (refined) Scott-Zhang interpolation is

$$\begin{aligned} Q_0 v = \sum _{x_j\in \mathcal{N}_h} Q_0v(x_j) \phi _j(x), \end{aligned}$$

(54)

where $\mathcal{N}_h$ is the set of all $C^0$-$P_{k+2}$ Lagrange nodes of triangulation $\mathcal{T }_h$.

Remark 8.1

If all corner nodes on a $(d-1)$-simplex $C_j$ have selected $C_j$ as their averaging patch in (51), then the solution of (52) is the $L^2$-projection of $v$ on $C_j$, i.e.,

$$\begin{aligned} Q_0 v(m_{j'}) = \int \limits _{C_j} \psi _{j'}^{[(d-1)D]} v(x) d x. \end{aligned}$$

(55)

Lemma 8.1

The Scott-Zhang interpolation operator (54) is well-defined, i.e., the linear systems of Eqs. (52) and (53) both have unique solutions.

Proof

For the linear system of Eq. (52), we change the $P_{k+2-d}^{d-1}$ basis functions ($\{p_i=1,x,...,x^k\}$ when $d=2$, or $\{p_i=1,x,y,x^2,xy,...,y^{k-1}\}$ when $d=3$) uniquely as linear combinations of the Lagrange basis functions on the subinterval $C_j^s$ ($d=2$) or the subtriangle $C_j^s$ ($d=3$), c.f. Fig. 5.

$$\begin{aligned} p_i = \sum _j c_{i,j} \phi _j^s, \quad i=1,2,...,\dim \left( P_{k+2-d}^{d-1}\right) . \end{aligned}$$

(56)

The nodal basis functions on a simplex $C_j$ and its subsimplex $C_j^s$ differ by a bubble functions:

$$\begin{aligned} \phi _j = \phi _{j'}^s \frac{ b(\mathbf x ) }{ b (\mathbf x _j)} , \end{aligned}$$

(57)

where $b(\mathbf x )$ is the bubble functions assuming 0 on the boundary of $C_j$. For example, c.f. Fig. 5, when $d=2$ and $C_j=[0,1]$,

$$\begin{aligned} b( x)&= x(1-x),\\ \phi _1^s(x)&= \frac{(x-2/4)(x-3/4)}{(1/4-2/4)(1/4-3/4)}, \\ \phi _2 (x)&= \phi _1^s(x) \frac{b(x)}{b(1/4)}. \end{aligned}$$

Fig. 5

Lagrange nodal basis $\phi _j^s$ on subinterval $(d=2)$ or subtriangle ($d=3$), c.f., (57)

By the change of basis, (56) and (57), the linear system (52) is equivalent to the following weighted-mass linear systems:

$$\begin{aligned} \sum _{m_{j'}\in C_j} Q_0 v(m_{j'})\int \limits _{C_j}\phi _{j'} \phi _{i}^s dx = \int \limits _{C_j} v \phi _{i}^s dx&- \sum _{j\in C_j} Q_0 v(c_j) \int \limits _{C_j} \phi _{j} \phi _{i}^s dx,\\&i=1,2,...,\dim \left( P_{k+2-d}^{d-1}\right) .\nonumber \end{aligned}$$

(58)

The coefficient matrix in (58) is the mass matrix on simplex $C_j$ (c.f. Fig. 5) with a positive weight:

$$\begin{aligned} a_{i,j'} = \int \limits _{C_j} \phi _{i}^s \phi _{j'} d\mathbf x = \int \limits _{C_j} \phi _{i}^s \phi _{j}^s w(\mathbf x ) d\mathbf x , \end{aligned}$$

where

$$\begin{aligned} w(\mathbf x ) = \frac{b(\mathbf x )}{b(\mathbf x _{j'})} >0 \quad \text {in interior}(C_j). \end{aligned}$$

As the Lagrange basis $\{ \phi _i^s \}$ (on the subsimplex) are linearly independent, the mass (with weight) matrix in (58) is invertible, and the equivalent linear system (52) has a unique solution too. For example, when $d=2$ and $k=2$, the coefficient matrix in (58) and its inverse are

$$\begin{aligned} \left( \begin{array}{l@{\quad }l@{\quad }l} 152/315 &{} -16/63 &{}8/63 \\ -4/21 &{} 18/35 &{} -4/21 \\ 8/63 &{} -16/63 &{} 152/315\\ \end{array}\right) ^{-1} =\left( \begin{array}{l@{\quad }l@{\quad }l} 2655/1024 &{} 75/64 &{} -225/1024 \\ 225/256 &{} 45/16 &{} 225/256 \\ -225/1024 &{} 75/64 &{} 2655/1024\\ \end{array}\right) . \end{aligned}$$

By the same argument, lifting the space dimension by 1, we can show that (53) has a unique solution too. In fact, the system (53) when $d=2$ is the same system (52) with $d=3$ there. $\square$

Lemma 8.2

If $d=2$, the Scott-Zhang interpolation operator (54) preserves the volume mass of order $k-1$ and the face mass of order $k$, i.e.,

$$\begin{aligned} \int \limits _{T} (v - Q_0 v)p_i dx&= 0 \quad \forall T\in \mathcal{T }_h, \ p_i \in P_{k-1} (T), \end{aligned}$$

(59)

$$\begin{aligned} \int \limits _{E} (v - Q_0 v)p_i dx&= 0 \quad \forall E\in \mathcal{E }_h, \ p_i \in P_{k}(E), \end{aligned}$$

(60)

where $\mathcal{E }_h$ is the set of edges in triangulation $\mathcal{T }_h$.

Proof

By the construction (52), we have

$$\begin{aligned} \int \limits _{E} (v - Q_0 v)p_i dx = \int \limits _{E} \Bigg ( v - \sum _{m_{j'}\in E} Q_0v(m_{j'}) \phi _{j'}- \sum _{c_{j }\in E} Q_0v(c_{j }) \phi _{j }\Bigg )p_i dx=0. \end{aligned}$$

That is (60). Here, because we have two missing dof (degrees of freedom) at the two end points of each edge, the polynomial degree in mass preservation is reduced by two, from $(k+2)$ to $k$. Similarly, (59) follows (53). Here, the polynomial degree reduction is three as each triangle has three edges where the interpolation values are not free (not determined by (53)). $\square$

Lemma 8.3

If $d=3$, the Scott-Zhang interpolation operator (54) preserves the volume mass of order $k-2$ and the face mass of order $k-1$, i.e.,

$$\begin{aligned} \int \limits _{T} (v - Q_0 v)p_i dx&= 0\quad \forall T\in \mathcal{T }_h, \ p_i \in P_{k-2} (T), \end{aligned}$$

(61)

$$\begin{aligned} \int \limits _{E} (v - Q_0 v)p_i dx&= 0 \quad \forall E\in \mathcal{E }_h, \ p_i \in P_{k-1} (E), \end{aligned}$$

(62)

where $\mathcal{E }_h$ is the set of face triangles in the tetrahedral grid $\mathcal{T }_h$.

Proof

As each triangle $E$ has three edges, where the interpolation is not determined possibly by function value on neighboring triangles, we lose dof’s on the three edges in the interpolation. That is, we lose three orders in face-mass conservation in (52). (62) is simply another expression of (52), as in the proof of Lemma 8.2. By (53), (61) follows. Here the polynomial-degree deduction in mass conservation is 4, due to 4 face-triangles each tetrahedron. $\square$

Remark 8.2

By Lemmas 8.2 and 8.3, the mass preservation on element and on faces is one order higher than the requirement (16), when $d=2$. When $d=3$, (61) and (62) imply (16).

Theorem 8.1

The Scott-Zhang interpolation operator (54) is of optimal order in approximation, i.e., when $k\ge -1$,

$$\begin{aligned} \Vert v-Q_0v\Vert +h\Vert v-Q_0v\Vert _1 \le C h^{k+3}\Vert v\Vert _{k+3} \quad \forall v \in H^{k+3}(\Omega ). \end{aligned}$$

(63)

Further, when $k\ge 0$,

$$\begin{aligned} \left( \sum _{T\in \mathcal{T }_h} h^2 |v-Q_0v |_{H^2(T)}^2 \right) ^{1/2} \le C h^{k+3} \Vert v\Vert _{k+3} \quad \forall v \in \cap H^{k+3}(\Omega ). \end{aligned}$$

(64)

Proof

The Scott-Zhang operator preserves a degree $(k+2)$ polynomial on the star union of an element $T$,

$$\begin{aligned} S_T = \cup _{\overline{T'}\cap T \ne \emptyset } \overline{T'},\quad T, T' \in \mathcal{T }_h. \end{aligned}$$

That is,

$$\begin{aligned} Q_0 v = v \quad \forall v \in P_{k+2}^d(S_T), \end{aligned}$$

(65)

when $k\ge -1$. (65) is shown in three steps. First, by (51), when $v\in P_{k+2}$, the dual basis defines

$$\begin{aligned} Q_0 v(c_j) = v(c_j), \end{aligned}$$

(66)

at all corner nodes. In the second step, by (66) and (58), (52) holds for $p_i=\psi _{i}^s$ too where $m_{i'}$ are all middle nodes on $C_j$. That is,

$$\begin{aligned} \int \limits _{C_j} \Bigg (v-\sum _{j\in C_j} Q_0(c_j)\phi _j -\sum _{m_{j'}} Q_0 v(m_{j'}) \phi _{j'}\Bigg ) \phi _i^s dx = 0. \end{aligned}$$

(67)

By (66), as $v\in P_{k+2}$,

$$\begin{aligned} v-\sum _{j\in C_j} Q_0(c_j)\phi _j = v_c b(\mathbf x ) \quad \hbox { for some } v_c \in P_{k+2-d}^{d-1}, \end{aligned}$$

where $b(\mathbf x )$ is a bubble function, cf. (57). For the middle node basis functions, we have also, c.f. (57),

$$\begin{aligned} \phi _{j'}= \phi _{j''}^s b(\mathbf x )/b(\mathbf x _{j'}) \quad \forall m_{j'}\in C_j. \end{aligned}$$

Thus (67) implies, where $\sum _{m_{j'}} Q_0 v(m_{j'}) \phi _{j'}=v_m b(\mathbf x )$ for some $v_m\in P_{k+2-d}^{d-1}$,

$$\begin{aligned}&\int \limits _{C_j} (v_b-v_m) \phi _i^s b(\mathbf x ) d\mathbf x = 0, \\&\Rightarrow v_b-v_m\equiv 0 \quad \hbox {and } \ v=Q_0 v \quad \hbox { on } \ C_j. \end{aligned}$$

Therefore, at all middle nodes,

$$\begin{aligned} Q_0 v(m_{j'}) = v(m_{j'}). \end{aligned}$$

(68)

In the third step, by (53) and the same argument in the second step,

$$\begin{aligned} Q_0 (i_{j''}) = v(i_{j''}), \end{aligned}$$

at all internal nodes. Thus $Q_0 v=v\in P_{k+2}$.

We then use the standard scaling argument (on the dual basis functions) and the Sobolev inequality, as in Theorem 3.1 of [26], it follows that

$$\begin{aligned} |Q_0 v |_{H^1(T) } \le C \Vert v\Vert _{H^1(S_T) } \quad \forall v\in H^1 (\Omega ). \end{aligned}$$

The above stability result leads directly to the optimal-order approximation (63), following the standard argument (i.e., by (65) and the existence of local Taylor polynomials, c.f. for example, [6]), as shown in Theorem 4.1 of [26]. We note again that the Scott-Zhang operator here is a refined version of the Scott-Zhang operator in [26]. After showing the local preservation of $P_{k+2}$ polynomials above, the proof of the theorem is the same as that in [26]. (64) is (4.4) in [26], with $p=q=2, m=2$, and $l=k+3$ there. $\square$