Novel Adaptive Hybrid Discontinuous Galerkin Algorithms for Elliptic Problems

Dongwook Shin; Yoongu Hwang; Eun-Jae Park

doi:10.1515/cmam-2020-0114

Publicly Available Published by De Gruyter June 29, 2021

Novel Adaptive Hybrid Discontinuous Galerkin Algorithms for Elliptic Problems

Dongwook Shin , Yoongu Hwang and Eun-Jae Park

From the journal Computational Methods in Applied Mathematics

https://doi.org/10.1515/cmam-2020-0114

Abstract

In this work, we develop novel adaptive hybrid discontinuous Galerkin algorithms for second-order elliptic problems. For this, two types of reliable and efficient, modulo a data-oscillation term, and fully computable a posteriori error estimators are developed: the first one is a simple residual type error estimator, and the second is a flux reconstruction based error estimator of a guaranteed type for polynomial approximations of any degree by using a simple postprocessing. These estimators can achieve high-order accuracy for both smooth and nonsmooth problems even with high-order approximations. In order to enhance the performance of adaptive algorithms, we introduce 𝐾-means clustering based marking strategy. The choice of marking parameter is crucial in the performance of the existing strategy such as maximum and bulk criteria; however, the optimal choice is not known. The new strategy has no unknown parameter. Several numerical examples are given to illustrate the performance of the new marking strategy along with our estimators.

Keywords: A Posteriori Error Estimation; Adaptive Finite Element Method; Hybrid Discontinuous Galerkin; Error Estimates; Marking Strategy

MSC 2010: 65N12; 65N30

1 Introduction

Various natural and physical phenomena can be described by partial differential equations. In order to solve these equations, a number of numerical methods have been introduced and developed. In particular, discontinuous Galerkin (DG) and mixed finite element methods are locally conservative and important for a certain class of problems, e.g., Darcy flows in porous media and convection diffusion equations. These methods allow to achieve high-order accuracy with high-order approximations. However, convergence rates are influenced by the regularity of the problem.

An adaptive algorithm aims to overcome slow convergence for nonsmooth solutions or in the presence of boundary and internal layers. The algorithm searches elements whose error is large, and it produces an adaptive mesh to selectively refine these elements. In this procedure, the error decreases quickly in terms of degrees of freedom. In order to design adaptive algorithms, many researchers have studied a posteriori error estimates and introduced a posteriori error estimators for various problems [43, 1, 3, 38, 8, 11, 10, 13, 14, 15, 29, 30, 31]. In the conforming finite element method, guaranteed upper bounds can be obtained using the concept of Prager and Synge [37, 42]. And we can find a general functional framework and the idea of using a local residual equilibration procedure in [38] and [1], respectively. In the nonconforming methods, a posteriori error estimates have been derived for Crouzeix–Raviart, mixed and DG methods (see [18, 2, 33, 32, 12, 34, 47]). On the other hand, many researchers make an effort to construct fully computable a posteriori error estimators which are completely free of unknown constants, lead to guaranteed upper bounds (reliability) and represent local lower bounds (efficiency) of the error [9, 2, 44, 16, 46].

A number of hybrid DG (HDG) methods have been introduced and developed by many researchers in order to reduce the total degrees of freedom via static condensation. In particular, Cockburn and his coworkers suggested HDG methods [17] via mixed formulation. They help achieve superconvergence through postprocessing in the sense of Stenberg [41] with the additional degrees of freedom. On the other hand, Jeon and Park introduced low-order cell boundary element methods [24] and then proposed generalized cell boundary element (GCBE) methods [25, 26, 27] which can be viewed as an HDG method based on the primal formulation without introducing the flux variable. The resulting method is particularly simple, and it needs neither stabilization nor penalty terms. This method also enjoys a local mass conservation property and continuity of the average flux on each interior edge for even-degree polynomial approximations. The number of degrees of freedom of this method is much less than standard conforming FE method’s with high-order ( k ≥ 3 ) approximations because the numerical trace is continuous and degrees of freedom are located only on the skeleton. The interface stabilized FE method proposed by Wells [45] is a symmetric interior penalty HDG (S-HDG) method, which also reduces lots of degrees of freedom with the continuous trace.

In this paper, we introduce a unified variational formulation which covers the GCBE, incomplete interior penalty HDG (I-HDG), non-symmetric interior penalty HDG (N-HDG) and S-HDG methods. The HDG method under consideration is classified as the HDG_C and HDG_D methods by choosing the continuous trace and discontinuous trace, respectively. Then we construct two types of a posteriori error estimators: a residual type estimator and a flux reconstruction based estimator. For the flux reconstruction based estimator, we exploit the flux reconstruction via the Raviart–Thomas–Nédélec space. The residual type estimator, on the other hand, does not need the flux reconstruction. These estimators are fully computable, reliable and efficient modulo a data-oscillation term. Proofs of reliability and efficiency are inspired by the work of Ern and Vohralík [20].

It is well known that the performance of an adaptive algorithm is affected by the choice of the marking strategy. In this work, we introduce a new marking strategy based on the unsupervised machine learning, particularly, 𝐾-means clustering. We compare the performance of the new marking strategy with the existing ones such as average, maximum and bulk criteria. It should be noted that the performance of maximum and bulk criteria depends on the choice of the marking parameter; however, the optimal choice is not known. The new strategy has no unknown marking parameter. While inappropriate choices of 𝐾, the number of clusters, yield poor results in general applications of 𝐾-means clustering, the adaptivity works well with 𝐾-means clustering with K = 2 . Indeed, our numerical experiments show that the proposed strategy outperforms the existing marking strategies for low-order polynomials. For high-order cases, the proposed strategy still shows better performance than the existing marking criteria in general. In the case when the optimal marking parameter is known, the existing one is marginally better. Thus, our new strategy can be very useful in practice. To the best of the authors’ knowledge, no machine learning technique applied to the marking strategy has been presented in the literature so far.

The remainder of the paper is organized as follows. In the next section, we briefly introduce HDG methods and state their properties. In Section 3, we propose and analyze two types of a posteriori error estimates. Several numerical experiments are presented to test the performance of error estimators in Section 4. We compare the residual type with the flux reconstruction based estimators and low- with high-order approximations. Concluding remarks are given in the last section.

2 Methods

In order to simplify our discussion, let us consider the following model problem with homogeneous Dirichlet boundary condition: (2.1)

- ∇ ⋅ ( ρ ⁢ ∇ ⁡ u ) = f in ⁢ Ω ,

u = 0 on ⁢ ∂ ⁡ Ω ,

where ρ = ρ ⁢ ( x ) is a piecewise constant function satisfying 0 < ρ m ≤ ρ ⁢ ( x ) ≤ ρ M , f ∈ L 2 ⁢ ( Ω ) , and Ω is a simply connected and bounded polygonal domain.

Let T h be the shape regular triangulation of Ω with h = max T ∈ T h ⁡ h T , where h T is a diameter of an element 𝑇. In order to ease our analysis, we assume that T h is generated so that 𝜌 is a constant in each T ∈ T h . The set of edges in T h is denoted by E h and the set of interior edges by E h i . The skeleton of T h is denoted by K h which is the union of all edges, i.e., K h = ⋃ e ∈ E h e . In this paper, the abbreviation a ≲ b is used for the inequality a ≤ C ⁢ b with a generic positive constant which is independent of ℎ.

For a set D ⊂ R 2 , H s ⁢ ( D ) is the usual Sobolev space with s ≥ 0 , and ∥ ⋅ ∥ s , D , | ⋅ | s , D represent the norm and the semi-norm of H s ⁢ ( D ) , respectively. The broken Sobolev space is defined as H s ⁢ ( T h ) = ∏ T ∈ T h H s ⁢ ( T ) , with the norm ∥ ⋅ ∥ s , h := ( ∑ T ∈ T h ∥ ⋅ ∥ s , T 2 ) 1 / 2 . Abbreviations ( ⋅ , ⋅ ) D and ⟨ ⋅ , ⋅ ⟩ ∂ ⁡ D represent L 2 inner products on 𝐷 and ∂ ⁡ D , respectively.

Let P k ⁢ ( T ) and P k ⁢ ( e ) be the function spaces of standard Lagrange polynomial functions of order k ≥ 1 on 𝑇 and 𝑒, respectively. Then

V h k := { v : v | T ∈ P k ⁢ ( T ) , T ∈ T h } ,

and

Λ h k , D := { μ : μ | e ∈ P k ⁢ ( e ) , e ∈ E h } , Λ h k , C := { μ ∈ C ⁢ ( K h ) : μ ∈ Λ h k , D } .

We will use the notation Λ h k , ∙ if both Λ h k , C and Λ h k , D are available. The function spaces Λ h , 0 k , D and Λ h , 0 k , C represent subspaces of Λ h k , D and Λ h k , C with zero trace on ∂ ⁡ Ω , respectively. Following Brezzi and Fortin [6] or Roberts and Thomas [39], the Raviart–Thomas–Nédélec (RTN) finite-dimensional subspace of H ⁢ ( div , Ω ) is defined as

RTN k ⁢ ( T h ) := { v h ∈ H ⁢ ( div , Ω ) : v h | T ∈ RTN k ⁢ ( T ) ⁢ for all ⁢ T ∈ T h } ,

where RTN k ⁢ ( T ) := [ P k ⁢ ( T ) ] 2 + x ⁢ P k ⁢ ( T ) is the RTN space on an element 𝑇. Note that v ∈ RTN k ⁢ ( T ) satisfies ∇ ⋅ v h ∈ P k ⁢ ( T ) for all T ∈ T h and v h ⋅ n e ∈ P k ⁢ ( e ) for all e ∈ E h .

We now introduce a unified version of variational formulation: find u h ∈ V h k and λ h ∈ Λ h , 0 k , ∙ such that, for v h ∈ V h k and μ h ∈ Λ h , 0 k , ∙ ,

(2.2) B ⁢ ( u h , λ h ; v h , μ h ) = L ⁢ ( v h ) ,

where

B ⁢ ( u , λ ; v , μ ) = ∑ T ∈ T h { ( ρ ⁢ ∇ ⁡ u , ∇ ⁡ v ) T - ⟨ σ T ⁢ ( u , λ ) , v - μ ⟩ ∂ ⁡ T - S ⁢ ⟨ u - λ , ρ ⁢ ∂ n ⁡ v ⟩ ∂ ⁡ T } , L ⁢ ( v ) = ( f , v ) Ω .

Here, the numerical flux σ T ⁢ ( u , λ ) is defined as

σ T ⁢ ( u , λ ) = ρ ⁢ ∂ n ⁡ u | T - ρ ⁢ κ h T β ⁢ ( u - λ ) | T .

The parameter S determines the I-HDG ( S = 0 ), N-HDG ( S = - 1 ) and S-HDG ( S = 1 ) methods, and κ ≥ 0 , β ≥ 1 depend on the methods. For the S-HDG method, 𝛽 is chosen to be 1 in general. The penalty parameter 𝜅 has to be sufficiently large for the stability [45]. The penalty parameter depends on the polynomial degree and the shape of the corresponding element, but it is independent of ℎ. On the other hand, 𝜅 can be any non-negative constant for the N-HDG method. In this case, the L 2 convergence rate is influenced by the choice of 𝛽 (see [40]). For example, 𝛽 is chosen to be 3 to obtain the optimal L 2 convergence. For the I-HDG method, 𝜅 has to be sufficiently large similar to the S-HDG method, and 𝛽 facilitates the optimal convergence in the L 2 norm similar to the N-HDG method. These results can be obtained from [40] with a simple modification.

Each HDG method is divided into two types by choosing a trace variable. For the continuous trace, a vertex node is related to more than two vertex nodes of neighboring elements, while an edge node has relevance to two edge nodes of neighboring elements on the interior skeleton. However, for the discontinuous one, a vertex node and an edge node correlate to two edge and two vertex nodes of neighboring elements, respectively. Note that the number of degrees of freedom for the continuous trace is much less than that for the discontinuous trace because end points of an edge are shared for the continuous trace. If κ > 0 , both continuous and discontinuous traces can be applied to the HDG methods regardless of 𝑆. However, the continuous trace has to be chosen in the case S = - 1 and κ = 0 (the GCBE method [25]). Note that we do not consider the case S ≠ - 1 and κ = 0 because the I-HDG and S-HDG methods are not stable without the penalty term. Figure 1 shows degrees of freedom for the HDG methods with the continuous and discontinuous traces.

Figure 1

Degrees of freedom for quadratic HDG with continuous (left) and discontinuous traces (right).

Lemma 2.1

Lemma 2.1 (Weak Solution)

There exists a unique solution u ∈ H 0 1 ⁢ ( Ω ) satisfying

(2.3) ( ρ ⁢ ∇ ⁡ u , ∇ ⁡ v ) Ω = ( f , v ) Ω for all ⁢ v ∈ H 0 1 ⁢ ( Ω ) .

Proof

By the Lax–Milgram theorem, there exists a unique solution of (2.3). ∎

Let us denote the flux jump by ⟦ σ ( u h , λ h ) ⟧ e = σ T ( u h , λ h ) + σ T ′ ( u h , λ h ) on e = ∂ ⁡ T ∩ ∂ ⁡ T ′ . With this definition, we state the following local conservation properties.

Proposition 2.2

Let u h ∈ V h k , λ h ∈ Λ h k , D be the approximate solution and trace which satisfy (2.2). Then the local conservation property

(2.4) - ∫ ∂ ⁡ T σ T ⁢ ( u h , λ h ) ⁢ d s = - ∫ T ∇ ⋅ ( ρ ⁢ ∇ ⁡ u h ) ⁢ d x + ∫ ∂ ⁡ T ρ ⁢ κ h T β ⁢ ( u h - λ h ) ⁢ d s = ∫ T f ⁢ d x

and the average flux continuity

(2.5) ∫ e ⟦ σ ( u h , λ h ) ⟧ e d s = 0

hold for all T ∈ T h and e ∈ E h i .

Proof

It is easy to see that (2.4) holds by taking v = 1 in 𝑇, v = 0 in T h ∖ T and μ = 0 in (2.2). The average flux continuity (2.5) can be obtained with v = 0 , μ = 1 on 𝑒 and μ = 0 on K h ∖ e in (2.2). ∎

In the case that λ h ∈ Λ h k , C , the local conservation property (2.4) holds, but the average flux continuity (2.5) does not hold. In this case, the total numerical flux is continuous, i.e.,

(2.6) ∑ e ∈ E h i ∫ e ⟦ σ ( u h , λ h ) ⟧ e μ h d s = ∑ T ∈ T h ∫ ∂ ⁡ T σ T ( u h , λ h ) μ h d s = 0 ,

by taking v = 0 and μ h ∈ Λ h , 0 k , C in (2.2).

3 A Posteriori Error Estimate

Prior to a posteriori error estimates, we introduce an a priori error estimate for the HDG methods. Let | | | v | | | be the energy norm defined by

(3.1) | | | v | | | := ( ∑ T ∈ T h | | | v | | | T 2 ) 1 / 2 , | | | v | | | T 2 := ∥ ρ 1 / 2 ⁢ ∇ ⁡ v ∥ 0 , T 2 + κ h T β ⁢ ∥ ρ 1 / 2 ⁢ ( v - μ ) ∥ 0 , ∂ ⁡ T 2 ,

where 𝜇 is the trace variable with respect to 𝑣. Note that the energy norm is represented as

| | | v | | | := ( ∑ T ∈ T h ∥ ρ 1 / 2 ⁢ ∇ ⁡ v ∥ 0 , T 2 ) 1 / 2

for the GCBE method. With these definitions, we have

and

(3.3) | | | u - u h | | | ≲ h k ⁢ ∥ u ∥ k + 1 , Ω ,

where u ∈ H 0 1 ⁢ ( Ω ) ∩ H k + 1 ⁢ ( Ω ) and u h are the weak and approximate solutions to (2.2), respectively. Proofs of (3.2) and (3.3) can be found in [25, 45]. Note that the super-penalty parameter 𝛽 does not affect the coercivity (3.2) and the error estimate in the energy norm (3.3). Since we focus on the energy error to analyze the a posteriori error estimate, we fix β = 1 in the remainder of the paper.

The HDG solution of (2.2) does not belong to H 0 1 ⁢ ( Ω ) . Thus, we need its correction, called the potential reconstruction, in order to analyze the a posteriori error estimate.

Definition 3.1

Definition 3.1 (Potential Reconstruction)

Let u h be the approximate solution of (2.2). Then the function s h is called the potential reconstruction constructed from u h , which satisfies s h ∈ H 0 1 ⁢ ( Ω ) .

In order to construct a potential reconstruction, we define I av by the average operator such that

I av ⁢ ( u h ) ⁢ ( a ) = 1 | T a | ⁢ ∑ T ∈ T a u h | T ⁢ ( a ) ,

where T a denotes all the elements sharing a node 𝒂 and | T a | represents their number. Note that we have I av ⁢ ( u h ) ⁢ ( a ) = u h ⁢ ( a ) when the node 𝒂 is in the interior of an element, and I av ⁢ ( u h ) ⁢ ( a ) = 0 when the node 𝒂 is on the boundary ∂ ⁡ Ω . Then the potential reconstruction s h is constructed as follows:

(3.4) s h := I av ⁢ ( u h ) .

Note that s h ∈ H 0 1 ⁢ ( Ω ) by the definition of average operator.

3.1 A Residual Type A Posteriori Error Estimate

In this section, a residual type error estimator is constructed and analyzed. In order to ease analysis, we introduce projection operators

π k , C : H 0 1 ⁢ ( Ω ) → Λ h , 0 k , C and π k , D : H 0 1 ⁢ ( Ω ) → Λ h , 0 k , D

as follows. Let π k , D be the edgewise L 2 orthogonal projection onto Λ h , 0 k , D , and let I C k : Λ h , 0 k , D → Λ h , 0 k , C be the averaging operator (cf. (3.4)). Then, for k ≥ 3 , we define π k , C ⁢ ϕ | e := I C k ⁢ ( π k , D ⁢ ϕ ) | e + c ⁢ b ⁢ ( x ) , where b ∈ H 0 1 ⁢ ( e ) is a polynomial bubble function on 𝑒. Here, the constant 𝑐 is determined for a given piecewise polynomial ⟦ σ T ( u h , λ h ) ⟧ ∈ P k ( e ) to satisfy

( π k , C ϕ , ⟦ σ T ( u h , λ h ) ⟧ ) e = ( ϕ , ⟦ σ T ( u h , λ h ) ⟧ ) e

on each e ∈ K h ∖ ∂ ⁡ Ω . By construction, π k , C ⁢ ϕ is well-defined for k ≥ 3 . Some modifications are required for k = 2 . Indeed, note that ( b ⁢ ( x ) , w h ) e = 0 if w h is linear passing through the midpoint on 𝑒 due to the symmetricity of a quadratic bubble function. In this case, we can set

I C k ⁢ ( π k , D ⁢ v ) | e = π k , D ⁢ v on ⁢ e ,

I C k ⁢ ( π k , D ⁢ v ) | e ′ = I C k ⁢ ( π k , D ⁢ v ) | e at the intersection point of ⁢ e ⁢ and ⁢ e ′ .

In practice, ( b ⁢ ( x ) , w h ) e ≠ 0 except for a pathological case where ⟦ σ T ( u h , λ h ) ⟧ = a ( x - x M ) for some a ≠ 0 . Here, x M is the midpoint on the edge.

Theorem 3.2

Theorem 3.2 (A Residual Type A Posteriori Error Estimate)

Let 𝑢 be the solution of (2.1), and let

( u h , λ h ) ∈ V h k × Λ h , 0 k , ∙

be the approximate solution pair of (2.2). Let s h be the potential reconstruction defined by (3.4). Then the following holds for all T ∈ T h :

(3.5) | | | u - u h | | | 2 ≤ ∑ T ∈ T h { ( η RE , T + C t , ∂ ⁡ T κ 1 / 2 η tr , T ) 2 + η NC , T 2 + η tr , T 2 } = : η T 2 ,

with the element residual estimator (3.6)

(3.6a) η RE , T := C P , T ⁢ h T ⁢ ρ - 1 / 2 ⁢ ∥ f + ∇ ⋅ ( ρ ⁢ ∇ ⁡ u h ) ∥ 0 , T ,

the nonconformity estimator

(3.6b) η NC , T := ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u h - s h ) ∥ 0 , T

and the trace estimator

(3.6c) η tr , T := ρ 1 / 2 ⁢ κ 1 / 2 h T 1 / 2 ⁢ ∥ u h - λ h ∥ 0 , ∂ ⁡ T .

Proof

By the Riesz representation theorem, there exists the unique function s ∈ H 0 1 ⁢ ( Ω ) such that

(3.7) ( ρ ⁢ ∇ ⁡ s , ∇ ⁡ v ) Ω = ( ρ ⁢ ∇ h ⁡ u h , ∇ ⁡ v ) Ω for all ⁢ v ∈ H 0 1 ⁢ ( Ω ) ,

where ∇ h denotes the broken gradient operator with respect to T h . Using (3.7), we have

∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - u h ) ∥ 0 , Ω 2 = ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - s + s - u h ) ∥ 0 , Ω 2 = ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - s ) ∥ 0 , Ω 2 + ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( s - u h ) ∥ 0 , Ω 2 + 2 ⁢ ( ρ ⁢ ∇ ⁡ ( s - u h ) , ∇ ⁡ ( u - s ) ) Ω .

Since 𝑢 and 𝑠 belong to the space H 0 1 ⁢ ( Ω ) , we have 2 ⁢ ( ρ ⁢ ∇ ⁡ ( s - u h ) , ∇ ⁡ ( u - s ) ) Ω = 0 . Hence

(3.8) ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - u h ) ∥ 0 , Ω 2 = ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - s ) ∥ 0 , Ω 2 + ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( s - u h ) ∥ 0 , Ω 2 .

Step 1: ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( s - u h ) ∥ 0 , Ω 2 ≤ ∑ T ∈ T h η NC , T 2 . For any function v ∈ H 0 1 ⁢ ( Ω ) , we have property (3.8), i.e.,

∥ ρ 1 / 2 ⁢ ∇ ⁡ ( v - u h ) ∥ 0 , Ω 2 = ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( v - s ) ∥ 0 , Ω 2 + ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( s - u h ) ∥ 0 , Ω 2 .

Then

∥ ρ 1 / 2 ⁢ ∇ ⁡ ( s - u h ) ∥ 0 , Ω 2 = ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( v - u h ) ∥ 0 , Ω 2 - ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( v - s ) ∥ 0 , Ω 2 ≤ ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( v - u h ) ∥ 0 , Ω 2

for all v ∈ H 0 1 ⁢ ( Ω ) . Since s h ∈ H 0 1 ⁢ ( Ω ) , we have

∥ ρ 1 / 2 ⁢ ∇ ⁡ ( s - u h ) ∥ 0 , Ω 2 ≤ ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( s h - u h ) ∥ 0 , Ω 2 = ∑ T ∈ T h ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( s h - u h ) ∥ 0 , T 2 .

Step 2: ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - s ) ∥ 0 , Ω 2 ≤ ∑ T ∈ T h ( η RE , T + C t , ∂ ⁡ T ⁢ κ 1 / 2 ⁢ η tr , T ) 2 . With the definition of the dual norm, we have

(3.9) ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - s ) ∥ = sup ϕ ∈ H 0 1 ⁢ ( Ω ) ∥ ρ 1 / 2 ⁢ ∇ ⁡ ϕ ∥ = 1 ⁡ ( ρ ⁢ ∇ ⁡ ( u - s ) , ∇ ⁡ ϕ ) = sup ϕ ∈ H 0 1 ⁢ ( Ω ) ∥ ρ 1 / 2 ⁢ ∇ ⁡ ϕ ∥ = 1 ⁡ ( ρ ⁢ ∇ ⁡ ( u - u h ) , ∇ ⁡ ϕ ) .

Using integration by parts, local conservation property (2.4), flux continuity (2.6), the projection operators π k , ∙ , the Poincaré inequality and the trace inequality, we have

( ρ ⁢ ∇ ⁡ ( u - u h ) , ∇ ⁡ ϕ ) Ω = ∑ T ∈ T h { ( f + ∇ ⋅ ( ρ ⁢ ∇ ⁡ u h ) , ϕ ) T - ⟨ σ T ⁢ ( u h , λ h ) , ϕ ⟩ ∂ ⁡ T - ρ ⁢ κ h T ⁢ ⟨ u h - λ h , ϕ ⟩ ∂ ⁡ T } = ∑ T ∈ T h { ( f + ∇ ⋅ ( ρ ⁢ ∇ ⁡ u h ) , ϕ - ϕ T ) T - ⟨ σ T ⁢ ( u h , λ h ) , π k , ∙ ⁢ ϕ ⟩ ∂ ⁡ T - ρ ⁢ κ h T ⁢ ⟨ u h - λ h , ϕ - ϕ T ⟩ ∂ ⁡ T } ≤ ∑ T ∈ T h { ∥ f + ∇ ⋅ ( ρ ⁢ ∇ ⁡ u h ) ∥ 0 , T ⁢ ∥ ϕ - ϕ T ∥ 0 , T + ρ ⁢ κ h T ⁢ ∥ u h - λ h ∥ 0 , ∂ ⁡ T ⁢ ∥ ϕ - ϕ T ∥ 0 , ∂ ⁡ T } ≤ ∑ T ∈ T h { C P , T ⁢ h T ⁢ ρ - 1 / 2 ⁢ ∥ f + ∇ ⋅ ( ρ ⁢ ∇ ⁡ u h ) ∥ 0 , T + C t , ∂ ⁡ T ⁢ κ 1 / 2 ⁢ ρ 1 / 2 ⁢ κ 1 / 2 h T 1 / 2 ⁢ ∥ u h - λ h ∥ 0 , ∂ ⁡ T } ⁢ ∥ ρ 1 / 2 ⁢ ∇ ⁡ ϕ ∥ 0 , T ,

where ϕ T is the mean value of 𝜙 over 𝑇 given by ϕ T := 1 | T | ⁢ ∫ T ϕ ⁢ d x . Here, C P , T can be evaluated as 1 π in view of the convexity of 𝑇 (cf. [4, 36]) and C t , ∂ ⁡ T 2 = 6 ⁢ h T ⁢ | ∂ ⁡ T | | T | , where | ∂ ⁡ T | is the perimeter length of 𝑇 and | T | is the area of 𝑇 (see [21]). Using the Cauchy–Schwarz inequality and the fact that ∥ ρ 1 / 2 ⁢ ∇ ⁡ ϕ ∥ = 1 yields

∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - s ) ∥ 2 ≤ ∑ T ∈ T h ( η RE , T + C t , ∂ ⁡ T ⁢ κ 1 / 2 ⁢ η tr , T ) 2 .

Therefore, the proof is complete by using the definition of the energy norm (3.1) and summing the results of step 1 and step 2. ∎

Estimator (3.5) does not cover the linear approximation of the HDG methods with the continuous trace because the projection operator π 1 , C does not exist. On the other hand, the estimator covers the GCBE method and all the HDG methods with the discontinuous trace. Especially, for the GCBE method ( S = - 1 , κ = 0 , k ≥ 2 ), the error estimator has a simple form

(3.10) | | | u - u h | | | 2 ≤ ∑ T ∈ T h { η RE , T 2 + η NC , T 2 } .

For all T ∈ T h , let us define E T by the set of all the edges in E h sharing at least a vertex with 𝑇. Then the jump estimator can be defined as

| u h | J , T := { ∑ e ∈ E T h e - 1 ⁢ ∥ [ u h ] ∥ e 2 } 1 / 2 ,

where [ ϕ ] := ϕ | T ℓ - ϕ | T r denotes the jump of 𝜙 across the interface boundary e = ∂ ⁡ T ℓ ∩ ∂ ⁡ T r . On each 𝑒,

∥ [ u h ] ∥ e 2 = ∥ u h | T ℓ - u h | T r ∥ e 2 ≤ ∥ u h | T ℓ - λ h ∥ e 2 + ∥ u h | T r - λ h ∥ e 2 .

Thus, we have

(3.11) | u h | J , T 2 ≤ ∑ e ∈ E T h e - 1 ⁢ ( ∥ u h | T ℓ - λ h ∥ e 2 + ∥ u h | T r - λ h ∥ e 2 ) ≲ ∑ T ^ ∈ T T κ h T ^ ⁢ ∥ ρ 1 / 2 ⁢ ( u h - λ h ) ∥ 0 , ∂ ⁡ T ^ 2 ,

where T T is the set of all elements in T h sharing at least an edge e ∈ ∂ ⁡ T .

We now introduce the approximation property of the averaging operator. This result can be proved as in [28, 7].

Lemma 3.3

Lemma 3.3 (Averaging Operator)

Let u h ∈ V h k , k ≥ 1 . Then

(3.12) ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u h - I av ⁢ ( u h ) ) ∥ 0 , T ≲ ρ 1 / 2 ⁢ | u h | J , T

for all T ∈ T h .

Also, we have efficiency of the element residual following [43] with a simple modification. Here, we assume that 𝑓 is a piecewise polynomial function in V h k .

Lemma 3.4

Lemma 3.4 (Efficiency of the Element Residuals)

Let T ∈ T h . Then

(3.13) h T ⁢ ∥ f + ∇ ⋅ ( ρ ⁢ ∇ ⁡ u h ) ∥ 0 , T ≲ ρ 1 / 2 ⁢ ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - u h ) ∥ 0 , T .

Then we have the following theorem for efficiency.

Theorem 3.5

Theorem 3.5 (Local Efficiency of Estimate (3.5))

Let 𝑢 be the weak solution of (2.3) and ( u h , λ h ) ∈ V h k × Λ h , 0 k , ∙ the approximate solution pair of (2.2). Let η T be the local estimator given by (3.5). Then

η T 2 ≲ ∑ T ∈ T T | | | u - u h | | | T 2 .

Proof

Lemma 3.4 and the definitions of η RE , T and η tr , T yield

( η RE , T + C t , ∂ ⁡ T ⁢ κ 1 / 2 ⁢ η tr , T ) 2 ≲ ( ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - u h ) ∥ 0 , T + κ 1 / 2 ⁢ η tr , T ) 2 ≲ ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - u h ) ∥ 0 , T 2 + η tr , T 2 .

Since the potential reconstruction is constructed to be the same as the average operator (3.4), (3.12) can be directly applied,

(3.14) η NC , T 2 ≲ ρ ⁢ | u h | J , T 2 .

Combining the above estimates implies

( η RE , T + C t , ∂ ⁡ T ⁢ κ 1 / 2 ⁢ η tr , T ) 2 + η NC , T 2 + η tr , T 2 ≲ ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - u h ) ∥ 0 , T 2 + η tr , T 2 + ρ ⁢ | u h | J , T 2 .

Note that | u - u h | J , T 2 = | u h | J , T since [ u ] = 0 for all e ∈ E h . By combining the above estimate and (3.11), we have (3.4). ∎

For the GCBE method, we have the following local efficiency whose proof is straightforward by (3.12), (3.13).

Corollary 3.6

Corollary 3.6 (Local Efficiency of Estimate (3.10))

Let 𝑢 be the weak solution of (2.3) and ( u h , λ h ) ∈ V h k × Λ h , 0 k , C the approximate solution pair for the GCBE method. Let η RE , T and η NC , T be the local estimators given by (3.6). Then

(3.15) η RE , T 2 + η NC , T 2 ≲ ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - u h ) ∥ 0 , T 2 + ρ ⁢ | u h | J , T 2 .

Here, (3.15) does not satisfy general local efficiency result. Instead, we have the following equivalence results satisfying both reliability and local efficiency (see [20]):

∑ T ∈ T h ( | | | u - u h | | | T 2 + ρ ⁢ | u - u h | J , T 2 ) ≤ ∑ T ∈ T h ( η T 2 + ρ ⁢ | u h | J , T 2 ) , η T 2 + ρ ⁢ | u h | J , T 2 ≲ | | | u - u h | | | T 2 + ρ ⁢ | u - u h | J , T 2 .

3.2 Flux Reconstruction Based on the A Posteriori Error Estimate

In the previous section, the the effectivity index or estimate (3.5) depends on the penalty parameter 𝜅. Thus, a residual type error estimator is not a robust estimator. In this section, we suggest a robust estimator for the HDG methods with the discontinuous trace by using a flux reconstruction.

Definition 3.7

Definition 3.7 (Flux Reconstruction)

Let ( u h , λ h ) ∈ V h k × Λ h , 0 k , D be the approximate solutions of (2.2). Then the function σ h is called the equilibrated flux reconstruction constructed from u h and λ h , which satisfies (3.16)

(3.16a) σ h ∈ H ⁢ ( div , Ω ) ,

(3.16b) ( - ∇ ⋅ σ h , v h ) T = ( f , v h ) T for all ⁢ v h ∈ V h k - 1 ⁢ and all ⁢ T ∈ T h .

From (3.16b), it is clear that ( - ∇ ⋅ σ h , 1 ) T = ( f , 1 ) T by taking v h = 1 in 𝑇. Using integration by parts and (2.2), we have

( σ h , ∇ ⁡ v h ) T - ⟨ σ h ⋅ n , v h ⟩ ∂ ⁡ T = ( ρ ⁢ ∇ ⁡ u h , ∇ ⁡ v h ) T - ⟨ σ T ⁢ ( u h , λ h ) , v h ⟩ ∂ ⁡ T - S ⁢ ⟨ u h - λ h , ρ ⁢ ∂ n ⁡ v h ⟩ ∂ ⁡ T

for all v h ∈ V h k - 1 . Thus, we can construct the equilibrated flux reconstruction by using the above equation, (2.5) and (3.16a), (3.17)

(3.17a) ⟨ σ h ⋅ n e , v h ⟩ e = ⟨ ρ ⁢ ∂ n ⁡ u h , v h ⟩ e - κ h T ⁢ ⟨ ρ ⁢ ( u h - λ h ) , v h ⟩ e ,

(3.17b) ( σ h , ∇ ⁡ v h ) T = ( ρ ⁢ ∇ ⁡ u h , ∇ ⁡ v h ) T - S ⁢ ⟨ u h - λ h , ρ ⁢ ∂ n ⁡ v h ⟩ ∂ ⁡ T

for all v h ∈ V h k - 1 , e ∈ E h and T ∈ T h .

Proposition 3.8

Proposition 3.8 (Flux Reconstruction for the HDG Method)

Let ( u h , λ h ) ∈ V h k × Λ h , 0 k , D be the approximate solution pair of (2.2). The equilibrated flux reconstruction σ h ∈ RTN k - 1 ⁢ ( T h ) is defined for all T ∈ T h , e ∈ E T and q h ∈ P k - 1 ⁢ ( e ) , (3.18)

(3.18a) ⟨ σ h ⋅ n e , q h ⟩ e = ⟨ ρ ⁢ ∂ n ⁡ u h , q h ⟩ e - κ h T ⁢ ⟨ ρ ⁢ ( u h - λ h ) , q h ⟩ e ,

and for all r h ∈ [ P k - 2 ⁢ ( T ) ] 2 ,

(3.18b) ( σ h , r h ) T = ( ρ ⁢ ∇ ⁡ u h , r h ) T - S ⁢ ⟨ u h - λ h , ρ ⁢ r h ⋅ n T ⟩ ∂ ⁡ T .

Then σ h satisfies (3.16).

Clearly, σ h satisfies (3.16a) with the definition of RTN k - 1 ⁢ ( T h ) . The remaining property of σ h (3.16b) can be proved by choosing q h = v h and r h = ∇ ⁡ v h in (3.17). From (3.4) and Proposition 3.8, we have the potential reconstruction s h and the equilibrated flux reconstruction σ h . Then we can use the following a posteriori error estimate.

Theorem 3.9

Theorem 3.9 (A Posteriori Error Estimate Based on the Flux Reconstruction)

Let 𝑢 be the weak solution of (2.3) and ( u h , λ h ) ∈ V h k × Λ h , 0 k , D the approximate solution pair of (2.2). Let s h be the potential reconstruction defined by (3.4) and σ h the equilibrated flux reconstruction constructed by Proposition 3.8. Then the following holds for all T ∈ T h :

| | | u - u h | | | 2 ≤ ∑ T ∈ T h { ( η R , T + η F , T ) 2 + η NC , T 2 + η tr , T 2 } = : η T 2 ,

with the residual estimator (3.19)

(3.19a) η R , T := C P , T ⁢ h T ⁢ ρ - 1 / 2 ⁢ ∥ f + ∇ ⋅ σ h ∥ 0 , T ,

the flux estimator

(3.19b) η F , T := ∥ ρ 1 / 2 ⁢ ∇ ⁡ u h - ρ - 1 / 2 ⁢ σ h ∥ 0 , T ,

and η NC , T and η tr , T are the same as (3.6b) and (3.6c), respectively.

Proof

The proof of the theorem basically follows the argument of the proof of Theorem 3.2. Step 1 is the same as before. The remaining part is to show the following.

Step 2: ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - s ) ∥ 0 , Ω 2 ≤ ∑ T ∈ T h ( η R , T + η F , T ) 2 . Using (2.3), (3.9) and (3.16), we have

( ρ ⁢ ∇ ⁡ ( u - u h ) , ∇ ⁡ ϕ ) Ω = ∑ T ∈ T h ( - ∇ ⋅ ( ρ ⁢ ∇ ⁡ u - σ h ) , ϕ ) T + ∑ T ∈ T h ( ρ - 1 / 2 ⁢ σ h - ρ 1 / 2 ⁢ ∇ ⁡ u h , ρ 1 / 2 ⁢ ∇ ⁡ ϕ ) T ≤ ∑ T ∈ T h ∥ f + ∇ ⋅ σ h ∥ 0 , T ⁢ ∥ ϕ - ϕ T ∥ 0 , T + ∑ T ∈ T h ∥ ρ 1 / 2 ⁢ ∇ ⁡ u h - ρ - 1 / 2 ⁢ σ h ∥ 0 , T ⁢ ∥ ρ 1 / 2 ⁢ ∇ ⁡ ϕ ∥ 0 , T ≤ ∑ T ∈ T h ( C P , T ⁢ h T ⁢ ρ - 1 / 2 ⁢ ∥ f + ∇ ⋅ σ h ∥ 0 , T + ∥ ρ 1 / 2 ⁢ ∇ ⁡ u h - ρ - 1 / 2 ⁢ σ h ∥ 0 , T ) ⁢ ∥ ρ 1 / 2 ⁢ ∇ ⁡ ϕ ∥ 0 , T .

Using the Cauchy–Schwarz inequality and the fact that ∥ ρ 1 / 2 ⁢ ∇ ⁡ ϕ ∥ = 1 yields

∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - s ) ∥ 2 ≤ sup ϕ ∈ H 0 1 ⁢ ( Ω ) ∥ ρ 1 / 2 ⁢ ∇ ⁡ ϕ ∥ = 1 ⁡ ∑ T ∈ T h ( η R , T + η F , T ) 2 ⁢ ∑ T ∈ T h ∥ ρ 1 / 2 ⁢ ∇ ⁡ ϕ ∥ 0 , T 2 = ∑ T ∈ T h ( η R , T + η F , T ) 2 .

Since u - λ = 0 on ∂ ⁡ T , we have

| | | u - u h | | | 2 = ∑ T ∈ T h ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - u h ) ∥ 0 , T 2 + ∑ T ∈ T h κ h T ⁢ ∥ ρ 1 / 2 ⁢ ( u h - λ h ) ∥ 0 , ∂ ⁡ T 2 ≤ ∑ T ∈ T h ( η R , T + η F , T ) 2 + ∑ T ∈ T h η NC , T 2 + ∑ T ∈ T h η tr , T 2 . ∎

We now define a projection Π h ℓ : RTN ℓ + 1 ⁢ ( T ) → [ P ℓ ⁢ ( T ) ] 2 such that

(3.20) ( Π h ℓ ⁢ v h , r h ) T = ( v h , r h ) T for all ⁢ r h ∈ [ P ℓ ⁢ ( T ) ] 2 ,

for v h ∈ RTN ℓ + 1 ⁢ ( T ) . Clearly, (3.20) is uniquely solvable. Using (3.20), scaling arguments and the Piola transformation, we have

(3.21) ∥ v h ∥ 0 , T 2 ≲ h T ⁢ ∑ e ∈ E T ∥ v h ⋅ n e ∥ 0 , e 2 + ∥ Π h ℓ ⁢ v h ∥ 0 , T 2

for all T ∈ T h and v h ∈ RTN ℓ + 1 ⁢ ( T ) . This result can be found in [19].

Theorem 3.10

Theorem 3.10 (Efficiency of the Estimate of Theorem 3.9)

Let 𝑢 be the weak solution of (2.3), and let

( u h , λ h ) ∈ V h k × Λ h , 0 k , D

be the approximate solution of (2.2). Let η R , T , η F , T , η NC , T and η tr , T be the local estimators given by (3.19). Then

( η R , T + η F , T ) 2 + η NC , T 2 + η tr , T 2 ≲ ∑ T ∈ T T | | | u - u h | | | T 2 .

Proof

Let T ∈ T h be fixed. Then the triangle inequality, the inverse estimate and (3.13) yield that

η R , T = C P , T ⁢ ρ - 1 / 2 ⁢ h T ⁢ ∥ f + ∇ ⋅ ( ρ ⁢ ∇ ⁡ u h ) - ∇ ⋅ ( ρ ⁢ ∇ ⁡ u h ) + ∇ ⋅ σ h ∥ 0 , T , ≤ C P , T ⁢ ρ - 1 / 2 ⁢ h T ⁢ ∥ f + ∇ ⋅ ( ρ ⁢ ∇ ⁡ u h ) ∥ 0 , T + C P , T ⁢ ρ - 1 / 2 ⁢ h T ⁢ ∥ ∇ ⋅ ( ρ ⁢ ∇ ⁡ u h - σ h ) ∥ 0 , T ≲ C P , T ⁢ ρ - 1 / 2 ⁢ h T ⁢ ∥ f + ∇ ⋅ ( ρ ⁢ ∇ ⁡ u h ) ∥ 0 , T + ∥ ρ 1 / 2 ⁢ ∇ ⁡ u h - ρ - 1 / 2 ⁢ σ h ∥ 0 , T ≲ ∥ ρ 1 / 2 ⁢ ∇ ⁡ ( u - u h ) ∥ T + η F , T .

Let us denote v h = ( ρ ⁢ ∇ ⁡ u h - σ h ) | T . Clearly, v h ∈ RTN k - 1 ⁢ ( T ) and v h ⋅ n e ∈ P k - 1 ⁢ ( e ) . By taking v h ⋅ n e in (3.18a) and using the Cauchy–Schwarz inequality,

(3.22) ∥ v h ⋅ n e ∥ 0 , e 2 = κ h T ⁢ ⟨ ρ ⁢ ( u h - λ h ) , v h ⋅ n e ⟩ e ≤ ρ 1 / 2 ⁢ κ h T ⁢ ∥ ρ 1 / 2 ⁢ ( u h - λ h ) ∥ 0 , e ⁢ ∥ v h ⋅ n e ∥ 0 , e .

By (3.18b) and (3.20),

( Π h k - 2 ⁢ v h , r h ) T = ( v h , r h ) T = S ⁢ ⟨ u h - λ h , ρ ⁢ r h ⋅ n T ⟩ ∂ ⁡ T .

Then we have the following estimate with r h = Π h k - 2 ⁢ v h and the trace inequality:

(3.23) ∥ Π h k - 2 ⁢ v h ∥ 0 , T 2 ≤ ρ 1 / 2 ⁢ ∥ ρ 1 / 2 ⁢ ( u h - λ h ) ∥ 0 , ∂ ⁡ T ⁢ ∥ Π h k - 2 ⁢ v h ⋅ n T ∥ 0 , ∂ ⁡ T ≲ max ⁡ ( 1 , κ - 1 / 2 ) ⁢ ρ 1 / 2 ⁢ κ 1 / 2 h T 1 / 2 ⁢ ∥ ρ 1 / 2 ⁢ ( u h - λ h ) ∥ 0 , ∂ ⁡ T ⁢ ∥ Π h k - 2 ⁢ v h ∥ 0 , T .

Combining (3.21), (3.22) and (3.23), we have

η F , T 2 ≲ ( max ⁡ ( 1 , κ - 1 ) + κ ) ⁢ κ h T ⁢ ∥ ρ 1 / 2 ⁢ ( u h - λ h ) ∥ 0 , ∂ ⁡ T 2 .

Combining the above estimates, (3.11), (3.14) and the trace estimator, we complete the proof. ∎

4 Numerical Experiments

In this section, the performance of error estimators is tested with several numerical experiments. We compare their results between low- and high-order approximations, and between continuous and discontinuous traces with the estimators given in Theorem 3.2 and 3.9. We present convergence plots in terms of the number of degrees of freedom. In this case, the expected optimal rate of convergence is k / 2 for order 𝑘 polynomial approximation because the number of degrees of freedom is proportional to h - 2 in the uniform triangulation. The stabilization parameter 𝜅 is chosen appropriately to make the HDG methods stable.

The approximate solution u h and trace λ h are obtained by solving (2.2) on adaptive meshes which are obtained as follows.

Algorithm

Algorithm (Adaptive Algorithm)

Set the initial triangulation T 0 .
Solve (2.2), and compute error estimators for given T ℓ .
Compute a subset M ℓ ⊆ T ℓ for a given marking strategy.
Design the triangulation T ℓ + 1 to refine the elements in M ℓ . Hanging nodes are eliminated with additional refinements. Go to (2).

Before we state commonly used marking strategies such as average, maximum and bulk criteria, we introduce a new marking strategy based on unsupervised learning. There are many methods introduced in papers such as [23, 5] for the unsupervised classification problem. The 𝐾-means algorithm is the most widely used because of its simple implementation (cf. [22, 35]). Usually, one of the demerits of the 𝐾-means clustering is the choice of the number of clusters, 𝐾. However, we have only two clusters ( K = 2 ) for the marking strategy. Since the 𝐾-means algorithm gives the centroid of each partition, the cluster whose centroid is larger will be marked. Our 𝐾-means algorithm is the following.

Algorithm

Algorithm (𝐾-Means Clustering)

Compute max T ⁡ η T , and set two initial seeds, 0 and max T ⁡ η T .
Create two clusters by pairing every training data to the nearest seed.
Compute new seeds as the centroid of each cluster.
Repeat (2) and (3) until the tolerance criterion is satisfied.

Now we state our 𝐾-means criterion along with the standard marking strategies.

𝐾-means criterion: all elements in the cluster whose centroid is larger are marked.
Average criterion: all elements for which
η T > ( ∑ T ∈ T h η T ) / | T h |
are marked, where | T h | is the number of elements in T h .
Maximum criterion: for given θ M ∈ ( 0 , 1 ) , all elements for which
η T > ( 1 - θ M ) ⁢ max T ∈ T h ⁡ η T
are marked.
Bulk criterion: for given θ M ∈ ( 0 , 1 ) and sorted elements so that η T 1 ≥ η T 2 ≥ ⋯ ≥ η T | T h | , the first 𝑘 elements for which
∑ i = 1 k η T i 2 ≥ θ M ⁢ ∑ i = 1 | T h | η T i 2 > ∑ i = 1 k - 1 η T i 2
are marked.

For the maximum and bulk criteria, the parameter θ M determines how many elements are refined at each level. If θ M ≈ 0 , few elements are marked. On the other hand, if θ M ≈ 1 , the triangulation T h will be refined globally.

We start with the Poisson equation defined in the L-shaped domain.

Example 4.1

Consider the domain Ω = ( - 1 , 1 ) 2 ∖ ( [ 0 , 1 ] × [ - 1 , 0 ] ) and the diffusivity ρ = 1 . The source term is f = 0 with the exact solution

(4.1) u ⁢ ( r , θ ) = r 2 / 3 ⁢ sin ⁡ ( 2 3 ⁢ θ ) .

$Figure 2 Adaptive meshes with θ M = 0.5 \theta_{M}=0.5 : residual type estimator of GCBE method for k = 2 k=2 (top-left) and k = 6 k=6 (top-right), and flux reconstruction based estimator of N-HDG method for k = 2 k=2 (bottom-left) and k = 6 k=6 (bottom-right).$

Figure 2

Adaptive meshes with θ M = 0.5 : residual type estimator of GCBE method for k = 2 (top-left) and k = 6 (top-right), and flux reconstruction based estimator of N-HDG method for k = 2 (bottom-left) and k = 6 (bottom-right).

The exact solution of Example 4.1 has a corner singularity at the origin. It is well known that the optimal rate of convergence is 1 3 with respect to the number of degrees of freedom by using uniform refinement. Since the solution is smooth apart from the origin, we can expect that the error is large near the origin. Thus, those elements near the origin are refined intensively by an adaptive algorithm as shown in Figure 2. This figure also shows that low-order approximations tend to be overrefined in comparison to high-order approximations. For each polynomial order, adaptive meshes obtained from the residual type and the flux reconstruction based error estimators are similar to each other. Note that adaptive meshes of other methods are similar to those in Figure 2, but we do not present them here.

Figure 3

Effectivity indices with respect to 𝜅 for the HDG methods with the residual type estimator.

Figure 4

Effectivity indices with respect to 𝜅 for the HDG methods with the flux reconstruction based estimator.

Numerical experiments were performed with various penalty parameters ( κ = 1 , 10, 50, 100, 200, 400, 600, 800, 1000) in order to check the robustness for the corresponding estimator. The results of these experiments are displayed by using the linear interpolation in Figures 3 and 4. Figure 3 shows that effectivity indices for the residual type estimator depend on the penalty parameter 𝜅. They are much larger than 1 if 𝜅 is large. On the other hand, effectivity indices for the flux reconstruction based estimator are close to 1 regardless of 𝜅. The effectivity index slightly increases when 𝜅 becomes larger, but the difference is not significant. In Figure 4, the result for the N-HDG_C method is displayed to compare the residual type with the flux reconstruction based estimators. From these numerical experiments, we can say the robustness with respect to 𝜅 is satisfied for the flux reconstruction based estimator, but it is not satisfied for the residual type estimator. However, the residual type estimator for the GCBE method is robust because there is no penalty term.

$Figure 5 Convergence history (left) and effectivity index (right) for GCBE method with θ M = 0.5 \theta_{M}=0.5 .$

Figure 5

Convergence history (left) and effectivity index (right) for GCBE method with θ M = 0.5 .

$Figure 6 Convergence history (left) and effectivity index (right) for N-HDG method with θ M = 0.5 \theta_{M}=0.5 .$

Figure 6

Convergence history (left) and effectivity index (right) for N-HDG method with θ M = 0.5 .

Figure 5 shows the convergence of the error | | | u - u h | | | + | u - u h | J and the effectivity index for the GCBE method. Here, each marker represents the error in level 5 ⁢ ℓ , ℓ = 0 , 1 , … . We can check that the optimal rate of convergence is recovered by the adaptive algorithm in this case. Although the residual type estimator helps recover the optimal convergence rate and design adaptive meshes, effectivity indices are rather large for high-order approximations. In the undisplayed results obtained from other HDG methods with continuous traces, the optimal convergence rates are recovered. However, their effectivity indices depend on the penalty parameter 𝜅 as shown in Figure 3.

On the other hand, the flux reconstruction based estimator gives better performances than the residual type estimator. Similar to the GCBE method, the optimal rates of convergence are observed in Figure 6 for the N-HDG method with κ = 1 . Also, similar behavior is observed for the I-HDG and S-HDG methods, but we do not display them here. For the flux reconstruction based estimator, effectivity indices are especially worth noting. These values are close to 1 when the number of degrees of freedom is sufficiently large. In other words, the estimator and the actual error are almost the same. Thus, this estimator can be an alternative to the unknown error.

Figure 7

Comparisons of 𝐾-means with average criteria ( k = 6 ).

In this example, in order to compare the performance of marking strategies, several figures and tables are presented. For all figures, each marker denotes the value in level 5 ⁢ ℓ , ℓ = 0 , 1 , 2 , … . When the 𝑦-axis represents the time (s), the value indicates the computing time accumulated from the initial level. Adaptive algorithms are performed until the stopping criteria, e.g., η < 5 × 10 - 3 for k = 1 and η < 5 × 10 - 9 for k = 6 , are satisfied.

Firstly, the 𝐾-means clustering and average criteria are compared since these criteria are independent of the marking parameter θ M . Figure 7 shows how much time it takes for the estimator to reach the stopping criterion. From the numerical results, the average criterion takes around four times longer than the 𝐾-means clustering. The average criterion tends to mark large number of elements, including where the error is relatively small, and the number of degrees of freedom increases remarkably fast. On the other hand, the degrees of freedom for the 𝐾-means clustering is much less than that for the average criterion although the number of loops is larger. Thus, we can conclude that the 𝐾-means clustering criterion is much more efficient than the average criterion.

Figure 8
Comparison of 𝐾-means with maximum and bulk criteria (

k
=
6

k=6

).

Figure 8

Comparison of 𝐾-means with maximum and bulk criteria ( k = 6 ).

Figure 8 shows the performance of each marking strategy with k = 6 . From the figure, we can say that the 𝐾-means marking strategy is the most efficient for θ M = 0.1 . For other choices of θ M , the results vary depending on the numerical schemes. The 𝐾-means strategy is the best one for θ M = 0.1 and 0.3, but it shows quite similar performance to the maximum criterion for θ M = 0.5 and 0.7. The bulk strategy is the most efficient for θ M = 0.5 and 0.7, but computing time for three strategies is comparable. From the results, the rule of thumb is to choose θ M = 0.5 for the HDG methods. Figure 9 shows the results obtained from the linear approximation ( k = 1 ). The 𝐾-means strategy is more efficient than others when θ M is small, but the bulk strategy is marginally better than others when θ M is large. The rule of thumb is θ M = 0.7 for the HDG method.

Figure 9
Comparison of 𝐾-means with maximum and bulk criteria (

k
=
1

k=1

).

Figure 9

Comparison of 𝐾-means with maximum and bulk criteria ( k = 1 ).

Table 1

Correlation among variables for maximum and bulk criteria (HDG).

k = 6

	maximum				bulk

θ M	𝐿	dof	𝜂	time (s)	𝐿	dof	𝜂	time (s)
0.1	89	58751	4.48e−9	73.73	100	57974	4.58e−9	77.23
0.3	43	63441	3.69e−9	30.14	48	58695	4.74e−9	31.80
0.5	39	59605	4.77e−9	25.53	38	59157	4.74e−9	23.05
0.7	38	59150	4.76e−9	25.05	37	70840	3.59e−9	25.12

k = 1

	maximum				bulk

θ M	𝐿	dof	𝜂	time (s)	𝐿	dof	𝜂	time (s)
0.1	86	245282	4.97e−3	286.58	47	265754	4.77e−3	105.78
0.3	25	280258	4.64e−3	64.70	21	279012	4.66e−3	48.68
0.5	18	281860	4.63e−3	39.60	15	363648	4.11e−3	43.92
0.7	17	294758	4.55e−3	40.86	12	387338	4.24e−3	37.98

Table 1 represents the results for the HDG_D when the stopping criterion is satisfied, in terms of the number of loops (𝐿), degrees of freedom (dof), the value for the estimator (𝜂) and computing time (s) corresponding to θ M for the maximum and bulk criteria. Here, the result for the 𝐾-means criterion is omitted because it is independent of θ M . For this criterion, L = 39 , dof = 61509 , η = 4.18 ⁢ e - 09 , time = 25.24 with k = 6 , and L = 18 , dof = 335170 , η = 4.23 ⁢ e - 03 , time = 47.49 with k = 1 . From the numerical results, we notice that the computing time is significantly affected by the number of loops on the fine grid. Thus, for small computational cost, it is advantageous to reduce the number of loops on the fine grid in the adaptive algorithm. For the HDG methods with k = 6 , the choice of θ M = 0.5 shows better performance than other choices. When k = 1 , θ M = 0.7 is the best one for the bulk criterion.

From the results of Figures 8–9 and Table 1, we observe that the best choice of marking parameter θ M depends on the polynomial degree. However, there is no theoretical result for the optimal choice of the marking parameter. Thus, numerical experiments are performed with a proper choice of θ M . In this case, if we choose inappropriate θ M , it may take quite long to satisfy stopping criteria. Therefore, the 𝐾-means strategy is useful from the point of view of computation because it is independent of θ M , and it shows comparable results if it is not better.

Figure 10

Convergence history (left) and the effectivity indices (right).

Figure 10 shows the convergence rate and the effectivity index defined by η / | | | u - u h | | | . From the results, the optimal convergence rates with respect to the number of degrees of freedom are observed for all numerical schemes by adaptive algorithms with the 𝐾-means strategy. Note that the exact solution (4.1) has a corner singularity at the origin, and it is well known that the optimal rate of convergence is 1 3 with respect to the number of degrees of freedom on the uniform mesh. As the number of degrees of freedom is increasing, the effectivity index is close to 1. Thus, we can say that our estimator can be an alternative to the non-computable exact error.

Example 4.2

Consider the domain Ω = ( - 1 , 1 ) 2 divided into four subdomains Ω i , where

Ω 1 = ( 0 , 1 ) 2 , Ω 2 = ( - 1 , 0 ) × ( 0 , 1 ) , Ω 3 = ( - 1 , 0 ) 2 and Ω 4 = ( 0 , 1 ) × ( - 1 , 0 ) .

The diffusion coefficient 𝜌 is defined as ρ = 100 in Ω 1 , Ω 3 and ρ = 1 in Ω 2 , Ω 4 . The source term is f = 0 with the exact solution

(4.2) u ⁢ ( r , θ ) = r α ⁢ ( a i ⁢ sin ⁡ ( α ⁢ θ ) + b i ⁢ cos ⁡ ( α ⁢ θ ) ) ,

in each Ω i . Coefficients a i , b i depending on Ω i and 𝛼 are given as follows: α = 0.12690207 ,

a 1 = 0.1 , b 1 = 1 , a 2 = - 9.60396040 , b 2 = 2.96039604 ,

a 3 = - 0.48035487 , b 3 = - 0.88275659 , a 4 = 7.70156488 , b 4 = - 6.45646175 .

Figure 11
The exact solution (top), linear approximate solution (bottom left) and sixth-order approximate solution (bottom right) of Example 4.2 (N-HDG).

Figure 11

The exact solution (top), linear approximate solution (bottom left) and sixth-order approximate solution (bottom right) of Example 4.2 (N-HDG).

The exact solution has a singularity at the origin and relatively large changes on the common boundaries between neighboring subdomains. On the other hand, this solution is smooth on the remaining part. Thus, we can expect that mesh structure from adaptive algorithm has elements refined intensively near the origin and possibly near the common boundaries. In Figure 11, we display the exact solution and approximate solutions. We can see that sixth-order approximation produces better results than linear approximation, when their approximate solutions have almost the same number of degrees of freedom.

$Figure 12 Convergence history (left) and effectivity index (right) for GCBE method with θ M = 0.5 \theta_{M}=0.5 .$

Figure 12

Convergence history (left) and effectivity index (right) for GCBE method with θ M = 0.5 .

$Figure 13 Convergence history (left) and effectivity index (right) for N-HDG method with θ M = 0.5 \theta_{M}=0.5 .$

Figure 13

Convergence history (left) and effectivity index (right) for N-HDG method with θ M = 0.5 .

Figure 12 shows the convergence history for the residual type estimator for the GCBE method. Each marker represents the error in the level 10 ⁢ ℓ , ℓ = 0 , 1 , … . Similar to Example 4.1, the optimal convergence rates are recovered for various polynomial approximation k ≥ 2 . For effectivity indices, they increase as the polynomial order becomes larger. The results of the flux reconstruction based estimator for the N-HDG_D method are displayed in Figure 13. Here, the penalty parameter 𝜅 is chosen to be 1. The optimal convergence rates are observed regardless of the polynomial order. Also, effectivity indices are close to 1 when the number of degrees of freedom is sufficiently large. If the number of degrees of freedom is small, effectivity indices are larger than those in Example 4.1, but they are comparable if the number of degrees of freedom is sufficiently large.

In our (undisplayed) numerical results, we have adaptive meshes which are refined intensively near the origin for both the residual type and the flux reconstruction based estimators. The residual type estimator for the HDG methods depends on the penalty parameter 𝜅, but the flux reconstruction based estimator does not depend on it for Example 4.2. Also, performance of I-HDG_D and S-HDG_D with proper 𝜅 is comparable to that of the N-HDG_D method.

Figure 14

Comparisons of 𝐾-means with average criteria ( k = 6 ).

In this example, for Figures 14–16, each marker denotes the value in level 10 ⁢ ℓ , ℓ = 0 , 1 , 2 , … . As in Example 4.1, when the 𝑦-axis represents the time (s), the value indicates the computing time accumulated from the initial level. Adaptive algorithms are performed until the stopping criteria, e.g., η < 3 × 10 - 2 for k = 1 and η < 5 × 10 - 6 for k = 6 , are satisfied.

Figure 14 shows the results for the 𝐾-means and average criteria. Similar to Example 4.1, the average criterion takes around four times longer than the 𝐾-means strategy. Also, the number of degrees of freedom for the average criterion is around five times more than that for the 𝐾-means strategy. This is due to the small average because only few elements of the large number of elements have significantly large value of the error estimator.

Figure 15
Comparison of 𝐾-means with maximum and bulk criteria (

k
=
6

k=6

).

Figure 15

Comparison of 𝐾-means with maximum and bulk criteria ( k = 6 ).

Figure 16
Comparison of 𝐾-means with maximum and bulk criteria (

k
=
1

k=1

).

Figure 16

Comparison of 𝐾-means with maximum and bulk criteria ( k = 1 ).

Performance of three marking strategies is compared in Figures 15–16. Figure 15 shows the results obtained from the 𝐾-means, maximum and bulk criteria for k = 6 . For the HDG method, the 𝐾-means strategy is the most efficient for θ M = 0.1 and 0.5, but the bulk criterion is the best one for θ M = 0.3 . For θ M = 0.7 , three strategies show comparable results. Figure 16 shows the results for k = 1 . Here, the 𝐾-means strategy generally shows better performance than other criteria. From the results in Figures 15–16, we observe that the best choice of θ M is 0.3 for the HDG method, while θ M = 0.5 for k = 6 and θ M = 0.7 for k = 1 in Example 4.1.

Figure 17 shows the rate of convergence and the effectivity index for the HDG method. Note that the exact solution (4.2) has a singularity at the origin, and the regularity of the solution is low. Also, relatively large changes arise on the common boundaries between neighboring subdomains because of the large difference of the diffusion coefficient 𝜌. Similar to Example 4.1, the optimal convergence rates with respect to the number of degrees of freedom are observed for all numerical schemes by adaptive algorithms with the 𝐾-means strategy. As running adaptive algorithm, the fact that effectivity goes to 1 tells us the error estimator 𝜂 can be an alternative to exact energy error when the exact solution is unknown.

Figure 17

Convergence history (left) and the effectivity indices (right).

5 Conclusions

In this paper, we have developed the residual type and the flux reconstruction based a posteriori error estimates for a class of HDG methods. The estimators are fully computable, reliable and efficient modulo a data-oscillation term. Moreover, they provide optimal rates of convergence in terms of the number of degrees of freedom for any given HDG approximations. Next, we have proposed a new marking strategy which is based on the 𝐾-means clustering. As can be seen in the Figures 8, 9, 15, 16 and Table 1, the choice of marking parameter is crucial in the performance of the existing strategy such as maximum and bulk criteria; however, the optimal choice is not known a priori. The new strategy has no unknown parameters. Our numerical experiments show that the proposed strategy outperforms the existing ones for low-order polynomials. For high-order cases, the proposed strategy still shows better performance than the existing marking criteria in general.

We summarize our findings of hybrid discontinuous Galerkin methods on various perspectives. At first, the S-HDG methods are stable with sufficiently large 𝜅, while the N-HDG methods are stable regardless of 𝜅. The S-HDG methods are more effective in terms of solver (cf. [40]). Second, the HDG methods with continuous traces can reduce lots of degrees of freedom, but the average flux continuity does not hold. The optimal convergence rate is recovered with the residual type error estimator. Note that its effectivity index depends on the penalty parameter 𝜅. Since the GCBE method, a special case of the N-HDG_C method, has no penalty term, the residual type error estimator is simple and robust in terms of 𝜅. Third, the HDG methods with discontinuous traces satisfy the average flux continuity. The flux reconstruction based error estimator can be applied to these methods. Although computational cost is higher than the residual type estimator, we can achieve the effectivity index close to 1 when the number of degrees of freedom is sufficiently large. Lastly, high-order methods have advantages in comparison to low-order methods. We can achieve high-order accuracy with high-order methods even for nonsmooth solutions via adaptivity.

Funding source: National Research Foundation of Korea

Award Identifier / Grant number: NRF-2015R1A5A1009350

Award Identifier / Grant number: NRF-2019R1A2C2090021

Funding statement: The first author was supported by 2020 Research Fund of Myongji University. The third author was supported by NRF-2015R1A5A1009350 and 2019R1A2C2090021.

References

[1] M. Ainsworth and J. T. Oden, A Posteriori Error Estimation in Finite Element Analysis, Pure Appl. Math. (New York), Wiley-Interscience, New York, 2000. 10.1002/9781118032824Search in Google Scholar

[2] M. Ainsworth and R. Rankin, Fully computable bounds for the error in nonconforming finite element approximations of arbitrary order on triangular elements, SIAM J. Numer. Anal. 46 (2008), no. 6, 3207–3232. 10.1137/07070838XSearch in Google Scholar

[3] I. Babuška and T. Strouboulis, The Finite Element Method and its Reliability, Oxford University, New York, 2001. 10.1093/oso/9780198502760.001.0001Search in Google Scholar

[4] M. Bebendorf, A note on the Poincaré inequality for convex domains, Z. Anal. Anwendungen 22 (2003), no. 4, 751–756. 10.4171/ZAA/1170Search in Google Scholar

[5] Y. Bengio, A. C. Courville and P. Vincent, Unsupervised feature learning and deep learning: A review and new perspectives, preprint (2012), https://arxiv.org/abs/1206.5538. Search in Google Scholar

[6] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer Ser. Comput. Math. 15, Springer, New York, 1991. 10.1007/978-1-4612-3172-1Search in Google Scholar

[7] E. Burman and A. Ern, Continuous interior penalty h ⁢ p -finite element methods for advection and advection-diffusion equations, Math. Comp. 76 (2007), no. 259, 1119–1140. 10.1090/S0025-5718-07-01951-5Search in Google Scholar

[8] Z. Cai, V. Carey, J. Ku and E.-J. Park, Asymptotically exact a posteriori error estimators for first-order div least-squares methods in local and global L2 norm, Comput. Math. Appl. 70 (2015), 648–659. 10.1016/j.camwa.2015.05.010Search in Google Scholar

[9] C. Carstensen and S. A. Funken, Constants in Clément-interpolation error and residual based a posteriori error estimates in finite element methods, East-West J. Numer. Math. 8 (2000), no. 3, 153–175. Search in Google Scholar

[10] C. Carstensen, D. Gallistl and J. Hu, A posteriori error estimates for nonconforming finite element methods for fourth-order problems on rectangles, Numer. Math. 124 (2013), no. 2, 309–335. 10.1007/s00211-012-0513-5Search in Google Scholar

[11] C. Carstensen, J. Gedicke and E.-J. Park, Numerical experiments for the Arnold–Winther mixed finite elements for the Stokes problem, SIAM J. Sci. Comput. 34 (2012), no. 4, A2267–A2287. 10.1137/100802906Search in Google Scholar

[12] C. Carstensen, D. Kim and E.-J. Park, A priori and a posteriori pseudostress-velocity mixed finite element error analysis for the Stokes problem, SIAM J. Numer. Anal. 49 (2011), no. 6, 2501–2523. 10.1137/100816237Search in Google Scholar

[13] C. Carstensen, N. Nataraj and A. K. Pani, Comparison results and unified analysis for first-order finite volume element methods for a Poisson model problem, IMA J. Numer. Anal. 36 (2016), no. 3, 1120–1142. 10.1093/imanum/drv050Search in Google Scholar

[14] C. Carstensen and E.-J. Park, Convergence and optimality of adaptive least squares finite element methods, SIAM J. Numer. Anal. 53 (2015), no. 1, 43–62. 10.1137/130949634Search in Google Scholar

[15] C. Carstensen, E.-J. Park and P. Bringmann, Convergence of natural adaptive least squares finite element methods, Numer. Math. 136 (2017), no. 4, 1097–1115. 10.1007/s00211-017-0866-xSearch in Google Scholar

[16] E. T. Chung, E.-J. Park and L. Zhao, Guaranteed a posteriori error estimates for a staggered discontinuous Galerkin method, J. Sci. Comput. 75 (2018), no. 2, 1079–1101. 10.1007/s10915-017-0575-8Search in Google Scholar

[17] B. Cockburn, J. Gopalakrishnan and R. Lazarov, Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems, SIAM J. Numer. Anal. 47 (2009), no. 2, 1319–1365. 10.1137/070706616Search in Google Scholar

[18] E. Dari, R. Duran, C. Padra and V. Vampa, A posteriori error estimators for nonconforming finite element methods, RAIRO Modél. Math. Anal. Numér. 30 (1996), no. 4, 385–400. 10.1051/m2an/1996300403851Search in Google Scholar

[19] A. Ern, S. Nicaise and M. Vohralík, An accurate H ⁢ ( div ) flux reconstruction for discontinuous Galerkin approximations of elliptic problems, C. R. Math. Acad. Sci. Paris 345 (2007), no. 12, 709–712. 10.1016/j.crma.2007.10.036Search in Google Scholar

[20] A. Ern and M. Vohralík, A unified framework for a posteriori error estimation in elliptic and parabolic problems with application to finite volumes, Finite Volumes for Complex Applications VI. Problems & Perspectives. Vol. 1, 2, Springer Proc. Math. 4, Springer, Heidelberg (2011), 821–837. 10.1007/978-3-642-20671-9_85Search in Google Scholar

[21] R. Eymard, T. Gallouët and R. Herbin, Finite volume methods, Handbook of Numerical Analysis. Vol. VII, North-Holland, Amsterdam (2000), 713–1020. 10.1016/S1570-8659(00)07005-8Search in Google Scholar

[22] J. A. Hartigan and M. A. Wong, Algorithm AS 136: A k-means clustering algorithm, J. Roy. Statist. Soc. Ser. C 28 (1979), 100–108. 10.2307/2346830Search in Google Scholar

[23] A. K. Jain, M. N. Murty and P. J. Flynn, Data clustering: A review, ACM Comput. Surv. (CSUR) 31 (1999), 264–323. 10.1145/331499.331504Search in Google Scholar

[24] Y. Jeon and E.-J. Park, Nonconforming cell boundary element methods for elliptic problems on triangular mesh, Appl. Numer. Math. 58 (2008), no. 6, 800–814. 10.1016/j.apnum.2007.02.010Search in Google Scholar

[25] Y. Jeon and E.-J. Park, A hybrid discontinuous Galerkin method for elliptic problems, SIAM J. Numer. Anal. 48 (2010), no. 5, 1968–1983. 10.1137/090755102Search in Google Scholar

[26] Y. Jeon and E.-J. Park, New locally conservative finite element methods on a rectangular mesh, Numer. Math. 123 (2013), no. 1, 97–119. 10.1007/s00211-012-0477-5Search in Google Scholar

[27] Y. Jeon, E.-J. Park and D. Sheen, A hybridized finite element method for the Stokes problem, Comput. Math. Appl. 68 (2014), no. 12, 2222–2232. 10.1016/j.camwa.2014.08.005Search in Google Scholar

[28] O. A. Karakashian and F. Pascal, A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems, SIAM J. Numer. Anal. 41 (2003), no. 6, 2374–2399. 10.1137/S0036142902405217Search in Google Scholar

[29] D. Kim and E.-J. Park, A posteriori error estimator for expanded mixed hybrid methods, Numer. Methods Partial Differential Equations 23 (2007), no. 2, 330–349. 10.1002/num.20178Search in Google Scholar

[30] D. Kim and E.-J. Park, A posteriori error estimators for the upstream weighting mixed methods for convection diffusion problems, Comput. Methods Appl. Mech. Engrg. 197 (2008), no. 6–8, 806–820. 10.1016/j.cma.2007.09.009Search in Google Scholar

[31] D. Kim and E.-J. Park, A priori and a posteriori analysis of mixed finite element methods for nonlinear elliptic equations, SIAM J. Numer. Anal. 48 (2010), no. 3, 1186–1207. 10.1137/090747002Search in Google Scholar

[32] K. Y. Kim, A posteriori error analysis for locally conservative mixed methods, Math. Comp. 76 (2007), no. 257, 43–66. 10.1090/S0025-5718-06-01903-XSearch in Google Scholar

[33] K. Y. Kim, A posteriori error estimators for locally conservative methods of nonlinear elliptic problems, Appl. Numer. Math. 57 (2007), no. 9, 1065–1080. 10.1016/j.apnum.2006.09.010Search in Google Scholar

[34] K.-Y. Kim, Guaranteed a posteriori error estimator for mixed finite element methods of elliptic problems, Appl. Math. Comput. 218 (2012), no. 24, 11820–11831. 10.1016/j.amc.2012.04.084Search in Google Scholar

[35] J. MacQueen, Some methods for classification and analysis of multivariate observations, Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability (Berkeley 1965/66), University of California, Berkeley (1967), 281–297. Search in Google Scholar

[36] L. E. Payne and H. F. Weinberger, An optimal Poincaré inequality for convex domains, Arch. Rational Mech. Anal. 5 (1960), 286–292. 10.1007/BF00252910Search in Google Scholar

[37] W. Prager and J. L. Synge, Approximations in elasticity based on the concept of function space, Quart. Appl. Math. 5 (1947), 241–269. 10.1090/qam/25902Search in Google Scholar

[38] S. Repin, A Posteriori Estimates for Partial Differential Equations, Radon Ser. Comput. Appl. Math. 4, Walter de Gruyter, Berlin, 2008. 10.1515/9783110203042Search in Google Scholar

[39] J. E. Roberts and J.-M. Thomas, Mixed and hybrid methods, Handbook of Numerical Analysis. Vol. II, North-Holland, Amsterdam (1991), 523–639. 10.1016/S1570-8659(05)80041-9Search in Google Scholar

[40] D.-W. Shin, Y. Jeon and E.-J. Park, A hybrid discontinuous Galerkin method for advection-diffusion-reaction problems, Appl. Numer. Math. 95 (2015), 292–303. 10.1016/j.apnum.2014.11.003Search in Google Scholar

[41] R. Stenberg, A family of mixed finite elements for the elasticity problem, Numer. Math. 53 (1988), no. 5, 513–538. 10.1007/BF01397550Search in Google Scholar

[42] J. L. Synge, The Hypercircle in Mathematical Physics: A Method for the Approximate Solution of Boundary Value Problems, Cambridge University, New York, 1957. 10.1063/1.3060143Search in Google Scholar

[43] R. Verfürth, A Review of a Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques, Teubner-Wiley, Stuttgart, 1996. Search in Google Scholar

[44] M. Vohralík, Two types of guaranteed (and robust) a posteriori estimates for finite volume methods, Finite Volumes for Complex Applications V, ISTE, London (2008), 649–656. Search in Google Scholar

[45] G. N. Wells, Analysis of an interface stabilized finite element method: The advection-diffusion-reaction equation, SIAM J. Numer. Anal. 49 (2011), no. 1, 87–109. 10.1137/090775464Search in Google Scholar

[46] L. Zhao and E.-J. Park, Fully computable bounds for a staggered discontinuous Galerkin method for the Stokes equations, Comput. Math. Appl. 75 (2018), no. 11, 4115–4134. 10.1016/j.camwa.2018.03.018Search in Google Scholar

[47] L. Zhao, E.-J. Park and D.-W. Shin, A staggered DG method of minimal dimension for the Stokes equations on general meshes, Comput. Methods Appl. Mech. Engrg. 345 (2019), 854–875. 10.1016/j.cma.2018.11.016Search in Google Scholar

Received: 2021-06-02

Accepted: 2021-06-08

Published Online: 2021-06-29

Published in Print: 2021-10-01

Novel Adaptive Hybrid Discontinuous Galerkin Algorithms for Elliptic Problems

Abstract

1 Introduction

2 Methods

Lemma 2.1 (Weak Solution)

Proof

Proof

3 A Posteriori Error Estimate

Definition 3.1 (Potential Reconstruction)

3.1 A Residual Type A Posteriori Error Estimate

Theorem 3.2 (A Residual Type A Posteriori Error Estimate)

Proof

Lemma 3.3 (Averaging Operator)

Lemma 3.4 (Efficiency of the Element Residuals)

Theorem 3.5 (Local Efficiency of Estimate (3.5))

Proof

Corollary 3.6 (Local Efficiency of Estimate (3.10))

3.2 Flux Reconstruction Based on the A Posteriori Error Estimate

Definition 3.7 (Flux Reconstruction)

Proposition 3.8 (Flux Reconstruction for the HDG Method)

Theorem 3.9 (A Posteriori Error Estimate Based on the Flux Reconstruction)

Proof

Theorem 3.10 (Efficiency of the Estimate of Theorem 3.9)

Proof

4 Numerical Experiments

Algorithm (Adaptive Algorithm)

Algorithm (𝐾-Means Clustering)

5 Conclusions

References

Journal and Issue

Articles in the same Issue