    We consider an adaptive finite element method with arbitrary but fixed polynomial degree p > 1, where adaptivity is driven by an edge-based residual error estimator. Based on the modified maximum criterion from [Diening et al, Found. Comput. Math. 16, 2016], we propose a goal-oriented adaptive algorithm and prove that it is instance optimal. Numerical experiments underline our theoretical findings.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

### 1.1. Rate optimality vs. instance optimality of AFEM

For an elliptic PDE with sought solution , the adaptive finite element method (AFEM) iterates the loop

 (1) \textsl\slSOLVE⟶\textslESTIMATE⟶\textsl\slMARK⟶\textsl\slREFINE

to successively adapt an initial mesh to the possible singularities of the sought solution . This leads to a sequence of meshes and corresponding discrete solutions being -piecewise polynomials of degree . In the last decades, the mathematical understanding of AFEM has matured. We refer to [Dör96, MNS00, BDD04, Ste07, CKNS08, FFP14] for some milestones of the analysis of convergence with optimal algebraic rates, to [MSV08, Sie11] for a general theory on (plain) convergence, and to [CFPP14] for an abstract approach for optimal convergence rates.

The works [Ste07, CKNS08, FFP14, CFPP14] aim for rate optimality of adaptive algorithms, i.e., they prove that the usual adaptive loop (1), based on (quasi-) minimal Dörfler marking [Dör96] for the step MARK, leads to optimal decay of the total error

 (2) error(TH):=|||u−uH|||+oscH(TH),

being the sum of energy error plus data oscillation terms (or, equivalently, the error estimator, since ). Let denote the number of elements of a finite set. Then, rate optimality reads as follows: For all , it holds that

 (3) ∥u∥As<∞⟹∃C>0∀ℓ∈N0:error(Tℓ)≤C(#Tℓ)−s,

where corresponds to certain approximation properties of (which can be characterized in terms of Besov regularity [BDDP02, GM14, Gan17]).

Unlike this, the first work [BDD04] on AFEM with convergence rates implicitly proved even instance optimality of AFEM, i.e., the total error on an adaptive mesh is quasi-optimal with respect to all refinements of , which have essentially the same number of elements: It holds that

 (4) ∃C>1∀ℓ∈N0∀TH∈refine(T0):[C#TH≤#Tℓ⟹error(Tℓ)≤Cerror(TH)].

The key argument for the proof of (4) in [BDD04] is an additional coarsening step in the adaptive loop (1). While [BDD04] also employed the Dörfler marking criterion [Dör96] for the step MARK, the work [DKS16] proposed a modified maximum criterion to single out edges for refinement. For P1-AFEM for the 2D Poisson problem, [DKS16] then proved instance optimality (4) of their adaptive strategy without resorting to an additional coarsening step. In [KS16], their analysis was extended to AFEM with non-conforming P1 elements for the Poisson problem and the Stokes system in 2D. We stress that any instance-optimal AFEM (4) is, in particular, rate optimal (3).

While standard adaptivity aims to approximate the PDE solution by some discrete approximation in the energy norm, a goal-oriented adaptive finite element method (GOAFEM) aims only to approximate by , where is the so-called goal functional or quantity of interest.

In the present paper, we consider the following problem: Let be a polygonal Lipschitz domain, which is resolved by the initial mesh , where is admissible in the sense of [BDD04, Ste08]. Given and , the (linear) goal functional reads

 (5) G(u):=∫Ωgu−g⋅∇udx,

where is the unique solution to

 (6a) −div(A∇u) =f+divf in Ω, (6b) u =0 on Γ\coloneqq∂Ω.

For technical reasons, we assume that the diffusion matrix is -piecewise constant and is symmetric and positive definite. Moreover, we assume that the restrictions are smooth for all .

Convergence and rate-optimality of GOAFEM has been addressed in [MS09, BET11, HPZ15, HP16, FGH16, FPZ16]. The key idea of the argument is to let be the unique solution to the dual problem

 (7a) −div(A∇u∗) =g+divg in Ω, (7b) u∗ =0 on Γ.

Throughout, quantities associated with the dual problem are indexed by an asterisk. We note that the (symmetric) primal problem (6) and the dual problem (7) coincide up to the right-hand side. Let be the problem-induced energy norm. For FEM approximations and , standard duality arguments (together with the Galerkin orthogonality) lead to

 (8) |G(u)−G(uℓ)|≤|||u−uℓ||||||u∗−u∗ℓ|||≤[|||u−uℓ|||+oscℓ(Tℓ)][|||u∗−u∗ℓ|||+osc∗ℓ(Tℓ)];

see, e.g., [MS09, FPZ16]. Therefore, GOAFEM aims to control and steer the product of the total errors

 (9) error(TH)error∗(TH):=[|||u−uH|||+oscH(TH)][|||u∗−u∗H|||+osc∗H(TH)].

While [BET11, HPZ15, HP16] focus on linear convergence of GOAFEM, the works [MS09, FGH16, FPZ16] also prove rate optimality. All works employ variants of the Dörfler marking criterion [Dör96]: The seminal work [MS09] employs (quasi-) minimal Dörfler marking separately for the primal and the dual problem, and then uses the smaller set for MARK. Instead, [BET11] proposes a (quasi-) minimal Dörfler marking for some combined estimator. Both strategies guarantee rate optimality for the product of the total errors

 (10) ∥u∥As+∥u∗∥At<∞⟹∃C>0∀ℓ∈N0:error(Tℓ)error∗(Tℓ)≤C(#Tℓ)−(s+t),

where the possible algebraic rate now depends on the approximability properties of the primal and dual problem; see [MS09, FGH16, FPZ16].

### 1.3. Instance-optimal GOAFEM

The new GOAFEM algorithm can briefly be outlined as follows: SOLVE computes the FEM solution to the primal problem (6) and to the dual problem (7). ESTIMATE computes the corresponding residual error estimators and . MARK employs the modified maximum strategy from [DKS16] to obtain two sets of marked edges, namely with respect to and with respect to . With , we then define , where and are arbitrary up to . Finally, REFINE employs 2D newest vertex bisection (NVB) to generate the coarsest mesh , where all edges in have been bisected once.

The main result of the present work states that the proposed GOAFEM is instance optimal (4) with respect to the total-error product, i.e.,

 (11) ∃C>1∀ℓ∈N0∀TH,TH∗∈refine(T0): \mathclap[C max{#TH,#TH∗}≤#Tℓ⟹error(Tℓ)error∗(Tℓ)≤Cerror(TH)error∗(TH∗)].

Again, we note that this implies, in particular, rate optimality (3). On a technical side, we note that the seminal work [DKS16] is restricted to the lowest-order FEM discretization , while the present analysis also allows higher (but fixed) polynomial degrees . In this sense, the present work provides also the technical tools to generalize the instance-optimal AFEM of [DKS16] from to general but fixed .

### 1.4. Outline

The remainder of this paper is organized as follows: In Section 2, we give a precise formulation of the modules SOLVE, ESTIMATE, MARK, and REFINE of the adaptive loop (1). In particular, we state the modified maximum criterion (Algorithm 2.4) from [DKS16] as well as our extension to GOAFEM (Algorithm 2.4). Then, we thoroughly formulate our GOAFEM algorithm (Algorithm 2.6) and state our main result that the proposed GOAFEM is instance optimal (Theorem 2.6). Section 3 reviews the proof of the seminal work [DKS16] in an abstract framework. Section 4 collects the technical results to generalize the seminal work [DKS16] from lowest-order FEM to arbitrary polynomial degree (Theorem 2.5). Thereafter, Section 5 gives the proof of Theorem 2.6. Some numerical experiments in Section 6 conclude this work and empirically compare the instance-optimal GOAFEM algorithm from the present work with the rate-optimal GOAFEM strategies from [MS09, BET11, FPZ16].

### 1.5. General notation

In all results, the involved constants (as well as their dependencies) are stated explicitly. In proofs, however, we write to abbreviate up to a multiplicative constant which is clear from the context. Moreover, we write if both estimates, and , hold.

## 2. Main result

Before stating our main result, we discuss the particular modules of the adaptive loop (1) and fix the necessary notation.

### 2.1. Refine

A mesh is a conforming triangulation of into non-degenerate compact triangles . The edges of are denoted by . The set of interior edges of is denoted by , i.e., for each , there exist unique such that . The set of vertices of is denoted by . We define the patches

 (12) TH(ω):={T∈TH:T∩ω≠∅}for all ω⊂¯¯¯¯Ω.

For vertices , we abbreviate . For neighbors , we also consider the reduced edge patch

 (13) TredH(E):={T∈TH:E⊂∂T}={T,T′}for E=T∩T′∈EH.

Similarly, we define the area associated to a set of triangles by

 (14) Ω(UH):=⋃T∈UHT⊆¯¯¯¯Ωfor all UH⊆TH.

For mesh-refinement, we employ an edge-based variant of newest vertex bisection (NVB) [Ste08]. We suppose that the initial mesh is admissible in the sense of [BDD04, Ste08]: For all neighbors , the joint edge is the reference edge of if and only if it is also the reference edge of . While this assumption is unnecessary for the NVB algorithm [KPP13], it provides additional structure which is crucial in the instance-optimality analysis of [DKS16].

For a mesh and a set of marked edges, let be the coarsest NVB refinement of such that all edges have been bisected. Moreover, we write , if can be obtained by finitely many steps of NVB refinement. Then, is the set of all possible NVB refinements of . We note that NVB leads to uniformly shape-regular meshes in the sense of

 (15) Cmesh:=supTH∈TmaxT∈THdiam(T)2|T|<∞,

where is the area of a triangle .

### 2.2. Solve

As usual, the primal problem (6) is understood in weak form. The Lax–Milgram lemma guarantees existence and uniqueness of such that

 (16) a(u,v):=∫ΩA∇u⋅∇vdx=∫Ωfv−f⋅∇vdx=:F(v)for all v∈H10(Ω).

We define the energy norm and note that on . Let and . For the discretization of (16), define the space of -piecewise polynomials

 (17) Pp(TH):={v∈L2(Ω):∀T∈THv|T is a polynomial of degree ≤p}

as well as the conforming FEM spaces

 (18) \SSp(TH)\coloneqqPp(TH)∩C(Ω)=Pp(TH)∩H1(Ω)% and\SSp0(TH)\coloneqqSp(TH)∩H10(Ω).

Again, the Lax–Milgram lemma proves the existence and uniqueness of such that

 (19) a(uH,vH)=F(vH)for all vH∈\SSp0(TH).

The dual problem (7) is solved analogously with in (16) and (19) being replaced by from (5). The Lax–Milgram lemma guarantees existence and uniqueness of the dual solution and its FEM approximation .

### 2.3. Estimate

For a posteriori error estimation, we employ an edge-based residual error estimator. Let be the length of an edge . For the primal problem (6) with , we define

 (20) ηH(E)2\coloneqq|E|∥[[(A∇uH+f)⋅ν]]∥2L2(E)+∑T∈TredH(E)|T|∥f+div(A∇uH+f)∥2L2(T)% for all E∈EH,

where

is a normal vector on

and denotes the jump across . For the dual problem (7) with , we define analogously

 (21) η∗H(E)2\coloneqq|E|∥[[(A∇u∗H+g)⋅ν]]∥2L2(E)+∑T∈TredH(E)|T|∥g+div(A∇u∗H+g)∥2L2(T)for all E∈EH.

With this notation, we define

 (22) ηH(UH):=(∑E∈UHηH(E)2)1/2and η∗H(UH):=(∑E∈UHη∗H(E)2)1/2for all UH⊆EH.

With the -orthogonal projections and , where , the data oscillations read

 (23) oscH(T)2\coloneqq|T|∥(1−ΠT)(f+divf)∥2L2(T)+∑E∈EHE⊂∂T|T|1/2∥(1−ΠE)[[f⋅ν]]∥2L2(E).

Note that for the volume term simply reads . The oscillations for the dual problem are defined analogously with and instead of and , respectively. We note that

 (24) oscH(T)≲ηH(E) and osc∗H(T)≲η∗H(E)for all T∈TH,E∈EH with E⊆∂T,

where the hidden constant depends only on from (15). For a subset , we define and analogously to (22). We note that

 (25) C−1rel|||u−uH||| ≤ηH(EH)≤Ceff[|||u−uH|||+oscH(TH)], (26) C−1rel|||u∗−u∗H|||

where the reliability constant depends only on from (15), while the efficiency constant depends additionally on .

### 2.4. Mark

Let . We define the tail of an edge by

 (27) tailH(E)\coloneqqEH∖EH,E,whereTH,E:=refine(TH,{E}),

i.e., the tail consists of all edges, which have to be refined to ensure conformity of the triangulation if is bisected. Moreover, we define

 (28) tailH(UH):=⋃E∈UHtailH(E)for all UH⊆EH.

With these definitions, we recall the following modified maximum criterion from [DKS16], which leads to an instance-optimal AFEM.

Algorithm 1 (Modified maximum criterion).
Input:
Edges with indicators , marking parameter .
Output: Set of marked edges.

1:   and
2:
3:  while  do
4:     pick and update
5:     compute
6:     if  then
7:
8:     end if
9:  end while

We refer to [DKS16] for a recursive implementation of Algorithm 2.4, which has linear costs. In case of GOAFEM, we need a different marking step that takes the primal as well as the dual problem into account. To this end, recall the Gauss brackets for . We note that the following algorithm is slightly more general than the strategy outlined in Section 1.3 of the introduction (where ).

Algorithm 2 (Modified maximum criterion for GOAFEM).
Input:
Edges , indicators and , marking parameters and .
Output: Set of marked edges.

1:  generate by Algorithm 2.4
2:  generate by Algorithm 2.4
3:  choose and
4:  define
5:  pick with
6:  choose

### 2.5. Instance-optimal AFEM

The work [DKS16] analyzes the following instance of the adaptive loop (1), which turns out to be instance-optimal; see Theorem 2.5.

Algorithm 3 (Instance-optimal AFEM).
Input:
Initial mesh , polynomial degree , marking parameter .
Output: Meshes , discrete solutions , and estimators for all .

1:  for all  do
2:     compute FEM solution
3:     compute indicators
4:     generate by Algorithm 2.4
5:     employ NVB to generate
6:  end for

For , the following theorem is the main result of [DKS16]. Our analysis below implies that the result remains true for arbitrary polynomial degrees .

Theorem 4.  Let the initial mesh be admissible in the sense of [BDD04]. Let and . Then, the AFEM Algorithm 2.5 for the primal problem (6) is instance optimal with respect to the total error, i.e.,

 (29) ∃C>1∀ℓ∈N0∀TH∈refine(T0): \mathclap(C# (TH∖T0)≤#(Tℓ∖T0)⟹|||u−uℓ|||2+oscℓ(Tℓ)2≤C[|||u−uH|||2+oscH(TH)2]).

The constant depends only on , , , and the data , , .

Remark 5. We note that elementary calculation shows that, for all ,

 #TH−#T0≤#(TH∖T0)≤#TH≤(#T0)(#TH−#T0+1)≤2(#T0)(#TH−#T0);

see, e.g., [BHP17, Lemma 22]. Hence, in (29) can, in fact, be replaced by (at the cost that the constant in (4) will additionally depend on ). Therefore, the statement of Theorem 2.5 is equivalent to the introductory statement of instance optimality (4) in Section 1.1.

### 2.6. Instance-optimal GOAFEM

As outlined in the introduction, the main idea behind GOAFEM is the duality-based estimate

 |G(u)−G(uℓ)| =|a(u−uℓ,u∗)| =|a(u−uℓ,u∗−u∗ℓ)|≤ |||u−uℓ||||||u∗−u∗ℓ|||,

The formal statement of our GOAFEM algorithm reads as follows:

Algorithm 6 (Instance-optimal GOAFEM).
Input:
Initial mesh , polynomial degree , marking parameters and .
Output: Meshes , discrete solutions , estimators and goal quantities for all .

1:  for all  do
2:     compute FEM solutions and
3:     compute indicators and
4:     generate by Algorithm 2.4
5:     employ NVB to generate
6:  end for

The following theorem is the main result of this work. We stress that the theorem involves the adaptively generated mesh for the primal and the dual error and compares it with arbitrary meshes and , where is used for the primal error and is used for the dual error.

Theorem 7.  Let the initial mesh be admissible in the sense of [BDD04]. Let and as well as . Let be the sequence of meshes generated by Algorithm 2.6. Then, the AFEM Algorithm 2.5 is instance optimal with respect to the product of total errors, i.e.,

 (30) (Cmax{#(TH∖T0),#(TH∗∖T0)}≤#(Tℓ∖T0)⟹[|||u−uℓ|||2+oscℓ(Tℓ)2][|||u∗−u∗ℓ|||2+osc∗ℓ(Tℓ)2]≤C[|||u−uH|||2+oscH(Tℓ)2][|||u∗−u∗H∗|||2+osc∗H∗(Tℓ)2]).

The constant depends only on , , , and the data , , , , .

Remark 8. Note that the natural statement of instance-optimality for GOAFEM in the sense of (4) and (9) would be: such that

 (C#(TH∖T0)≤#(Tℓ∖T0) ⟹[|||u−uℓ|||2+oscℓ(Tℓ)2][|||u∗−u∗ℓ|||2+osc∗ℓ(Tℓ)2] ≤C[|||u−uH|||2+oscH(Tℓ)2][|||u∗−u∗H|||2+osc∗H(Tℓ)2]).

Our Theorem 2.6, however, is stronger. There, the mesh for the right-hand side can be chosen for both factors independently.

## 3. Abstract result on instance optimality

This section aims to review the proof of [DKS16, Theorem 7.3] in an abstract framework. For arbitrary , the tuple is a diamond, if

• are meshes for all

• with finest common coarsening and coarsest common refinement

• such that the areas are pairwise disjoint for all .

Note that exist (and are unique), since newest vertex bisection is a a binary refinement rule, where the order of refinements does not matter. This allows to write

 TH =m⋃j=1{T∈Tj:∀k∈{1,…,m}∀T′∈Tk(T⊆T′⟹T=T′)}, Th =m⋃j=1{T∈Tj:∀k∈{1,…,m}∀T′∈Tk(T′⊆T⟹T=T′)}.

Moreover, according to [DKS16], this also implies that

 (31) TH∖Th=m⋃j=1(Tj∖Th).

Diamonds are a means to couple the lattice structure of with an abstract energy

 (32) E:T→R≥0.

Only energies that are compatible with this structure are suitable to prove instance optimality. This is encoded in the following properties 15, where , , , are generic constants, is a computable edge-based estimator, and mark is an abstract marking strategy:

1. [label=(A0)]

2. Marking criterion: For all meshes with edges , the marking strategy guarantees that the marked edges satisfy that

 MH≠∅andη2H(tailH(MH))≥Cmark(#MH)maxE∈EHη2H(tailH(E)).
3. Monotonicity of energy: For all and all , it holds that

 0≤E(Th)≤E(TH).
4. Diamond estimate: For all diamonds