NISTIR 7668

A Collection of 2D Elliptic

Problems for Testing

Adaptive Algorithms

William F. Mitchell

U. S. Department of Commerce

Technology Administration

National Institute of Standards and Technology

Information Technology Laboratory

Gaithersburg, MD 20899 USA

February 2010

A Collection of 2D Elliptic Problems for Testing

Adaptive Algorithms

William F. Mitchell∗

Mathematical and Computational Sciences Division

National Institute of Standards and Technology

Gaithersburg, MD 20899-8910

Abstract

Adaptive grid refinement is a critical component of the improvements

that have recently been made in algorithms for the numerical solution of

partial differential equations (PDEs). The development of new algorithms

and computer codes for the solution of PDEs usually involves the use of

proof-of-concept test problems. 2D elliptic problems are often used as

the first test bed for new algorithms and codes. This paper contains

a set of twelve parametrized 2D elliptic test problems for adaptive grid

refinement algorithms and codes. The problems exhibit a variety of types

of singularities, near singularities, and other difficulties.

Keywords: adaptive grid refinement, elliptic partial differential equations,

test problems

1 Introduction

The numerical solution of partial differential equations (PDEs) is the most

compute-intensive part of a wide range of scientific and engineering applications. Consequently the development and application of faster and more accurate methods for solving partial differential equations has received much attention in the past fifty years. Self-adaptive methods to determine a quasi-optimal

grid are a critical component of the improvements. Although adaptive grid refinement techniques are now in widespread use in applications, they remain an

active field of research, particularly in the context of hp-adaptive techniques.

The development of new algorithms and computer codes for the solution

of PDEs usually involves the use of proof-of-concept test problems. Such test

problems have a variety of uses such as demonstrating that a new algorithm

is effective, verifying that a new code is correct in the sense of achieving the

theoretical order of convergence, and comparing the performance of different

∗Contribution of NIST, not subject to copyright.

1
algorithms and codes. Nearly every paper on algorithms for solving PDEs contains a numerical results section with one or more test problems.

2D elliptic problems are often used as the first test bed for new algorithms

and codes for solving PDEs. This paper contains a set of twelve 2D elliptic test

problems for adaptive grid refinement algorithms and codes. Most of the problems are taken from the numerical results section of papers in the adaptive grid

refinement literature. The problems exhibit a variety of types of singularities

(point and line singularities on the boundary and in the interior), near singularities (sharp peaks, boundary layers, and wave fronts), and other difficulties.

Most of the problems are parametrized to allow “easy” and “hard” variations

on the problem.

We primarily consider elliptic partial differential equations of the form

−

∂

∂x 

p(x, y)

∂u

∂x

−

∂

∂y 

q(x, y)

∂u

∂y 

+ r(x, y)u = f(x, y) in Ω (1)

u = gD(x, y) on ∂ΩD (2)

p(x, y)

∂u

∂x

∂y

∂s − q(x, y)

∂u

∂y

∂x

∂s + c(x, y)u = gN (x, y) on ∂ΩN (3)

with coefficient functions p, q, r and c, and right hand sides f, gD and gN ,

where Ω is a bounded, connected, polygonal, open region in R

2 with boundary

∂Ω = ∂ΩD ∪ ∂ΩN , ∂ΩD ∩ ∂ΩN = ∅. Differentiation with respect to s is with

respect to a counterclockwise parametrization of the boundary (x(s), y(s)) with

k(dx/ds dy/ds)k = 1. Equation 2 represents Dirichlet boundary conditions, and

Equation 3 represents natural boundary conditions if c = 0 and mixed boundary

conditions otherwise. We assume the data in Equations 1-3 satisfy the usual

ellipticity and regularity assumptions. Some of the test problems extend this

to a coupled system of two equations, and the inclusion of first order derivative

and mixed derivative terms. Many of the test problems use the special case of

Poisson’s Equation

−∇2u := −

∂

2u

∂x2

−

∂

2u

∂y2

= f(x, y)

and Laplace’s Equation, which is Poisson’s equation with f=0.

We use the terms singular and singularity rather loosely. We consider a

function to be singular (or to have a singularity) if it or a derivative of any

order is singular. More precisely, we consider it to be singular if there exists a

finite positive m such that the function does not lie in Hm, where Hm is the

usual Sobolev space [5] of functions whose derivatives of order m are square

integrable, and the usual extensions to noninteger m. We refer to the smallest

such m as the Sobolev regularity of the function.

2
Figure 1: The solution of the analytic problem with p = 10.

2 The Set of Test Problems

2.1 Analytic Solution

This is a well behaved problem with a smooth solution that has no trouble spots.

It can be used for seeing how an adaptive algorithm behaves in a context where

adaptivity isn’t really needed.

Equation: Poisson

Domain: unit square

Boundary conditions: Dirichlet

Solution: 24px

p

(1 − x)

pyp

(1 − y)

p

Parameters: p determines the degree of the polynomial solution. It should

be chosen to be large enough that the highest order finite elements to be used

will not give the exact solution.

The solution with p = 10 is shown in Figure 1, both as a color-mapped

image and as a surface in perspective. The other figures that show solutions

also present these two views.

2.2 Reentrant Corner

For elliptic partial differential equations, reentrant corners in the domain cause

a singularity in the solution. In particular, for a corner with an angle ω as shown

in Figure 2, the solution behaves like r

α where r is the distance from the corner

and α = π/ω. The solution is in H1+α− ∀ > 0 [11].

Equation: Laplace

Domain: (−1, 1) × (−1, 1) with a section removed from the clockwise side of

the positive x axis, as shown in Figure 2.

Boundary conditions: Dirichlet

Solution: r

α sin(αθ) where r =

p

x

2 + y2 and θ = tan−1

(y/x)

Parameters: ω determines the angle of the reentrant corner, and consequently α and the strength of the singularity. Varying ω can be used to study the

3
Figure 2: Domain for the reentrant corner problem.

Figure 3: The solution of the reentrant corner problem with ω = π + .01.

effect of the strength of the singularity on adaptive algorithms. Using ω = 3π/2

gives the infamous “L domain” problem used heavily in the adaptive refinement

community. With ω = 2π the domain is a square with a slit. A solution with ω

slightly larger than π is nearly linear.

The solutions for ω = π + 0.01, 5π/4, 3π/2, 7π/4, and 2π are shown in

Figures 3-7.

2.3 Linear Elasticity

Several papers [2, 3, 4, 5] use a problem from linear elasticity as an example.

This is a coupled system of two equations with a mixed derivative in the coupling

4
Figure 4: The solution of the reentrant corner problem with ω = 5π/4.

Figure 5: The solution of the reentrant corner problem with ω = 3π/2.

Figure 6: The solution of the reentrant corner problem with ω = 7π/4.

5
Figure 7: The solution of the reentrant corner problem with ω = 2π.

mode λ Q

1 0.5444837367825 0.5430755788367

2 0.9085291898461 -0.2189232362488

Table 1: Parameters for the linear elasticity problem.

term. In [2] the equations are given as







−E

1 − ν

2

1 − 2ν

∂

2u

∂x2

− E

1 − ν

2

2 − 2ν

∂

2u

∂y2

− E

1 − ν

2

(1 − 2ν)(2 − 2ν)

∂

2v

∂x∂y = Fx

−E

1 − ν

2

2 − 2ν

∂

2v

∂x2

− E

1 − ν

2

1 − 2ν

∂

2v

∂y2

− E

1 − ν

2

(1 − 2ν)(2 − 2ν)

∂

2u

∂x∂y = Fy

(4)

where u and v are the x and y displacements, E is Young’s Modulus, and ν is

Poisson’s ratio.

Two solutions are given in [2] in polar coordinates; a mode 1 solution

(

u =

1

2G

r

λ

[(κ − Q(λ + 1)) cos(λθ) − λ cos((λ − 2)θ)]

v =

1

2G

r

λ

[(κ + Q(λ + 1)) sin(λθ) + λ sin((λ − 2)θ)]

(5)

and a mode 2 solution

(

u =

1

2G

r

λ

[(κ − Q(λ + 1)) sin(λθ) − λ sin((λ − 2)θ)]

v = −

1

2G

r

λ

[(κ + Q(λ + 1)) cos(λθ) + λ cos((λ − 2)θ)]

(6)

where κ = 3 − 4ν, G = E/(2(1 + ν)), and λ and Q are constants. The solutions

have a point singularity at the origin and are in H1+λ− ∀ > 0 [3].

The domain is taken to be a square with a slit (a cracked plate) as in [5].

Some of the other references use an L-shaped domain.

Equation: coupled system of two equations given by Equation 4

Domain: (−1, 1) × (−1, 1) with a slit from (0, 0) to (1, 0)

6
Figure 8: The u component of the mode 1 solution of the linear elasticity problem.

Figure 9: The v component of the mode 1 solution of the linear elasticity problem.

Boundary conditions: Dirichlet

Solution: Two solutions as given in Equations 5 and 6.

Parameters: All the above references use ν = 0.3 and E = 1. The values

for λ and Q differ in the two solutions, and are given in Table 1. See [2]

for the derivation of these constants. With these solutions and parameters,

Fx = Fy = 0.

The solutions are shown in Figures 8-11.

2.4 Peak

This problem has an exponential peak in the interior of the domain. It is based

on Problem 10 in [13].

Equation: Poisson

Domain: Unit square.

Boundary conditions: Dirichlet

7
Figure 10: The u component of the mode 2 solution of the linear elasticity

problem.

Figure 11: The v component of the mode 2 solution of the linear elasticity

problem.

8
Figure 12: The solution of the peak problem with α = 1000, (xc, yc) = (0.5, 0.5).

Figure 13: The solution of the peak problem with α = 100000, (xc, yc) =

(0.51, 0.117).

Solution: e

−α((x−xc)

2+(y−yc)2

)

Parameters: (xc, yc) is the location of the peak, and α determines the

strength of the peak.

The solutions of two instances of this problem are shown in Figures 12 and

13. The first has a mild peak at the vertex (0.5,0.5) with α = 1000. The second

has a sharp peak at (0.51,0.117), which is not a vertex, with α = 100000.

2.5 Battery

This problem comes from [5], where it is attributed to Ivo Babuˇska and Sandia

National Laboratories. It features piecewise constant coefficient functions p

and q and right hand side f, and mixed boundary conditions. The solution

has multiple point singularities in the interior of the domain. The equation

models heat conduction in a battery with nonhomogeneous materials. The

domain is the rectangle shown in Figure 14. The numbered regions indicate the

areas of different material constants, with the constants given in Table 2. The

9
Figure 14: Domain for the battery problem.

k p q f

1 25.0 25.0 0.0

2 7.0 0.8 1.0

3 5.0 0.0001 1.0

4 0.2 0.2 0.0

5 0.05 0.05 0.0

Table 2: Piecewise constant coefficient functions for the battery problem.

location of the line segments that separate the regions are given in Table 3. The

coefficients of the mixed boundary conditions are given in Table 4. The solution

has singularities at the points where three or more materials meet. For any

 > 0 there exists coefficients such that the solution is in H1+. By observing

the rate of convergence with uniform h-refinement and comparing it with the

theoretical a priori error bound, we estimate that  is about 1/2 for the given

set of coefficients. The solution is shown in Figure 15.

Equation: −

∂

∂x 

p

∂u

∂x 

−

∂

∂y 

q

∂u

∂y 

= f with p, q and f given in Table 2

Domain: (0, 8.4) × (0, 24)

Boundary conditions: Mixed, with c and gN given in Table 4

Solution: Unknown.

Parameters: None.

10
x y

0.0 0.0

6.1 0.8

6.5 1.6

8.0 3.6

8.4 18.8

21.2

23.2

24.0

Table 3: The locations of the line segments that separate the regions of the

battery problem. For example, the line y = 24.0 is the top of the region in

Figure 14.

side c gN

left 0.0 0.0

top 1.0 3.0

right 2.0 2.0

bottom 3.0 1.0

Table 4: Boundary condition coefficients for the four sides of the battery problem

domain in Figure 14.

Figure 15: The solution of the battery problem.

11
Figure 16: The solution of the boundary layer problem with  = 10−1

.

Figure 17: The solution of the boundary layer problem with  = 10−3.

2.6 Boundary Layer

This problem comes from [1]. It has an O() boundary layer along the right and

top sides of the domain. It is a convection-diffusion equation with first order

terms.

Equation: −∇2u + 2 ∂u

∂x +

∂u

∂y = f

Domain: (−1, 1) × (−1, 1)

Boundary conditions: Dirichlet

Solution: (1 − e

−(1−x)/)(1 − e−(1−y)/) cos(π(x + y))

Parameters:  determines the strength of the boundary layer.

Figures 16 and 17 show solutions with a mild boundary layer from  = 10−1,

and a steep boundary layer from  = 10−3.

2.7 Boundary Line Singularity

Many papers [5, 6, 10, 12] use a 1D example with a singularity of the form x

α

at the left endpoint of the domain. The solution is in Hα+1/2− ∀ > 0 [6]. This

12
Figure 18: The solution of the boundary line singularity problem with α = 0.6.

can be extended to 2D by simply making the solution be constant in y. On the

unit square, the result is a solution that is singular along the left boundary.

Equation: Poisson

Domain: Unit square

Boundary conditions: Dirichlet

Solution: x

α

Parameters: α ≥ 1/2 determines the strength of the singularity. All of the

cited references use α = 0.6.

The solution with α = 0.6 is shown in Figure 18.

2.8 Oscillatory

This problem is inspired by the wave function that satisfies a Schr ̈odinger equation model of two interacting atoms [8]. It is highly oscillatory near the origin,

with the wavelength decreasing closer to the origin. The number of oscillations,

N, is determined by the parameter α =

1

Nπ . We use a Helmholtz equation for

this problem.

Equation: −∇2u −

1

(α+r)

4 u = f, where r =

p

x

2 + y2

Domain: Unit square

Boundary conditions: Dirichlet

Solution: sin( 1

α+r

)

Parameters: α determines the number of oscillations.

The solutions of a relatively easy problem with α =

1

10π

and a more difficult

problem with α =

1

50π

are shown in Figures 19 and 20. In the surface plot of

Figure 19 we have zoomed in on the origin to show the detail.

2.9 Wave Front

A commonly used example for testing adaptive refinement algorithms is Poisson’s equation with a solution that has a steep wave front in the interior of the

domain. Usually it is a circular wave front given by an arctangent, or sometimes

13
Figure 19: The solution of the oscillatory problem with α =

1

10π

.

Figure 20: The solution of the oscillatory problem with α =

1

50π

.

14
name α xc yc r0

mild 20 -.05 -.05 0.7

steep 1000 -.05 -.05 0.7

asymmetric 1000 1.5 0.25 0.92

well 50 0.5 0.5 0.25

Table 5: Parameters for the wave front problem.

Figure 21: The solution of the mild wave front problem.

a hyperbolic tangent. Parameters determine the steepness and location of the

wave front. With the arctangent wave front, there is a singularity at the center

of the circle. By observing the convergence rate with uniform h-refinement, we

estimate that the solution is in Hm with m ≈ 2 if the center of the circle is in

the closure of the domain, and smooth otherwise.

Equation: Poisson

Domain: Unit square

Boundary conditions: Dirichlet

Solution: tan−1(α(r − r0)) where r =

p

(x − xc)

2 + (y − yc)2

Parameters: (xc, yc) is the center of the circular wave front, r0 is the distance

from the wave front to the center of the circle, and α gives the steepness of the

wave front.

Four example solutions are shown in Figures 21-24, with the parameters

given in Table 5. In the first three we choose the center of the circle to be outside

the domain so that we are examining the performance on the wave front, not

the singularity. These solutions are characterized as a mild wave front, a steep

wave front, and a steep wave front that is not symmetric about the origin. In

the fourth example, the entire circle is inside the domain, resulting in a solution

that is a well with a mild singularity at the bottom.

15
Figure 22: The solution of the steep wave front problem.

Figure 23: The solution of the asymmetric wave front problem.

Figure 24: The solution of the wave front “well” problem.

16
Figure 25: The solution of the interior line singularity problem with α = 2.5

and β = 0.

Figure 26: The solution of the interior line singularity problem with α = 1.1

and β = 0.

2.10 Interior Line Singularity

Houston et al. [6] extend the 1D x

α problem in Section 2.7 to 2D by extending

the 1D domain to (-1,1), defining the solution to be 0 for x < 0, extending the

domain to 2D with y ∈ (−1, 1), and adding cos(πy/2). We extend this further

to allow a sloped line so that the singularity does not necessarily coincide with

element edges. There is no PDE in [6] since they are only evaluating their

estimate of the regularity. We use a Poisson equation with Dirichlet boundary

conditions. The solution is in Hα+1/2− ∀ > 0 [6].

Equation: Poisson

Domain: (−1, 1) × (−1, 1)

Boundary conditions: Dirichlet

Solution: (

cos(πy/2) x ≤ β(y + 1)

cos(πy/2) + (x − β(y + 1))α x > β(y + 1)

Parameters: α determines the strength of the singularity and β determines

the slope of the singularity line.

17
Figure 27: The solution of the interior line singularity problem with α = 1.5

and β = 0.6.

Solutions containing a mild singularity with α = 2.5, β = 0, a stronger

singularity with α = 1.1, β = 0, and a slanted line singularity with α = 1.5, β =

0.6 are shown in Figures 25-27.

2.11 Intersecting Interfaces

This problem comes from a paper by Kellogg [7] in which he studies Poisson

problems with intersecting interfaces. Two interfaces, given by the lines y = 0

and y = tan(φ)x, divide the plane into four regions for a given φ ∈ (0, π/2]. The

PDE coefficients p and q are a piecewise constant function taking the value pi

in the i

th region counterclockwise from the positive x-axis. Let ψ = π/2 − φ.

The solution is given in polar coordinates by

u = r

a1 μ(θ) (7)

where

μ(θ) =







cos((ψ − β1)a1) cos((θ − φ + α1)a1) 0 ≤ θ ≤ φ

cos(α1a1) cos((θ − π + β1)a1) φ <= θ <= π

cos(β1a1) cos(θ − π − α1)a1) π <= θ <= π + φ

cos((φ − α1)a1 cos(θ − φ − π − β1)a1) π + φ <= θ <= 2π

(8)

and where the numbers a1, α1 and β1 satisfy the relations







p1/p2 = − tan((ψ − β1)a1) cot(α1a1)

p2/p3 = − tan(α1a1) cot(β1a1)

p3/p4 = − tan(β1a1) cot((φ − α1)a1)

p4/p3 = − tan((φ − α1)a1) cot((ψ − β1)a1)

0 < a1 < π/ψ

max{0, 2φa1 − π} < 2a1α1 < min{2φa1, π}

max{0, π − 2ψa1} < −2a1β1 < min{π, 2π − 2ψa1}

(9)

18
Figure 28: The solution of the intersecting interfaces problem.

The solution has a discontinuous derivative along the interfaces, and an infinite

derivative at the origin. The solution is in H1+a1− ∀ > 0, so a1 can be chosen

to make the Sobolev regularity arbitrarily close to 1.

Equation: −p∇2u = 0 where p is piecewise constant

Domain: (−1, 1) × (−1, 1)

Boundary conditions: Dirichlet

Solution: given by Equations 7 and 8

Parameters: As in [9], we take θ = π/2, p1 = p3 = R, p2 = p4 = 1,

and a1 = 0.1. The conditions in Equation 9 can then be solved to obtain

R = 161.447 638 797 588 1, α1 = π/4 and β1 = −14.922 565 104 551 52.

The solution for the given parameters is shown in Figure 28.

2.12 Multiple Difficulties

In [12] one of the test cases involves both a singularity due to a reentrant corner

and a sharp gradient in the form of an arctangent wave front. In this problem

we combine four or five difficulties of different strengths into the same problem

by combining some of the features of the other test problems. It contains a

point singularity due to a reentrant corner, a circular wave front (which might

include a singularity at the center of the circle), a sharp peak, and a boundary

layer.

Equation: Poisson

Domain: L-shaped domain (−1, 1) × (−1, 1) \ (0, 1) × (−1, 0)

Boundary conditions: Dirichlet

Solution: r

(π/ω)

sin(θπ/ω) + tan−1(αw(

p

(x − xw)

2 + (y − yw)2 − r0)) +

e

−αp((x−xp)

2+(y−yp)2

) + e−(1+y)/

Parameters: This problem has the same parameters as the problems that

it combines. The boundary layer was placed on y = −1. It could instead be

placed on y = 1, x = −1 or x = 1, or any combination of them.

Figure 29 shows a solution where the wave front intersects the boundary

layer and corner singularity, and the peak is centered on the wave front. The

19
Figure 29: The solution of the problem with multiple difficulties.

following parameters were used. For the reentrant corner, ω = 3π/2. The wave

front is defined by (xw, yw) = (0, −3/4), r0 = 3/4 and αw = 200. The peak is

centered at (xp, yp) = (√

5/4, −1/4) with strength αp = 1000. The boundary

layer is given by  = 1/100.

References

[1] S. Adjerid, M. Aiffa, and J.E. Flaherty, Computational methods for singularly perturbed systems, Singular Perturbation Concepts of Differential

Equations (Providence) (J. Cronin and R.E. O’Malley, eds.), AMS, 1998.

[2] M. Ainsworth and B. Senior, Aspects of an adaptive hp-finite element

method: adaptive strategy, conforming approximation, and efficient solvers,

Comput. Methods Appl. Mech. Engrg. 150 (1997), 65–87.

[3] I. Babuˇska and M. Suri, The h-p version of the finite element method with

quasiuniform meshes, RAIRO Mod ́el. Math. Anal. Num ́er. 21 (1987), 199–

238.

[4] Ivo Babuˇska and Manil Suri, The p- and h-p versions of the finite element

method, an overview, Comput. Methods Appl. Mech. Engrg. 80 (1990),

5–26.

[5] L. Demkowicz, Computing with hp-adaptive finite elements, Volume 1, One

and two dimensional elliptic and Maxwell problems, Chapman & Hall/CRC,

Boca Raton, FL, 2007.

[6] P. Houston, B. Senior, and E. S ̈uli, Sobolev regularity estimation for hpadaptive finite element methods, Numerical Mathematics and Advanced

Appplications (Berlin) (F. Brezzi, A Buffa, S. Corsaro, and A. Murli, eds.),

Springer-Verlag, 2003, pp. 619–644.

20
[7] R.B. Kellogg, On the poisson equation with intersecting interfaces, Appl.

Anal. 4 (1975), 101–129.

[8] W. F. Mitchell and E. Tiesinga, Adaptive grid refinement for a model of two

confined and interacting atoms, Appl. Numer. Math. 52 (2005), 235–250.

[9] P. Morin, R. H. Nochetto, and K. G. Siebert, Data oscillation and convergence of adaptive FEM, SIAM J. Numer. Anal. 38 (2000), no. 2, 466–488.

[10] A. A. Novotny, J. T. Pereira, E. A. Fancello, and C. S. de Barcellos, A

fast hp adaptive finite element mesh design for 2D elliptic boundary value

problems, Comput. Methods Appl. Mech. Engrg. 190 (2000), 133–148.

[11] J. T. Oden and A. Patra, A parallel adaptive strategy for hp finite element

computations, Comput. Methods Appl. Mech. Engrg. 121 (1995), 449–470.

[12] W. Rachowicz, J. T. Oden, and L. Demkowicz, Toward a universal h-p

adaptive finite element strategy, Part 3. Design of h-p meshes, Comput.

Methods Appl. Mech. Engrg. 77 (1989), 181–212.

[13] J.R. Rice, E.N. Houstis, and W.R. Dyksen, A population of linear, second

order, elliptic partial differential equations on rectangular domains, Math.

Comp. 36 (1981), 475–484.

21