- Email: [email protected]
Weak Galerkin based a posteriori error estimates for second order elliptic interface problems on polygonal meshes
Journal of Computational and Applied Mathematics 361 (2019) 413–425
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Weak Galerkin based a posteriori error estimates for second order elliptic interface problems on polygonal meshes Lin Mu 1 Computer Science and Mathematics Division Oak Ridge National Laboratory, Oak Ridge, TN, 37831, USA
article
info
Article history: Received 17 August 2017 Received in revised form 20 March 2019 MSC: primary 65N15 65N30 secondary 35J50
a b s t r a c t In this paper, we present a posteriori error estimate of weak Galerkin (WG) finite element methods based on the second order elliptic interface problems. This estimate can be applied to polygonal meshes or meshes with hanging nodes. The reliability and efficiency of the designed error estimator has been proved in this work. Extensive numerical tests are performed to validate our algorithm. These results demonstrate the effectiveness of the adaptive mesh refinement guided by the proposed error estimator. © 2019 Elsevier B.V. All rights reserved.
Keywords: Weak Galerkin Finite element methods Second-order elliptic interface problems A posterior error estimate Polygonal meshes
1. Introduction The weak Galerkin (WG) method was first proposed by Wang and Ye in [1] as a general framework for solving partial differential equations (PDE). This method is a natural extension of finite element methods for the discontinuous functions, with finite element formulation as simple as the weak form of the corresponding PDE. The key in WG methods is to replace the classical derivatives by the weakly defined discrete derivatives. Then a parameter-free stabilizer is added to produce a stable, symmetric, and positive definite formulation. Advantages of WG methods include the flexibility of employing polygonal meshes and naturally way to design high-order numerical schemes by using polynomials with higher degree. There are two different error analysis approaches: one is a priori error estimate and the other one is a posteriori error estimate. A priori error analysis offers theoretical conclusions about the convergence order of the errors measured in certain norms with respect to mesh size h. Due to the unknown exact solution is included in the measurement, the a priori error is not actually computable. The a posteriori error estimate extracts a computable quantity as estimator which is equivalent to the interested quantity measured in the chosen norm. If the solutions contain singularity, a posteriori error estimator can be used to drive mesh refinement and thus generate the adaptive finite element (AFE) methods. Recently, AFE methods are widely used in computational science and engineering to obtain better accuracy with minimum effort when singularities are included in the problems. The savings are achieved through only refine the mesh at the needed places containing bigger errors. In the adaptive mesh methods, several refinement rules have been developed, including red and green refinement [2], longest edge bisection [3,4] and newest vertex bisection [5]. Since the refinement may produce hanging nodes in the mesh, one needs further refine the adjoint neighbors for keeping the E-mail address: [email protected]. 1 This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research. This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. https://doi.org/10.1016/j.cam.2019.04.026 0377-0427/© 2019 Elsevier B.V. All rights reserved.
414
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
Fig. 2.1. Depiction of a shape-regular polygonal element ABCDEFA.
shape regular requirement. The limitations of quad/triangular mesh in two dimensional space and unnecessary neighbor element refinement in AFE motivate our project of investigate the a posteriori error estimate on the polygonal meshes. In recent research, many numerical schemes have been developed and analyzed on general polytopal meshes [6–18]. However it is still hard to find a posteriori estimators work for general polygonal/polyhedral meshes [19–21]. Even for the polygonal and polyhedral finite element methods, due to the technique difficulties in analysis, most a posteriori error analysis is still based on simplicial meshes and strong assumptions [22–29]. In reference [19], the authors designed a polygonal error estimator for the virtual element method. Reference [30] developed a simple posteriori error estimator for weak Galerkin finite element approximation for Poisson equations. In this paper, we shall extend the a posteriori error estimate to study the second order elliptic interface problems. In this work, our diffusion coefficient is assumed to be piecewise constant, which may be different values at different sub-regions. For simplicity, we consider a simple model problem that seeks an unknown function u satisfying
−∇ · κ (x)∇ u = f , in Ω , u = 0, on ∂ Ω ,
(1.1) (1.2)
where Ω is a polytopal domain in Rd (polygonal or polyhedral domain for d = 2, 3). The diffusion coefficient κ (x) is positive and piecewise constant on the subdomains of Ω with possible large jumps across interface. Denote κmin and κmax as the minimum and maximum values of κ . The rest of this paper is organized as follows: Weak Galerkin finite element methods are introduced in Section 2. Section 3 introduces a proposed a posteriori error estimator and is contributed to analyze reliability and efficiency of the proposed estimator. Three numerical examples are present in Section 4. The paper ends with conclusions and future work plan as Section 5. 2. Weak Galerkin finite element schemes 2.1. Finite element spaces We will adopt the standard notation and definitions for Sobolev spaces [31]. In particular, given an open region B ⊂ Rd , denote the Sobolev space by W s,r (B) on B, where s is a real number. When r = 2, W s,2 (B) = H s (B) denotes a Hilbert space equipped with the norm ∥ · ∥s,B . When s = 0, H 0 (B) coincides with the standard Lebesgue L2 (B) space equipped with the inner product (·, ·)T and the L2 -norm ∥ · ∥B . If B = Ω , the subscript B will be omitted. We shall use similar notations for ∂ B. Besides, ⟨·, ·⟩∂ B denotes the inner product on the edge ∂ B. The space H(div, B) is defined for vector-valued functions, which have square integrable divergence on B, and the norm is defined as follows:
∥v∥H(div,B) := (∥v∥2B + ∥∇ · v∥2B )1/2 . Let Th be a partition of the domain Ω consisting of elements which are closed and simply connected polygons in two dimension or polyhedra in three dimension; see Fig. 2.1. Let Eh be the set of all edges or flat faces in Th , and Eh0 = Eh \∂ Ω be the set of all interior edges or flat faces. Denote by hT the diameter for every element T ∈ Th and h = maxT ∈Th hT the mesh size for Th . We need some shape regularity assumptions for the partition Th described as below (cf. [32,33]).
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
415
A1: Assume that there exist two positive constants ϱv and ϱe such that for every element T ∈ Th we have
ϱe hde −1 ≤ |e|
ϱv hdT ≤ |T |,
(2.1)
for all edges or flat faces of T . A2: Assume that there exists a positive constant κ such that for every element T ∈ Th we have
κ hT ≤ he
(2.2)
for all edges or flat faces e of T . A3: Assume that the mesh edges or faces are flat. We further assume that for every T ∈ Th , and for every edge/face e ∈ ∂ T , there exists a pyramid P(e, T , Ae ) contained in T such that its base is identical with e, its apex is Ae ∈ T , and its height is proportional to hT with a proportionality constant σe bounded away from a fixed positive number σ ∗ from below. In other words, the height of the pyramid is given by σe hT such that σe ≥ σ ∗ > 0. The pyramid is also assumed to stand up above the base e in the sense that the angle between the vector xe − Ae , for any xe ∈ e, and the outward normal direction of e is strictly acute by falling into an interval [0, θ0 ] with θ0 < π2 . A4: Assume that each T ∈ Th has a circumscribed simplex S(T ) that is shape regular and has a diameter hS(T ) proportional to the diameter of T ; i.e., hS(T ) ≤ γ∗ hT with a constant γ∗ independent of T . Furthermore, assume that each circumscribed simplex S(T ) interests with only a fixed and small number of such simplices for all other elements T ∈ Th . Assume the interface for coefficients κ (x) does not cut element T ∈ Th . For simplicity of notation, in the following we shall assume d = 2. The results can be extended to three dimensional space. For a scalar-valued function v , we define the curl operator by
∇ × v = (−
∂v ∂v , ). ∂ x2 ∂ x1
For a given integer k ≥ 1, Pk (T ) denotes the space of polynomials of degree k on element T and Pk (e) denotes the space of polynomials of degree k on the edge e. Let Vh be the weak Galerkin finite element space associated with Th defined as follows Vh = {v = {v0 , vb } : v0 |T ∈ Pk (T ), vb |e ∈ Pk (e), e ∈ Eh , T ∈ Th } and Vh0 = {v : v ∈ Vh , vb = 0 on ∂ Ω }. We would like to emphasize that v0 is piecewise polynomial on each element T and vb is piecewise polynomial on each edge e. 2.2. WG finite element approximation For any weak function v = {v0 , vb }, a weak gradient ∇w v ∈ [Pk−1 (T )]2 is defined on T as the unique solution satisfying (∇w v, τ )T = −(v0 , ∇ · τ )T + ⟨vb , τ · n⟩∂ T ,
∀τ ∈ [Pk−1 (T )]2 .
(2.3)
Now the following introduces two bilinear forms on Vh : ∀v, w ∈ Vh s(v, w ) =
∑
sT (v, w ),
T ∈Th
(κ∇w v, ∇w w ) =
∑
(κT ∇w v, ∇w w )T ,
T ∈Th 1 where sT (v, w ) = h− T ⟨κT (v0 − vb ), w0 − wb ⟩∂ T and κT = κ|T .
Weak Galerkin Algorithm 1. Find uh ∈ Vh0 satisfying the following equations: (κ∇w uh , ∇w v ) + s(uh , v ) = (f , v0 ), where (f , v0 ) :=
∑
∫ T ∈Th
T
f v0 dT .
∀ v ∈ Vh0 ,
(2.4)
416
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
The energy norm can be defined as,
|||v|||2 = (κ∇w v, ∇w v ) + s(v, v ).
(2.5)
The existence and uniqueness of numerical solution can be found in [32]. Furthermore, we cite the following theorem regarding a priori error estimate [32]. Theorem 2.1. Let uh ∈ Vh and u ∈ H k+1 (Ω ) be the weak Galerkin finite element solution and exact solution to problems (2.4) and (1.1)–(1.2). Then, there exists a constant C independent of mesh size h such that
|||u − uh ||| ≤ Chk ∥u∥k+1 , ∥u − uh ∥ ≤ Chk+1 ∥u∥k+1 .
(2.6) (2.7)
3. A posteriori error estimate Let fh be the L2 projection of f to Vh . Then a local estimator is defined as follows
ηT2 = sT (uh , uh ) + osc2 (f , T ),
(3.1)
where osc(f , T ) is a high order local data oscillation defined by osc2 (f , T ) = h2T ∥f − fh ∥2T . Define a global error estimator and data oscillation as
η2 =
∑
ηT2 ,
T ∈Th
osc(f , Th ) = 2
∑
osc(f , T )2 .
T ∈Th
Let T be an element with e as an edge. The following trace inequality will be also used in the following section. It is well known that there exists a constant C such that for any function g ∈ H 1 (T )
( 1 ) 2 2 ∥g ∥2e ≤ C h− T ∥g ∥T + hT ∥∇ g ∥T .
(3.2)
For each element T ∈ Th , denote by Q0 the L2 projection from L2 (T ) to Pk (T ), and by Qb the L2 projection from L2 (e) to Pk−1 (e). By an abuse of notation, we use here the same symbols for Q0 and Qb as global operators. Denote by Qh the L2 projection from [L2 (T )]2 to a local weak gradient space [Pk−1 (T )]2 . Define Qh u = {Q0 u, Qb u} ∈ Vh . The following lemma can be found in [32] Lemma 3.1. For any element T ∈ Th , we have the following commutative property for φ ∈ H 1 (T ),
∇w (Qh φ ) = Qh (∇φ ), ∇w φ = Qh (∇φ ).
(3.3) (3.4)
Theorem 3.2 (Reliability). Let uh ∈ Vh0 and u ∈ H01 (Ω ) be the solutions of (2.4) and (1.1)–(1.2) respectively. Then there exists a positive constant C independent of mesh size h such that:
|||u − uh ||| ≤ C η.
(3.5)
Proof. We shall apply Helmholtz decomposition first. It is well known [34] that for κ (∇w u − ∇w uh ) ∈ [L2 (Ω )]2 , there exist ψ ∈ H01 (Ω ) and φ ∈ H 1 (Ω ) such that
κ (∇w u − ∇w uh ) = κ∇ψ + ∇ × φ
(3.6)
and that
∥κ 1/2 (∇w u − ∇w uh )∥2 = ∥κ 1/2 ∇ψ∥2 + ∥κ −1/2 ∇ × φ∥2 . (3.4) and the definition of Qh implies: ∀q ∈ [Pk−1 (T )]2 , (κT ∇w u, q)T = (κT Qh (∇ u), q)T = (κT ∇ u, q)T .
(3.7)
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
417
From the above equality, and (3.6), one obtains,
∥κ 1/2 (∇w u − ∇w uh )∥2 =
∑
( κ T ( ∇ u − ∇ w uh ) , ∇ w u − ∇ w uh ) T
T ∈Th
=
∑
(∇ u − ∇w uh , κT (∇w u − ∇w uh ))T
(3.8)
T ∈Th
= (κ∇ψ, (∇ u − ∇w uh )) + (∇ × φ, ∇ u − ∇w uh ). By (3.3), the definition of Qh , and (f , Q0 ψ ) = (κ∇w uh , ∇w Qh ψ ) + s(uh , Qh ψ )
= (κ∇w uh , Qh (∇ψ )) + s(uh , Qh ψ ) = (κ∇w uh , ∇ψ ) + s(uh , Qh ψ )
(3.9)
it is derived that (κ∇w uh , ∇ψ ) = (f , Q0 ψ ) − s(uh , Qh ψ ).
(3.10)
Also, by the definition of Qb , we have 1/2
1 −1 2 h− T ∥κT (Q0 ψ − Qb ψ )∥∂ T = hT ⟨κT (Q0 ψ − Qb ψ ), Q0 ψ − Qb ψ⟩∂ T 1 = h− T ⟨κT (Q0 ψ − ψ ), Q0 ψ − Qb ψ⟩∂ T ( )1/2 ( )1/2 1/2 1/2 −1 −1 2 2 ≤ hT ∥κT (Q0 ψ − ψ )∥∂ T hT ∥κT (Q0 ψ − Qb )ψ∥∂ T
and thus by trace inequality and property of projection operator [33], 1/2
1/2
h−1 ∥κT (Q0 ψ − Qb ψ )∥2 ∂ T ≤ h−1 ∥κT (Q0 ψ − ψ )∥2∂ T 1/2
≤ C ∥κT ∇ψ∥2T .
(3.11)
Using (3.9), Cauchy–Schwarz inequality, inequality (3.11), and (3.7), we have (κ∇ψ, (∇ u − ∇w uh )) = (κ (∇ u − ∇w uh ), ∇ψ )
= (κ∇ u, ∇ψ ) − (κ∇w uh , ∇ψ ) = (f , ψ ) − (f , Q0 ψ ) + s(uh , Qh ψ ) ∑ = (f − fh , ψ − Q0 ψ ) + hT −1 ⟨κT (u0 − ub ), Q0 ψ − Qb ψ⟩∂ T T ∈Th
≤ C (osc(f , Th ) + s1/2 (uh , uh ))∥κ 1/2 ∇ψ∥ ≤ C (osc(f , Th ) + s1/2 (uh , uh ))∥κ 1/2 (∇w u − ∇w uh )∥. Since ∇ × φ ∈ H(div, Ω ) and u ∈
H01 (Ω ),
(3.12)
we have
( ∇ u − ∇ w uh , ∇ × φ ) = ( ∇ u, ∇ × φ ) − ( ∇ w uh , ∇ × φ )
= − ( ∇ w uh , ∇ × φ ) . It follows from (2.3), uh ∈ Vh0 and the integration by parts, (∇w uh , ∇ × φ )T = (∇w uh , Qh ∇ × φ )T
= = = =
⟨ub , Qh (∇ × φ ) · n⟩∂ T − (u0 , ∇ · (Qh ∇ × φ ))T ⟨ub , Qh (∇ × φ ) · n⟩∂ T − ⟨u0 , (Qh ∇ × φ ) · n⟩∂ T + (∇ u0 , Qh ∇ × φ )T ⟨ub − u0 , Qh (∇ × φ ) · n⟩∂ T + (∇ u0 , ∇ × φ )T ⟨ub − u0 , Qh (∇ × φ ) · n⟩∂ T + ⟨u0 , (∇ × φ ) · n⟩T
and ub is single valued on each edge e ∈ Eh , one have
∑ T ∈Th
( ∇ w uh , ∇ × φ ) T =
∑
(⟨ub − u0 , Qh (∇ × φ ) · n⟩∂ T + ⟨u0 − ub , ∇ × φ·n⟩∂ T ).
T ∈Th
The Cauchy–Schwarz inequality, the trace inequality (3.2) and the inverse inequality imply
⟨u0 − ub , Qh ∇ × φ · n⟩∂ T ≤ C ∥u0 − ub ∥∂ T ∥Qh ∇ × φ · n∥∂ T 1/2 −1/2 ≤ C sT (uh , uh )∥κT ∇ × φ∥T .
(3.13)
418
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
Using the inverse inequality and the fact ∇ · ∇ × φ = 0, we arrive at
∑
⟨ u0 − ub , ∇ × φ · n ⟩ ∂ T =
⟨u0 − ub , ∇ × φ · n⟩e
e⊂∂ T
∑
≤
1/2
−1/2
∥κT (u0 − ub )∥H 1/2 (e) ∥κT
∇ × φ · n∥H −1/2 (e)
e⊂∂ T
≤ C(
∑
−1/2
hT
1/2
−1/2
∥κT (u0 − ub )∥2e )1/2 ∥κT
∇ × φ∥H(div,T )
e⊂∂ T 1/2
−1/2
≤ C (sT (uh , uh ))∥κT
∇ × φ∥T .
Combining the above three estimates and taking a summation over T give (∇w uh , ∇ × φ ) ≤ C (s1/2 (uh , uh ))∥κ −1/2 ∇ × φ∥.
(3.14)
Using (3.8), (3.12), (3.13) (3.14) and (3.7), we have
∥k1/2 (∇w u − ∇w uh )∥ ≤ C η. Next we easily have s(u − uh , u − uh ) = s(uh , uh ) ≤ η2 . Thus we complete the proof.
□
Define the localized error on each element T as follows,
|||v|||2T = (∇w v, ∇w v )T + sT (v, v ). Since sT (v, v ) is only part of the local error |||v|||2T , then the following local lower bound is automatically obtained. Theorem 3.3 (Efficiency). The local estimator ηT as defined in (3.1) is bounded as the following inequality:
ηT2 ≤ |||u − uh |||2T + osc2 (f , T ).
(3.15)
4. Numerical example In this section, we shall validate the proposed a posteriori error estimate for several tests. First, we shall explore the convergence properties of errors measured in ||| · |||-norm (denoted by H 1 -error in tables and H 1 err in figures) and the estimator η. The effectivity of the estimator is defined as follows, Eff-index =
η . |||Qh u − uh |||
(4.1)
The polygon refinement will be generated by two different strategies:
• Mesh Type 1: connect the barycenter of each element to the mid-point of each edge corresponding to this element. This approach was also used in [19]. In this way, all the polygons will be refined to several quadrilaterals. An example can be found as Fig. 4.1(a). • Mesh Type 2: connect the barycenter of each element to the mid-point of each edge corresponding to this element, and mark the mid-points of the connecting edges. By connecting these marked mid-points, marked mid-points and the middle point of each edge, one can recursively refine the mesh. This approach was developed in [35]. An example can be found as Fig. 4.1(b). 4.1. Example 1 Let Ω = (0, 1)2 , assume the interface at x = 1/2, and the diffusion coefficient
κ (x) =
{
1, if x < 1/2, 2, if x ≥ 1/2,
and u = sin(2π x)2 cos(π y).
In this test, we shall check the effectivity index as defined in (4.1) for comparing our posteriori error estimator η and error |||Qh u − uh |||. For simplicity, we only consider linear WG element on the polygonal mesh. We start the WG simulation on the initial mesh as shown in Fig. 4.2(a), which is denoted as Mesh Level 1. Two different uniformly refinement techniques (Mesh Type 1 and Mesh Type 2) will be applied to generate a sequence of meshes. Two examples of Mesh Level 2 of two different refinement strategies can be found as Fig. 4.2(b)–(c). The error profiles, convergence results, and Eff-index are reported in Table 4.1. Because of the setting for exact solution in this problem, one can observe the errors, measured in ||| · |||-norm and estimator η, converge at the optimal rate in convergence at order O(h) with respect to mesh size h. Also, the last column presents the Eff-Index, which is shown as a constant, and thus confirms the equivalency of η and |||Qh u − uh |||.
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
419
Fig. 4.1. Two different refinements (a) Mesh Type 1; (b) Mesh Type 2.
Fig. 4.2. Example 1: (a) Mesh Level 1 for initial Mesh with #ELEM = 58, DOF = 486; (b) Mesh Level 2 of Mesh Type 1 with #ELEM = 288, DOF = 3528; (c) Mesh Level 2 of Mesh Type 2 with #ELEM = 346, DOF = 2814.
420
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425 Table 4.1 Example 1. Error profiles and effectivity index. Mesh Type 1 Mesh Mesh Mesh Mesh Mesh Mesh
Level Level Level Level Level
1 2 3 4 5
H 1 -error
Order
η
Order
Eff-Index
3.0999E+00 1.4237E+00 7.3099E−01 3.6861E−01 1.8475E−01
– 1.12 0.96 0.99 1.00
3.6210E+00 1.6143E+00 8.3385E−01 4.2212E−01 2.1191E−01
– 1.17 0.95 0.98 0.99
1.17 1.13 1.14 1.15 1.15
3.0999E+00 1.3387E+00 5.5922E−01 2.3261E−01 9.6777E−02
– 0.94 0.97 0.98 0.98
3.6210E+00 1.5357E+00 6.4293E−01 2.6818E−01 1.1185E−01
– 0.96 0.97 0.98 0.98
1.15 1.15 1.15 1.16
Mesh Type 2 Mesh Mesh Mesh Mesh Mesh
Level Level Level Level Level
1 2 3 4 5
Fig. 4.3. Example 2: Final refined mesh for k = 1 (a) Mesh Type 1 with #ELEM = 7977, DOF = 83254, and ∥∇ u0 − ∇ u∥ = 9.5310e − 3; (b) Mesh Type 2 with #ELEM = 10474, DOF = 90794 and ∥∇ u0 − ∇ u∥ = 9.398e − 3.
4.2. Example 2 In this test, we consider L-shape problem. Let Ω = (−1, 1)2 \(0, 1) × (−1, 0), κ = 1, and the exact solution is as follows: u = r 2/3 sin(2θ /3),
√
where r = x2 + y2 and θ = arctan(y/x). In this test, it is easy to check the right-hand function is f = 0. It is known that roughly speaking u ∈ H 2/3 (Ω ) and the singularity for this problem is at origin (0, 0). Because of the low regularity in this test, the uniform refinement cannot produce good numerical solutions with optimal rates in the simulation. We shall apply the adaptive finite element methods in solving this problem. A typical adaptive algorithm is described as follows: Solve → Estimate → Mark → Refine. In our numerical experiments, the following adaptive steps shall be performed: A. B. C. D. E.
(i)
(i)
We solve weak Galerkin numerical solution uh on a given polygonal mesh Th ; (i) Estimate the a posteriori error estimator η(i) and ηT ; (i) Mark the elements require refinement guided by the calculated error estimator ηT ; (i+1) Refine the marked elements, and then derive the updated polygonal mesh Th ; Repeat Steps A–D until the maximum iteration number is reached or a stopping criteria is satisfied.
Here, the Dörfler/bulk marking method with parameter θ = 0.2 will be used in the mark procedure. Weak Galerkin simulations with degrees k = 1, 2, 3 have been investigated. Assume u is the exact solution to (1.1)– (1.2) and uh = {u0 , ub } as the numerical solution to (2.4). In order to examine the performance of different weak Galerkin
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
421
Fig. 4.4. Example 2: Final refined mesh for k = 2 (a) Mesh Type 1 with #ELEM = 824, DOF = 10161, and ∥∇ u0 − ∇ u∥ = 9.133e − 3; (b) Mesh Type 2 with #ELEM = 825, DOF = 10263, and ∥∇ u0 − ∇ u∥ = 9.261e − 3.
Fig. 4.5. Example 2: Final refined mesh for k = 3 (a) Mesh Type 1 with #ELEM = 787, DOF = 14462, and ∥∇ u0 − ∇ u∥ = 9.159e − 3; (b) Mesh Type 2 with #ELEM = 794, DOF = 14652, and ∥∇ u0 − ∇ u∥ = 7.863e − 3.
elements, we denote
∥∇ u − ∇ u0 ∥ :=
(∑
∥∇ u − ∇ u0 ∥2T
)1/2 (4.2)
T ∈Th
and then measure the errors in this norm. Figs. 4.3–4.5 plot the final refined mesh for the similar accuracy ∥∇ u − ∇ u0 ∥ = 0.01 corresponding to different weak Galerkin elements. It can be seen from these figures that the estimators ηT can accurately detect the singularity and automatically add refinement focusing at (0, 0). Also, with the higher-order polynomial in the simulation, the error will converge faster and thus achieve the same accuracy within fewer adaptive steps than that of the lower-order polynomials. Furthermore, for the similar accuracy, higher order scheme requires fewer number in the element and better DOFs in the calculation than that of linear element. Fig. 4.6 plots the error profiles and convergence rates with respect to DOF for Mesh Type 1 and Mesh Type 2. All of them show that as we refine the mesh guided by the a posteriori error estimator ηT , the convergence rates with respect to WG element k = 1, 2, 3 are at O(DOF−0.5 ), O(DOF−1 ), and O(DOF−1.5 ), respectively. Moreover, the achieved convergence rate is the optimal rate in convergence.
422
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
Fig. 4.6. Example 2: Convergence test for Mesh Type 1 (left) and Mesh Type 2 (right).
Fig. 4.7. Example 3: Final refined mesh for k = 1 (a) Mesh Type 1 with #ELEM = 6723, DOF = 51471, and ∥∇ u0 − ∇ u∥ = 9.979e − 2; (b) Mesh Type 2 with #ELEM = 8386, DOF = 73486, and ∥∇ u0 − ∇ u∥ = 9.703e − 2.
4.3. Example 3 Next we will consider another numerical example of an elliptic interface problem from [23,36]. For Ω = (−1, 1)2 , let u(r , θ ) = r γ µ(θ ) be the solution to (1.1)–(1.2), and set
⎧ cos((π /2 − σ )γ ) · cos((θ − π/2 + ρ )γ ), ⎪ ⎪ ⎨ cos(ργ ) · cos((θ − π + σ )γ ), µ(θ ) = ⎪cos(σ γ ) · cos((θ − π − ρ )γ ), ⎪ ⎩ cos((π /2 − ρ )γ ) · cos((θ − 3π/2 − σ )γ ), with the numbers γ , ρ, σ satisfying the following relations
if 0 ≤ θ ≤ π/2 if π/2 ≤ θ ≤ π if π ≤ θ ≤ 3π/2 if 3π/2 ≤ θ ≤ 2π
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
423
Fig. 4.8. Example 3: Final refined mesh for k = 2 (a) Mesh Type 1 with #ELEM = 917, DOF = 11769, and ∥∇ u0 − ∇ u∥ = 9.998e − 2; (b) Mesh Type 2 with #ELEM = 1061, DOF = 14745, and ∥∇ u0 − ∇ u∥ = 9.587e − 2.
Fig. 4.9. Example 3: Final refined mesh for k = 3 (a) Mesh Type 1 with #ELEM = 611, DOF = 11502, and ∥∇ u0 − ∇ u∥ = 9.912e − 2; (b) Mesh Type 2 with #ELEM = 637, DOF = 12566, and ∥∇ u0 − ∇ u∥ = 9.992e − 2.
⎧ R := a1 /a2 = − tan((π/2 − σ )γ ) · cot(ργ ), ⎪ ⎪ ⎪ ⎪ ⎪ 1/R = − tan(ργ ) · cot(σ ρ ), ⎪ ⎪ ⎨ R = − tan(ργ ) · cot((π/2 − ρ )γ ), ⎪ 0 <γ <2 ⎪ ⎪ ⎪ ⎪ max {0, π γ − π } < 2γ ρ < min{πγ , π} ⎪ ⎪ ⎩ max{0, π − π γ } < −2γ σ < min{π, 2π − πγ }. We have tested the example in the worst-case scenarios, i.e., γ = 0.1, which produces a very singular solution u that is barely in H 1 . Here R = a1 /a2 ≈ 161.4476387975881, ρ = π/4, σ ≈ −14.92256510455152, and κ = a1 = R, in (0, 1)2 ∪ (−1, 0)2 and κ = a2 = 1 in Ω \((0, 1)2 ∪ (−1, 0)2 ). We start with 20 × 20 uniform rectangular mesh, and then employ the error estimator ηT to locate the trouble cells which need further refinement. Final refined meshes corresponding to k = 1, 2, 3 are plotted in Figs. 4.7–4.9 for similar
424
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
Fig. 4.10. Example 3: Error profiles and convergence results for weak Galerkin numerical schemes with element k = 1, 2, 3 for Mesh Type 1 (left) and Mesh Type 2 (right).
Fig. 4.11. Example 3: Error profiles and convergence results for weak Galerkin numerical schemes with element k = 1, 2, 3 for Mesh Type 2.
accuracy as ∥∇ u − ∇ u0 ∥ = 0.1. The results show that the refinement only concentrates at the origin and thus capture the singularity in this problem. One can observe that at the same accuracy, for the linear weak Galerkin element, one needs 4 or 5 times more DOF than that of weak Galerkin element k = 2. The saving in the computation cost of k = 3 is even more. That shows the advantage of using high-order numerical formulation linear schemes in the efficient adaptive numerical schemes. Figs. 4.10–4.11 plot the error profiles and convergence results for different weak Galerkin elements. All of them show the optimal rate in convergence, which are O(DOF−0.5 ), O(DOF−1 ), and O(DOF−1.5 ). Hence, these validate our theoretical conclusions. 5. Conclusions and future work In this paper, the a posteriori error estimate has been analyzed and applied to solve the second order elliptic interface problems. Theoretically, we proved the equivalence of the error estimator η and actual error measured in discrete H 1 norm. Numerically, several numerical experiments are conducted to validate our theoretical conclusions. It shows that η is efficient and reliable indicator to locate the singularity and thus guiding the local refinement for achieving optimal rate in convergence. Furthermore, two different mesh refinement strategies have been considered in the numerical examples. Both of them deal with the polygonal mesh locally very well. The refinement guided by indicator ηT will not propagate to the neighbor elements which do not need refinement.
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
425
In the future, our work will be extended to elliptic interface problems with non-body fitted mesh and used to design high-order numerical schemes by coupling advanced immersed finite element methods coupled with the approach developed in this paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36]
J. Wang, X. Ye, A weak Galerkin finite element method for second-order elliptic problems, J. Comput. Appl. Math. 241 (2013) 103–115. R.E. Bank, A.H. Sherman, A. Weiser, Refinement algorithms and data structures for regular local mesh refinement, Sci. Comput. (1983) 3–17. M.C. Rivara, Design and data structure for fully adaptive, multigrid finite element software, ACM Trans. Math. Software 10 (1984) 242–264. M.C. Rivara, Mesh refinement processes based on the generalized bisection of simplices, SIAM J. Numer. Anal. 21 (1984) 604–613. E.G. Sewell, Automatic Generation of Triangulations for Piecewise Polynomial Approximation (Ph.D. thesis), Purdue Univ., West Lafayette, Ind., 1972. L. Beirao da Veiga, K. Lipnikov, G. Manzini, Convergence analysis of the high-order mimetic finite difference method, Numer. Math. 113 (2009) 325–356. L. Beirao da Veiga, C. Lovadinab, D. Morac, A Virtual Element Method for elastic and inelastic problems on polytope meshes, Comput. Methods Appl. Mech. Eng. 1 (2015) 327–346. L. Beirao da Veiga, F. Brezzi, L. Marini, A. Russo, Mixed virtual element methods for general second order elliptic problems on polygonal meshes, ESAIM M2AN Math. Model. Numer. Anal. 50 (2016) 727–747. L. Beirao da Veiga, F. Brezzi, A. Cangiani, G. Manzini, D. Marini, A. Russo, Basic principles of virtual element methods, Math. Models Methods Appl. Sci. 23 (2013) 199–214. L. Beirao da Veiga, F. Brezzi, D. Marini, A. Russo, Virtual Element Method for general second order elliptic problems on polygonal meshes, Math. Models Methods Appl. Sci. 26 (2016) 729–750. L. Beirao da Veiga, K. Lipnikov, G. Manzini, The Mimetic Finite Difference Method for Elliptic Problems, in: MS & A. Modeling, Simulation and Applications, vol. 11, Springer, Cham, 2014. L. Beirao da Veiga, G. Manzini, A virtual element method with arbitrary regularity, IMA J. Numer. Anal. 34 (2014) 759–781. V. Gyrya1, K. Lipnikov1, G. Manzini, The arbitrary order mixed mimetic finite difference method for the diffusion equation, ESAIM M2AN Math. Model. Numer. Anal. 50 (2016) 851–877. D. Pietro, A. Ern, Hybrid high-order methods for variable-diffusion problems on general meshes, C. R. Math. 353 (2015) 31–34. D. Pietro, A. Ern, A hybrid high-order locking-free method for linear elasticity on general meshes, Comput. Methods Appl. Mech. Eng. 283 (2015) 1–21. D. Pietro, J. Droniou, A. Ern, A discontinuous-skeletal method for advection-diffusion-reaction on general meshes, SIAM J. Numer. Anal. 53 (2015) 2135–2157. D. Pietro, A. Ern, Arbitrary-order mixed methods for heterogeneous anisotropic diffusion on general meshes, IMA J. Numer. Anal. 37 (2017) 40–63. D. Pietro, J. Droniou, A hybrid high-order method for Leray-Lions elliptic equations on general meshes, Math. Comp. 307 (2017) 2159–2191. A. Cangiani, E. Georgoulis, T. Pryer, O. Sutton, A posteriori error estimates for the virtual element method, Numer. Math. (2017) 1–37. D. Pietro, R. Specogna, An a posteriori-driven adaptive mixed high-order method with application to electrostatics, J. Comput. Phys. 326 (2016) 35–55. M. Vohralk, S. Yousef, A simple a posteriori estimate on general polytopal meshes with applications to complex porous media flows, Comput. Methds Appl. Mech. Eng. 331 (2018) 728–760. C. Bi, V. Cinting, A posteriori error estimates of discontinuous Galerkin method for nonmonotone quasi-linear elliptic problems, J. Sci. Comput. 55 (2013) 659–687. Z. Cai, X. Ye, S. Zhang, Discontinuous Galerkin finite element methods for interface problems: a priori and a posteriori error estimations, SIAM J. Numer. Anal. 49 (2011) 1761–1787. L. Chen, J. Wang, X. Ye, A posteriori error estimates for weak Galerkin finite element methods for second order elliptic problem, J. Sci. Comput. 59 (2014) 496–511. O. Karakashian, F. Pascal, A posteriori error estimates for a discontinuous Galerkin approximation of second order elliptic problems, SIAM J. Numer. Anal. 45 (2003) 1777–1798. C. Lovadina, D. Marini, A posteriori error estimates for discontinuous Galerkin approximations of second order elliptic problems, J. Sci. Comput. 40 (2009) 340–359. D. Pietroa, E. Flauraudb, M. Vohralkc, S. Yousef, A posteriori error estimates, stopping criteria, and adaptivity for multiphase compositional Darcy flows in porous media, J. Comput. Phys. 276 (2014) 163–187. F. Wang, W. Han, J. Eichholze, X. Cheng, A posteriori error estimates for discontinuous Galerkin methods of obstacle problems, Nonlinear Anal. RWA 22 (2015) 664–679. J. Wang, Y. Wang, X. Ye, Unified a posteriori error estimator for finite element methods for the Stokes equations, Int. J. Numer. Anal. Model. 10 (2013) 551–570. H. Li, L. Mu, X. Ye, A posteriori error estimates for the weak Galerkin finite element methods on polytopal meshes, Commun. Comput. Phys. 26 (2019) 558–578. R.A. Adams, J.J.F. Fournier, Sobolev Spaces, Vol. 140, Academic Press, 2003. L. Mu, J. Wang, X. Ye, Weak Galerkin finite element methods on polytopal meshes, Int. J. Numer. Anal. Model. 12 (2015) 31–53. J. Wang, X. Ye, A Weak Galerkin mixed finite element method for second-order elliptic problems, Math. Comp. 83 (2014) 2101–2126. R. Duran, Mixed Finite Element Methods, Class Notes, 2007. M.J. Lai, G. Slavov, On recursive refinement of convex polygons, Comput. Aided Geom. Design 45 (2016) 83–90. R.B. Kellogg, On the Poisson equation with intersecting interfaces, Appl. Anal. 4 (1975) 101–129.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Weak Galerkin based a posteriori error estimates for second order elliptic interface problems on polygonal meshes Lin Mu 1 Computer Science and Mathematics Division Oak Ridge National Laboratory, Oak Ridge, TN, 37831, USA
article
info
Article history: Received 17 August 2017 Received in revised form 20 March 2019 MSC: primary 65N15 65N30 secondary 35J50
a b s t r a c t In this paper, we present a posteriori error estimate of weak Galerkin (WG) finite element methods based on the second order elliptic interface problems. This estimate can be applied to polygonal meshes or meshes with hanging nodes. The reliability and efficiency of the designed error estimator has been proved in this work. Extensive numerical tests are performed to validate our algorithm. These results demonstrate the effectiveness of the adaptive mesh refinement guided by the proposed error estimator. © 2019 Elsevier B.V. All rights reserved.
Keywords: Weak Galerkin Finite element methods Second-order elliptic interface problems A posterior error estimate Polygonal meshes
1. Introduction The weak Galerkin (WG) method was first proposed by Wang and Ye in [1] as a general framework for solving partial differential equations (PDE). This method is a natural extension of finite element methods for the discontinuous functions, with finite element formulation as simple as the weak form of the corresponding PDE. The key in WG methods is to replace the classical derivatives by the weakly defined discrete derivatives. Then a parameter-free stabilizer is added to produce a stable, symmetric, and positive definite formulation. Advantages of WG methods include the flexibility of employing polygonal meshes and naturally way to design high-order numerical schemes by using polynomials with higher degree. There are two different error analysis approaches: one is a priori error estimate and the other one is a posteriori error estimate. A priori error analysis offers theoretical conclusions about the convergence order of the errors measured in certain norms with respect to mesh size h. Due to the unknown exact solution is included in the measurement, the a priori error is not actually computable. The a posteriori error estimate extracts a computable quantity as estimator which is equivalent to the interested quantity measured in the chosen norm. If the solutions contain singularity, a posteriori error estimator can be used to drive mesh refinement and thus generate the adaptive finite element (AFE) methods. Recently, AFE methods are widely used in computational science and engineering to obtain better accuracy with minimum effort when singularities are included in the problems. The savings are achieved through only refine the mesh at the needed places containing bigger errors. In the adaptive mesh methods, several refinement rules have been developed, including red and green refinement [2], longest edge bisection [3,4] and newest vertex bisection [5]. Since the refinement may produce hanging nodes in the mesh, one needs further refine the adjoint neighbors for keeping the E-mail address: [email protected]. 1 This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research. This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. https://doi.org/10.1016/j.cam.2019.04.026 0377-0427/© 2019 Elsevier B.V. All rights reserved.
414
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
Fig. 2.1. Depiction of a shape-regular polygonal element ABCDEFA.
shape regular requirement. The limitations of quad/triangular mesh in two dimensional space and unnecessary neighbor element refinement in AFE motivate our project of investigate the a posteriori error estimate on the polygonal meshes. In recent research, many numerical schemes have been developed and analyzed on general polytopal meshes [6–18]. However it is still hard to find a posteriori estimators work for general polygonal/polyhedral meshes [19–21]. Even for the polygonal and polyhedral finite element methods, due to the technique difficulties in analysis, most a posteriori error analysis is still based on simplicial meshes and strong assumptions [22–29]. In reference [19], the authors designed a polygonal error estimator for the virtual element method. Reference [30] developed a simple posteriori error estimator for weak Galerkin finite element approximation for Poisson equations. In this paper, we shall extend the a posteriori error estimate to study the second order elliptic interface problems. In this work, our diffusion coefficient is assumed to be piecewise constant, which may be different values at different sub-regions. For simplicity, we consider a simple model problem that seeks an unknown function u satisfying
−∇ · κ (x)∇ u = f , in Ω , u = 0, on ∂ Ω ,
(1.1) (1.2)
where Ω is a polytopal domain in Rd (polygonal or polyhedral domain for d = 2, 3). The diffusion coefficient κ (x) is positive and piecewise constant on the subdomains of Ω with possible large jumps across interface. Denote κmin and κmax as the minimum and maximum values of κ . The rest of this paper is organized as follows: Weak Galerkin finite element methods are introduced in Section 2. Section 3 introduces a proposed a posteriori error estimator and is contributed to analyze reliability and efficiency of the proposed estimator. Three numerical examples are present in Section 4. The paper ends with conclusions and future work plan as Section 5. 2. Weak Galerkin finite element schemes 2.1. Finite element spaces We will adopt the standard notation and definitions for Sobolev spaces [31]. In particular, given an open region B ⊂ Rd , denote the Sobolev space by W s,r (B) on B, where s is a real number. When r = 2, W s,2 (B) = H s (B) denotes a Hilbert space equipped with the norm ∥ · ∥s,B . When s = 0, H 0 (B) coincides with the standard Lebesgue L2 (B) space equipped with the inner product (·, ·)T and the L2 -norm ∥ · ∥B . If B = Ω , the subscript B will be omitted. We shall use similar notations for ∂ B. Besides, ⟨·, ·⟩∂ B denotes the inner product on the edge ∂ B. The space H(div, B) is defined for vector-valued functions, which have square integrable divergence on B, and the norm is defined as follows:
∥v∥H(div,B) := (∥v∥2B + ∥∇ · v∥2B )1/2 . Let Th be a partition of the domain Ω consisting of elements which are closed and simply connected polygons in two dimension or polyhedra in three dimension; see Fig. 2.1. Let Eh be the set of all edges or flat faces in Th , and Eh0 = Eh \∂ Ω be the set of all interior edges or flat faces. Denote by hT the diameter for every element T ∈ Th and h = maxT ∈Th hT the mesh size for Th . We need some shape regularity assumptions for the partition Th described as below (cf. [32,33]).
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
415
A1: Assume that there exist two positive constants ϱv and ϱe such that for every element T ∈ Th we have
ϱe hde −1 ≤ |e|
ϱv hdT ≤ |T |,
(2.1)
for all edges or flat faces of T . A2: Assume that there exists a positive constant κ such that for every element T ∈ Th we have
κ hT ≤ he
(2.2)
for all edges or flat faces e of T . A3: Assume that the mesh edges or faces are flat. We further assume that for every T ∈ Th , and for every edge/face e ∈ ∂ T , there exists a pyramid P(e, T , Ae ) contained in T such that its base is identical with e, its apex is Ae ∈ T , and its height is proportional to hT with a proportionality constant σe bounded away from a fixed positive number σ ∗ from below. In other words, the height of the pyramid is given by σe hT such that σe ≥ σ ∗ > 0. The pyramid is also assumed to stand up above the base e in the sense that the angle between the vector xe − Ae , for any xe ∈ e, and the outward normal direction of e is strictly acute by falling into an interval [0, θ0 ] with θ0 < π2 . A4: Assume that each T ∈ Th has a circumscribed simplex S(T ) that is shape regular and has a diameter hS(T ) proportional to the diameter of T ; i.e., hS(T ) ≤ γ∗ hT with a constant γ∗ independent of T . Furthermore, assume that each circumscribed simplex S(T ) interests with only a fixed and small number of such simplices for all other elements T ∈ Th . Assume the interface for coefficients κ (x) does not cut element T ∈ Th . For simplicity of notation, in the following we shall assume d = 2. The results can be extended to three dimensional space. For a scalar-valued function v , we define the curl operator by
∇ × v = (−
∂v ∂v , ). ∂ x2 ∂ x1
For a given integer k ≥ 1, Pk (T ) denotes the space of polynomials of degree k on element T and Pk (e) denotes the space of polynomials of degree k on the edge e. Let Vh be the weak Galerkin finite element space associated with Th defined as follows Vh = {v = {v0 , vb } : v0 |T ∈ Pk (T ), vb |e ∈ Pk (e), e ∈ Eh , T ∈ Th } and Vh0 = {v : v ∈ Vh , vb = 0 on ∂ Ω }. We would like to emphasize that v0 is piecewise polynomial on each element T and vb is piecewise polynomial on each edge e. 2.2. WG finite element approximation For any weak function v = {v0 , vb }, a weak gradient ∇w v ∈ [Pk−1 (T )]2 is defined on T as the unique solution satisfying (∇w v, τ )T = −(v0 , ∇ · τ )T + ⟨vb , τ · n⟩∂ T ,
∀τ ∈ [Pk−1 (T )]2 .
(2.3)
Now the following introduces two bilinear forms on Vh : ∀v, w ∈ Vh s(v, w ) =
∑
sT (v, w ),
T ∈Th
(κ∇w v, ∇w w ) =
∑
(κT ∇w v, ∇w w )T ,
T ∈Th 1 where sT (v, w ) = h− T ⟨κT (v0 − vb ), w0 − wb ⟩∂ T and κT = κ|T .
Weak Galerkin Algorithm 1. Find uh ∈ Vh0 satisfying the following equations: (κ∇w uh , ∇w v ) + s(uh , v ) = (f , v0 ), where (f , v0 ) :=
∑
∫ T ∈Th
T
f v0 dT .
∀ v ∈ Vh0 ,
(2.4)
416
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
The energy norm can be defined as,
|||v|||2 = (κ∇w v, ∇w v ) + s(v, v ).
(2.5)
The existence and uniqueness of numerical solution can be found in [32]. Furthermore, we cite the following theorem regarding a priori error estimate [32]. Theorem 2.1. Let uh ∈ Vh and u ∈ H k+1 (Ω ) be the weak Galerkin finite element solution and exact solution to problems (2.4) and (1.1)–(1.2). Then, there exists a constant C independent of mesh size h such that
|||u − uh ||| ≤ Chk ∥u∥k+1 , ∥u − uh ∥ ≤ Chk+1 ∥u∥k+1 .
(2.6) (2.7)
3. A posteriori error estimate Let fh be the L2 projection of f to Vh . Then a local estimator is defined as follows
ηT2 = sT (uh , uh ) + osc2 (f , T ),
(3.1)
where osc(f , T ) is a high order local data oscillation defined by osc2 (f , T ) = h2T ∥f − fh ∥2T . Define a global error estimator and data oscillation as
η2 =
∑
ηT2 ,
T ∈Th
osc(f , Th ) = 2
∑
osc(f , T )2 .
T ∈Th
Let T be an element with e as an edge. The following trace inequality will be also used in the following section. It is well known that there exists a constant C such that for any function g ∈ H 1 (T )
( 1 ) 2 2 ∥g ∥2e ≤ C h− T ∥g ∥T + hT ∥∇ g ∥T .
(3.2)
For each element T ∈ Th , denote by Q0 the L2 projection from L2 (T ) to Pk (T ), and by Qb the L2 projection from L2 (e) to Pk−1 (e). By an abuse of notation, we use here the same symbols for Q0 and Qb as global operators. Denote by Qh the L2 projection from [L2 (T )]2 to a local weak gradient space [Pk−1 (T )]2 . Define Qh u = {Q0 u, Qb u} ∈ Vh . The following lemma can be found in [32] Lemma 3.1. For any element T ∈ Th , we have the following commutative property for φ ∈ H 1 (T ),
∇w (Qh φ ) = Qh (∇φ ), ∇w φ = Qh (∇φ ).
(3.3) (3.4)
Theorem 3.2 (Reliability). Let uh ∈ Vh0 and u ∈ H01 (Ω ) be the solutions of (2.4) and (1.1)–(1.2) respectively. Then there exists a positive constant C independent of mesh size h such that:
|||u − uh ||| ≤ C η.
(3.5)
Proof. We shall apply Helmholtz decomposition first. It is well known [34] that for κ (∇w u − ∇w uh ) ∈ [L2 (Ω )]2 , there exist ψ ∈ H01 (Ω ) and φ ∈ H 1 (Ω ) such that
κ (∇w u − ∇w uh ) = κ∇ψ + ∇ × φ
(3.6)
and that
∥κ 1/2 (∇w u − ∇w uh )∥2 = ∥κ 1/2 ∇ψ∥2 + ∥κ −1/2 ∇ × φ∥2 . (3.4) and the definition of Qh implies: ∀q ∈ [Pk−1 (T )]2 , (κT ∇w u, q)T = (κT Qh (∇ u), q)T = (κT ∇ u, q)T .
(3.7)
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
417
From the above equality, and (3.6), one obtains,
∥κ 1/2 (∇w u − ∇w uh )∥2 =
∑
( κ T ( ∇ u − ∇ w uh ) , ∇ w u − ∇ w uh ) T
T ∈Th
=
∑
(∇ u − ∇w uh , κT (∇w u − ∇w uh ))T
(3.8)
T ∈Th
= (κ∇ψ, (∇ u − ∇w uh )) + (∇ × φ, ∇ u − ∇w uh ). By (3.3), the definition of Qh , and (f , Q0 ψ ) = (κ∇w uh , ∇w Qh ψ ) + s(uh , Qh ψ )
= (κ∇w uh , Qh (∇ψ )) + s(uh , Qh ψ ) = (κ∇w uh , ∇ψ ) + s(uh , Qh ψ )
(3.9)
it is derived that (κ∇w uh , ∇ψ ) = (f , Q0 ψ ) − s(uh , Qh ψ ).
(3.10)
Also, by the definition of Qb , we have 1/2
1 −1 2 h− T ∥κT (Q0 ψ − Qb ψ )∥∂ T = hT ⟨κT (Q0 ψ − Qb ψ ), Q0 ψ − Qb ψ⟩∂ T 1 = h− T ⟨κT (Q0 ψ − ψ ), Q0 ψ − Qb ψ⟩∂ T ( )1/2 ( )1/2 1/2 1/2 −1 −1 2 2 ≤ hT ∥κT (Q0 ψ − ψ )∥∂ T hT ∥κT (Q0 ψ − Qb )ψ∥∂ T
and thus by trace inequality and property of projection operator [33], 1/2
1/2
h−1 ∥κT (Q0 ψ − Qb ψ )∥2 ∂ T ≤ h−1 ∥κT (Q0 ψ − ψ )∥2∂ T 1/2
≤ C ∥κT ∇ψ∥2T .
(3.11)
Using (3.9), Cauchy–Schwarz inequality, inequality (3.11), and (3.7), we have (κ∇ψ, (∇ u − ∇w uh )) = (κ (∇ u − ∇w uh ), ∇ψ )
= (κ∇ u, ∇ψ ) − (κ∇w uh , ∇ψ ) = (f , ψ ) − (f , Q0 ψ ) + s(uh , Qh ψ ) ∑ = (f − fh , ψ − Q0 ψ ) + hT −1 ⟨κT (u0 − ub ), Q0 ψ − Qb ψ⟩∂ T T ∈Th
≤ C (osc(f , Th ) + s1/2 (uh , uh ))∥κ 1/2 ∇ψ∥ ≤ C (osc(f , Th ) + s1/2 (uh , uh ))∥κ 1/2 (∇w u − ∇w uh )∥. Since ∇ × φ ∈ H(div, Ω ) and u ∈
H01 (Ω ),
(3.12)
we have
( ∇ u − ∇ w uh , ∇ × φ ) = ( ∇ u, ∇ × φ ) − ( ∇ w uh , ∇ × φ )
= − ( ∇ w uh , ∇ × φ ) . It follows from (2.3), uh ∈ Vh0 and the integration by parts, (∇w uh , ∇ × φ )T = (∇w uh , Qh ∇ × φ )T
= = = =
⟨ub , Qh (∇ × φ ) · n⟩∂ T − (u0 , ∇ · (Qh ∇ × φ ))T ⟨ub , Qh (∇ × φ ) · n⟩∂ T − ⟨u0 , (Qh ∇ × φ ) · n⟩∂ T + (∇ u0 , Qh ∇ × φ )T ⟨ub − u0 , Qh (∇ × φ ) · n⟩∂ T + (∇ u0 , ∇ × φ )T ⟨ub − u0 , Qh (∇ × φ ) · n⟩∂ T + ⟨u0 , (∇ × φ ) · n⟩T
and ub is single valued on each edge e ∈ Eh , one have
∑ T ∈Th
( ∇ w uh , ∇ × φ ) T =
∑
(⟨ub − u0 , Qh (∇ × φ ) · n⟩∂ T + ⟨u0 − ub , ∇ × φ·n⟩∂ T ).
T ∈Th
The Cauchy–Schwarz inequality, the trace inequality (3.2) and the inverse inequality imply
⟨u0 − ub , Qh ∇ × φ · n⟩∂ T ≤ C ∥u0 − ub ∥∂ T ∥Qh ∇ × φ · n∥∂ T 1/2 −1/2 ≤ C sT (uh , uh )∥κT ∇ × φ∥T .
(3.13)
418
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
Using the inverse inequality and the fact ∇ · ∇ × φ = 0, we arrive at
∑
⟨ u0 − ub , ∇ × φ · n ⟩ ∂ T =
⟨u0 − ub , ∇ × φ · n⟩e
e⊂∂ T
∑
≤
1/2
−1/2
∥κT (u0 − ub )∥H 1/2 (e) ∥κT
∇ × φ · n∥H −1/2 (e)
e⊂∂ T
≤ C(
∑
−1/2
hT
1/2
−1/2
∥κT (u0 − ub )∥2e )1/2 ∥κT
∇ × φ∥H(div,T )
e⊂∂ T 1/2
−1/2
≤ C (sT (uh , uh ))∥κT
∇ × φ∥T .
Combining the above three estimates and taking a summation over T give (∇w uh , ∇ × φ ) ≤ C (s1/2 (uh , uh ))∥κ −1/2 ∇ × φ∥.
(3.14)
Using (3.8), (3.12), (3.13) (3.14) and (3.7), we have
∥k1/2 (∇w u − ∇w uh )∥ ≤ C η. Next we easily have s(u − uh , u − uh ) = s(uh , uh ) ≤ η2 . Thus we complete the proof.
□
Define the localized error on each element T as follows,
|||v|||2T = (∇w v, ∇w v )T + sT (v, v ). Since sT (v, v ) is only part of the local error |||v|||2T , then the following local lower bound is automatically obtained. Theorem 3.3 (Efficiency). The local estimator ηT as defined in (3.1) is bounded as the following inequality:
ηT2 ≤ |||u − uh |||2T + osc2 (f , T ).
(3.15)
4. Numerical example In this section, we shall validate the proposed a posteriori error estimate for several tests. First, we shall explore the convergence properties of errors measured in ||| · |||-norm (denoted by H 1 -error in tables and H 1 err in figures) and the estimator η. The effectivity of the estimator is defined as follows, Eff-index =
η . |||Qh u − uh |||
(4.1)
The polygon refinement will be generated by two different strategies:
• Mesh Type 1: connect the barycenter of each element to the mid-point of each edge corresponding to this element. This approach was also used in [19]. In this way, all the polygons will be refined to several quadrilaterals. An example can be found as Fig. 4.1(a). • Mesh Type 2: connect the barycenter of each element to the mid-point of each edge corresponding to this element, and mark the mid-points of the connecting edges. By connecting these marked mid-points, marked mid-points and the middle point of each edge, one can recursively refine the mesh. This approach was developed in [35]. An example can be found as Fig. 4.1(b). 4.1. Example 1 Let Ω = (0, 1)2 , assume the interface at x = 1/2, and the diffusion coefficient
κ (x) =
{
1, if x < 1/2, 2, if x ≥ 1/2,
and u = sin(2π x)2 cos(π y).
In this test, we shall check the effectivity index as defined in (4.1) for comparing our posteriori error estimator η and error |||Qh u − uh |||. For simplicity, we only consider linear WG element on the polygonal mesh. We start the WG simulation on the initial mesh as shown in Fig. 4.2(a), which is denoted as Mesh Level 1. Two different uniformly refinement techniques (Mesh Type 1 and Mesh Type 2) will be applied to generate a sequence of meshes. Two examples of Mesh Level 2 of two different refinement strategies can be found as Fig. 4.2(b)–(c). The error profiles, convergence results, and Eff-index are reported in Table 4.1. Because of the setting for exact solution in this problem, one can observe the errors, measured in ||| · |||-norm and estimator η, converge at the optimal rate in convergence at order O(h) with respect to mesh size h. Also, the last column presents the Eff-Index, which is shown as a constant, and thus confirms the equivalency of η and |||Qh u − uh |||.
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
419
Fig. 4.1. Two different refinements (a) Mesh Type 1; (b) Mesh Type 2.
Fig. 4.2. Example 1: (a) Mesh Level 1 for initial Mesh with #ELEM = 58, DOF = 486; (b) Mesh Level 2 of Mesh Type 1 with #ELEM = 288, DOF = 3528; (c) Mesh Level 2 of Mesh Type 2 with #ELEM = 346, DOF = 2814.
420
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425 Table 4.1 Example 1. Error profiles and effectivity index. Mesh Type 1 Mesh Mesh Mesh Mesh Mesh Mesh
Level Level Level Level Level
1 2 3 4 5
H 1 -error
Order
η
Order
Eff-Index
3.0999E+00 1.4237E+00 7.3099E−01 3.6861E−01 1.8475E−01
– 1.12 0.96 0.99 1.00
3.6210E+00 1.6143E+00 8.3385E−01 4.2212E−01 2.1191E−01
– 1.17 0.95 0.98 0.99
1.17 1.13 1.14 1.15 1.15
3.0999E+00 1.3387E+00 5.5922E−01 2.3261E−01 9.6777E−02
– 0.94 0.97 0.98 0.98
3.6210E+00 1.5357E+00 6.4293E−01 2.6818E−01 1.1185E−01
– 0.96 0.97 0.98 0.98
1.15 1.15 1.15 1.16
Mesh Type 2 Mesh Mesh Mesh Mesh Mesh
Level Level Level Level Level
1 2 3 4 5
Fig. 4.3. Example 2: Final refined mesh for k = 1 (a) Mesh Type 1 with #ELEM = 7977, DOF = 83254, and ∥∇ u0 − ∇ u∥ = 9.5310e − 3; (b) Mesh Type 2 with #ELEM = 10474, DOF = 90794 and ∥∇ u0 − ∇ u∥ = 9.398e − 3.
4.2. Example 2 In this test, we consider L-shape problem. Let Ω = (−1, 1)2 \(0, 1) × (−1, 0), κ = 1, and the exact solution is as follows: u = r 2/3 sin(2θ /3),
√
where r = x2 + y2 and θ = arctan(y/x). In this test, it is easy to check the right-hand function is f = 0. It is known that roughly speaking u ∈ H 2/3 (Ω ) and the singularity for this problem is at origin (0, 0). Because of the low regularity in this test, the uniform refinement cannot produce good numerical solutions with optimal rates in the simulation. We shall apply the adaptive finite element methods in solving this problem. A typical adaptive algorithm is described as follows: Solve → Estimate → Mark → Refine. In our numerical experiments, the following adaptive steps shall be performed: A. B. C. D. E.
(i)
(i)
We solve weak Galerkin numerical solution uh on a given polygonal mesh Th ; (i) Estimate the a posteriori error estimator η(i) and ηT ; (i) Mark the elements require refinement guided by the calculated error estimator ηT ; (i+1) Refine the marked elements, and then derive the updated polygonal mesh Th ; Repeat Steps A–D until the maximum iteration number is reached or a stopping criteria is satisfied.
Here, the Dörfler/bulk marking method with parameter θ = 0.2 will be used in the mark procedure. Weak Galerkin simulations with degrees k = 1, 2, 3 have been investigated. Assume u is the exact solution to (1.1)– (1.2) and uh = {u0 , ub } as the numerical solution to (2.4). In order to examine the performance of different weak Galerkin
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
421
Fig. 4.4. Example 2: Final refined mesh for k = 2 (a) Mesh Type 1 with #ELEM = 824, DOF = 10161, and ∥∇ u0 − ∇ u∥ = 9.133e − 3; (b) Mesh Type 2 with #ELEM = 825, DOF = 10263, and ∥∇ u0 − ∇ u∥ = 9.261e − 3.
Fig. 4.5. Example 2: Final refined mesh for k = 3 (a) Mesh Type 1 with #ELEM = 787, DOF = 14462, and ∥∇ u0 − ∇ u∥ = 9.159e − 3; (b) Mesh Type 2 with #ELEM = 794, DOF = 14652, and ∥∇ u0 − ∇ u∥ = 7.863e − 3.
elements, we denote
∥∇ u − ∇ u0 ∥ :=
(∑
∥∇ u − ∇ u0 ∥2T
)1/2 (4.2)
T ∈Th
and then measure the errors in this norm. Figs. 4.3–4.5 plot the final refined mesh for the similar accuracy ∥∇ u − ∇ u0 ∥ = 0.01 corresponding to different weak Galerkin elements. It can be seen from these figures that the estimators ηT can accurately detect the singularity and automatically add refinement focusing at (0, 0). Also, with the higher-order polynomial in the simulation, the error will converge faster and thus achieve the same accuracy within fewer adaptive steps than that of the lower-order polynomials. Furthermore, for the similar accuracy, higher order scheme requires fewer number in the element and better DOFs in the calculation than that of linear element. Fig. 4.6 plots the error profiles and convergence rates with respect to DOF for Mesh Type 1 and Mesh Type 2. All of them show that as we refine the mesh guided by the a posteriori error estimator ηT , the convergence rates with respect to WG element k = 1, 2, 3 are at O(DOF−0.5 ), O(DOF−1 ), and O(DOF−1.5 ), respectively. Moreover, the achieved convergence rate is the optimal rate in convergence.
422
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
Fig. 4.6. Example 2: Convergence test for Mesh Type 1 (left) and Mesh Type 2 (right).
Fig. 4.7. Example 3: Final refined mesh for k = 1 (a) Mesh Type 1 with #ELEM = 6723, DOF = 51471, and ∥∇ u0 − ∇ u∥ = 9.979e − 2; (b) Mesh Type 2 with #ELEM = 8386, DOF = 73486, and ∥∇ u0 − ∇ u∥ = 9.703e − 2.
4.3. Example 3 Next we will consider another numerical example of an elliptic interface problem from [23,36]. For Ω = (−1, 1)2 , let u(r , θ ) = r γ µ(θ ) be the solution to (1.1)–(1.2), and set
⎧ cos((π /2 − σ )γ ) · cos((θ − π/2 + ρ )γ ), ⎪ ⎪ ⎨ cos(ργ ) · cos((θ − π + σ )γ ), µ(θ ) = ⎪cos(σ γ ) · cos((θ − π − ρ )γ ), ⎪ ⎩ cos((π /2 − ρ )γ ) · cos((θ − 3π/2 − σ )γ ), with the numbers γ , ρ, σ satisfying the following relations
if 0 ≤ θ ≤ π/2 if π/2 ≤ θ ≤ π if π ≤ θ ≤ 3π/2 if 3π/2 ≤ θ ≤ 2π
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
423
Fig. 4.8. Example 3: Final refined mesh for k = 2 (a) Mesh Type 1 with #ELEM = 917, DOF = 11769, and ∥∇ u0 − ∇ u∥ = 9.998e − 2; (b) Mesh Type 2 with #ELEM = 1061, DOF = 14745, and ∥∇ u0 − ∇ u∥ = 9.587e − 2.
Fig. 4.9. Example 3: Final refined mesh for k = 3 (a) Mesh Type 1 with #ELEM = 611, DOF = 11502, and ∥∇ u0 − ∇ u∥ = 9.912e − 2; (b) Mesh Type 2 with #ELEM = 637, DOF = 12566, and ∥∇ u0 − ∇ u∥ = 9.992e − 2.
⎧ R := a1 /a2 = − tan((π/2 − σ )γ ) · cot(ργ ), ⎪ ⎪ ⎪ ⎪ ⎪ 1/R = − tan(ργ ) · cot(σ ρ ), ⎪ ⎪ ⎨ R = − tan(ργ ) · cot((π/2 − ρ )γ ), ⎪ 0 <γ <2 ⎪ ⎪ ⎪ ⎪ max {0, π γ − π } < 2γ ρ < min{πγ , π} ⎪ ⎪ ⎩ max{0, π − π γ } < −2γ σ < min{π, 2π − πγ }. We have tested the example in the worst-case scenarios, i.e., γ = 0.1, which produces a very singular solution u that is barely in H 1 . Here R = a1 /a2 ≈ 161.4476387975881, ρ = π/4, σ ≈ −14.92256510455152, and κ = a1 = R, in (0, 1)2 ∪ (−1, 0)2 and κ = a2 = 1 in Ω \((0, 1)2 ∪ (−1, 0)2 ). We start with 20 × 20 uniform rectangular mesh, and then employ the error estimator ηT to locate the trouble cells which need further refinement. Final refined meshes corresponding to k = 1, 2, 3 are plotted in Figs. 4.7–4.9 for similar
424
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
Fig. 4.10. Example 3: Error profiles and convergence results for weak Galerkin numerical schemes with element k = 1, 2, 3 for Mesh Type 1 (left) and Mesh Type 2 (right).
Fig. 4.11. Example 3: Error profiles and convergence results for weak Galerkin numerical schemes with element k = 1, 2, 3 for Mesh Type 2.
accuracy as ∥∇ u − ∇ u0 ∥ = 0.1. The results show that the refinement only concentrates at the origin and thus capture the singularity in this problem. One can observe that at the same accuracy, for the linear weak Galerkin element, one needs 4 or 5 times more DOF than that of weak Galerkin element k = 2. The saving in the computation cost of k = 3 is even more. That shows the advantage of using high-order numerical formulation linear schemes in the efficient adaptive numerical schemes. Figs. 4.10–4.11 plot the error profiles and convergence results for different weak Galerkin elements. All of them show the optimal rate in convergence, which are O(DOF−0.5 ), O(DOF−1 ), and O(DOF−1.5 ). Hence, these validate our theoretical conclusions. 5. Conclusions and future work In this paper, the a posteriori error estimate has been analyzed and applied to solve the second order elliptic interface problems. Theoretically, we proved the equivalence of the error estimator η and actual error measured in discrete H 1 norm. Numerically, several numerical experiments are conducted to validate our theoretical conclusions. It shows that η is efficient and reliable indicator to locate the singularity and thus guiding the local refinement for achieving optimal rate in convergence. Furthermore, two different mesh refinement strategies have been considered in the numerical examples. Both of them deal with the polygonal mesh locally very well. The refinement guided by indicator ηT will not propagate to the neighbor elements which do not need refinement.
L. Mu / Journal of Computational and Applied Mathematics 361 (2019) 413–425
425
In the future, our work will be extended to elliptic interface problems with non-body fitted mesh and used to design high-order numerical schemes by coupling advanced immersed finite element methods coupled with the approach developed in this paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36]
J. Wang, X. Ye, A weak Galerkin finite element method for second-order elliptic problems, J. Comput. Appl. Math. 241 (2013) 103–115. R.E. Bank, A.H. Sherman, A. Weiser, Refinement algorithms and data structures for regular local mesh refinement, Sci. Comput. (1983) 3–17. M.C. Rivara, Design and data structure for fully adaptive, multigrid finite element software, ACM Trans. Math. Software 10 (1984) 242–264. M.C. Rivara, Mesh refinement processes based on the generalized bisection of simplices, SIAM J. Numer. Anal. 21 (1984) 604–613. E.G. Sewell, Automatic Generation of Triangulations for Piecewise Polynomial Approximation (Ph.D. thesis), Purdue Univ., West Lafayette, Ind., 1972. L. Beirao da Veiga, K. Lipnikov, G. Manzini, Convergence analysis of the high-order mimetic finite difference method, Numer. Math. 113 (2009) 325–356. L. Beirao da Veiga, C. Lovadinab, D. Morac, A Virtual Element Method for elastic and inelastic problems on polytope meshes, Comput. Methods Appl. Mech. Eng. 1 (2015) 327–346. L. Beirao da Veiga, F. Brezzi, L. Marini, A. Russo, Mixed virtual element methods for general second order elliptic problems on polygonal meshes, ESAIM M2AN Math. Model. Numer. Anal. 50 (2016) 727–747. L. Beirao da Veiga, F. Brezzi, A. Cangiani, G. Manzini, D. Marini, A. Russo, Basic principles of virtual element methods, Math. Models Methods Appl. Sci. 23 (2013) 199–214. L. Beirao da Veiga, F. Brezzi, D. Marini, A. Russo, Virtual Element Method for general second order elliptic problems on polygonal meshes, Math. Models Methods Appl. Sci. 26 (2016) 729–750. L. Beirao da Veiga, K. Lipnikov, G. Manzini, The Mimetic Finite Difference Method for Elliptic Problems, in: MS & A. Modeling, Simulation and Applications, vol. 11, Springer, Cham, 2014. L. Beirao da Veiga, G. Manzini, A virtual element method with arbitrary regularity, IMA J. Numer. Anal. 34 (2014) 759–781. V. Gyrya1, K. Lipnikov1, G. Manzini, The arbitrary order mixed mimetic finite difference method for the diffusion equation, ESAIM M2AN Math. Model. Numer. Anal. 50 (2016) 851–877. D. Pietro, A. Ern, Hybrid high-order methods for variable-diffusion problems on general meshes, C. R. Math. 353 (2015) 31–34. D. Pietro, A. Ern, A hybrid high-order locking-free method for linear elasticity on general meshes, Comput. Methods Appl. Mech. Eng. 283 (2015) 1–21. D. Pietro, J. Droniou, A. Ern, A discontinuous-skeletal method for advection-diffusion-reaction on general meshes, SIAM J. Numer. Anal. 53 (2015) 2135–2157. D. Pietro, A. Ern, Arbitrary-order mixed methods for heterogeneous anisotropic diffusion on general meshes, IMA J. Numer. Anal. 37 (2017) 40–63. D. Pietro, J. Droniou, A hybrid high-order method for Leray-Lions elliptic equations on general meshes, Math. Comp. 307 (2017) 2159–2191. A. Cangiani, E. Georgoulis, T. Pryer, O. Sutton, A posteriori error estimates for the virtual element method, Numer. Math. (2017) 1–37. D. Pietro, R. Specogna, An a posteriori-driven adaptive mixed high-order method with application to electrostatics, J. Comput. Phys. 326 (2016) 35–55. M. Vohralk, S. Yousef, A simple a posteriori estimate on general polytopal meshes with applications to complex porous media flows, Comput. Methds Appl. Mech. Eng. 331 (2018) 728–760. C. Bi, V. Cinting, A posteriori error estimates of discontinuous Galerkin method for nonmonotone quasi-linear elliptic problems, J. Sci. Comput. 55 (2013) 659–687. Z. Cai, X. Ye, S. Zhang, Discontinuous Galerkin finite element methods for interface problems: a priori and a posteriori error estimations, SIAM J. Numer. Anal. 49 (2011) 1761–1787. L. Chen, J. Wang, X. Ye, A posteriori error estimates for weak Galerkin finite element methods for second order elliptic problem, J. Sci. Comput. 59 (2014) 496–511. O. Karakashian, F. Pascal, A posteriori error estimates for a discontinuous Galerkin approximation of second order elliptic problems, SIAM J. Numer. Anal. 45 (2003) 1777–1798. C. Lovadina, D. Marini, A posteriori error estimates for discontinuous Galerkin approximations of second order elliptic problems, J. Sci. Comput. 40 (2009) 340–359. D. Pietroa, E. Flauraudb, M. Vohralkc, S. Yousef, A posteriori error estimates, stopping criteria, and adaptivity for multiphase compositional Darcy flows in porous media, J. Comput. Phys. 276 (2014) 163–187. F. Wang, W. Han, J. Eichholze, X. Cheng, A posteriori error estimates for discontinuous Galerkin methods of obstacle problems, Nonlinear Anal. RWA 22 (2015) 664–679. J. Wang, Y. Wang, X. Ye, Unified a posteriori error estimator for finite element methods for the Stokes equations, Int. J. Numer. Anal. Model. 10 (2013) 551–570. H. Li, L. Mu, X. Ye, A posteriori error estimates for the weak Galerkin finite element methods on polytopal meshes, Commun. Comput. Phys. 26 (2019) 558–578. R.A. Adams, J.J.F. Fournier, Sobolev Spaces, Vol. 140, Academic Press, 2003. L. Mu, J. Wang, X. Ye, Weak Galerkin finite element methods on polytopal meshes, Int. J. Numer. Anal. Model. 12 (2015) 31–53. J. Wang, X. Ye, A Weak Galerkin mixed finite element method for second-order elliptic problems, Math. Comp. 83 (2014) 2101–2126. R. Duran, Mixed Finite Element Methods, Class Notes, 2007. M.J. Lai, G. Slavov, On recursive refinement of convex polygons, Comput. Aided Geom. Design 45 (2016) 83–90. R.B. Kellogg, On the Poisson equation with intersecting interfaces, Appl. Anal. 4 (1975) 101–129.