This section contains the description of selected benchmarks. Contrary to regular examples, benchmarks typically do not have a significant physical or engineering motivation but they come with a known exact solution and thus they are a great resource for comparisons of various methods and adaptivity algorithms.
Git reference: Benchmark smooth-iso.
We show that it is a very bad idea to approximate smooth solutions using low-order elements.
Equation solved: Poisson equation
(1)
Domain of interest: Square
.
Right-hand side:
(2)
Boundary conditions: Zero Dirichlet.
Exact solution:
(3)
Code for the exact solution and the weak forms:
// Exact solution.
static double fn(double x, double y)
{
return sin(x)*sin(y);
}
static double fndd(double x, double y, double& dx, double& dy)
{
dx = cos(x)*sin(y);
dy = sin(x)*cos(y);
return fn(x, y);
}
// Boundary condition types.
BCType bc_types(int marker)
{
return BC_ESSENTIAL;
}
// Essential (Dirichlet) boundary conditions.
scalar essential_bc_values(int ess_bdy_marker, double x, double y)
{
return 0;
}
// Weak forms.
template<typename Real, typename Scalar>
Scalar bilinear_form(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *u, Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
return int_grad_u_grad_v<Real, Scalar>(n, wt, u, v);
}
template<typename Real>
Real rhs(Real x, Real y)
{
return 2*sin(x)*sin(y);
}
template<typename Real, typename Scalar>
Scalar linear_form(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
return int_F_v<Real, Scalar>(n, wt, rhs, v, e);
}
Solution:

Below we show meshes obtained using various types of adaptivity. Note the tremendous differences in their performance. The meshes do not correspond to the same level of accuracy since the low-order methods could not achieve the same error as hp-FEM. Therefore, compare not only the number of DOF but also the error level. Convergence graphs for all cases are shown at the end of this section.
Final mesh (h-FEM, p=1): 27469 DOF, error 0.39173795799476 %

Final mesh (h-FEM, p=2): 39185 DOF, error 0.0022127484879974 %

Final mesh (hp-FEM): 49 DOF, error 4.2775412425017e-05 %

DOF convergence graphs:

CPU time convergence graphs:

Git reference: Benchmark smooth-aniso-x.
We show that one should use (spatially as well as polynomially) anisotropic refinements for solutions containing anisotropy.
Equation solved: Poisson equation
(4)
Domain of interest: Square
.
Right-hand side:
(5)
Boundary conditions: Zero Dirichlet on the left and right edges, zero Neumann on the rest of the boundary.
Exact solution:
(6)
Solution:

Below we show meshes obtained using various types of adaptivity. Note the tremendous differences in their performance. The meshes do not correspond to the same level of accuracy since the low-order methods could not achieve the same error as hp-FEM. Therefore, compare not only the number of DOF but also the error level. Convergence graphs for all cases are shown at the end of this section.
Final mesh (h-FEM, p=1, isotropic refinements): 41033 DOF, error 0.22875054074711 %

Final mesh (h-FEM, p=1, anisotropic refinements): 39594 DOF, error 0.0039444224349215 %

Final mesh (h-FEM, p=2, isotropic refinements): 54627 DOF, error 0.0017755772528929 %

Final mesh (h-FEM, p=2, anisotropic refinements): 3141 DOF, error 9.3084842840514e-05 %

Final mesh (hp-FEM, isotropic refinements): 63 DOF, error = 3.6797337289125e-05 %

Final mesh (hp-FEM, anisotropic refinements): 14 DOF, error 3.6797337292196e-05 %, The color pattern means that the polynomial degrees are one and eight in the vertical and horizontal directions, respectively.

DOF convergence graphs:

CPU time convergence graphs:

Git reference: Benchmark smooth-aniso-y.
This example is very similar to the previous one, except now the solution is constant in the x-direction. It is good to have both to be able to check that anisotropic refinements work correctly.
Git reference: Benchmark lshape.
This is a standard adaptivity benchmark whose exact solution is smooth but contains singular gradient in a re-entrant corner.
Equation solved: Laplace equation
(7)
Domain of interest:

Exact solution:
(8)
where
and
.
Code for the exact solution, bundary conditions, and weak forms:
// Exact solution.
static double fn(double x, double y)
{
double r = sqrt(x*x + y*y);
double a = atan2(x, y);
return pow(r, 2.0/3.0) * sin(2.0*a/3.0 + M_PI/3);
}
static double fndd(double x, double y, double& dx, double& dy)
{
double t1 = 2.0/3.0*atan2(x, y) + M_PI/3;
double t2 = pow(x*x + y*y, 1.0/3.0);
double t3 = x*x * ((y*y)/(x*x) + 1);
dx = 2.0/3.0*x*sin(t1)/(t2*t2) + 2.0/3.0*y*t2*cos(t1)/t3;
dy = 2.0/3.0*y*sin(t1)/(t2*t2) - 2.0/3.0*x*t2*cos(t1)/t3;
return fn(x, y);
}
// Boundary condition types.
BCType bc_types(int marker)
{
return BC_ESSENTIAL;
}
// Essential (Dirichlet) boundary condition values.
scalar essential_bc_values(int ess_bdy_marker, double x, double y)
{
return fn(x, y);
}
// Bilinear form corresponding to the Laplace equation.
template<typename Real, typename Scalar>
Scalar bilinear_form(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *u, Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
return int_grad_u_grad_v<Real, Scalar>(n, wt, u, v);
}
Solution:

Final mesh (h-FEM with linear elements):

Final mesh (h-FEM with quadratic elements):

Final mesh (hp-FEM):

DOF convergence graphs:

CPU time convergence graphs:

Git reference: Benchmark layer.
This example has a smooth solution that exhibits a steep internal layer inside the domain.
Equation solved: Poisson equation
(9)
Domain of interest: Unit square
.
Right-hand side:
(10)
Exact solution:
(11)
where
is a parameter (slope of the layer). With larger
, this problem
becomes difficult for adaptive algorithms, and at the same time the advantage of
adaptive
-FEM over adaptive low-order FEM becomes more significant. We will
use
in the following.
Code for the exact solution and the weak forms:
// Exact solution.
static double fn(double x, double y)
{
return atan(SLOPE * (sqrt(sqr(x-1.25) + sqr(y+0.25)) - M_PI/3));
}
static double fndd(double x, double y, double& dx, double& dy)
{
double t = sqrt(sqr(x-1.25) + sqr(y+0.25));
double u = t * (sqr(SLOPE) * sqr(t - M_PI/3) + 1);
dx = SLOPE * (x-1.25) / u;
dy = SLOPE * (y+0.25) / u;
return fn(x, y);
}
// Boundary condition types.
BCType bc_types(int marker)
{
return BC_ESSENTIAL;
}
// Essential (Dirichlet) boundary condition values.
scalar essential_bc_values(int ess_bdy_marker, double x, double y)
{
return fn(x, y);
}
// Bilinear form for the Poisson equation.
template<typename Real, typename Scalar>
Scalar bilinear_form(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *u, Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
return int_grad_u_grad_v<Real, Scalar>(n, wt, u, v);
}
template<typename Real>
Real rhs(Real x, Real y)
{
Real t2 = sqr(y + 0.25) + sqr(x - 1.25);
Real t = sqrt(t2);
Real u = (sqr(M_PI - 3.0*t)*sqr(SLOPE) + 9.0);
return 27.0/2.0 * sqr(2.0*y + 0.5) * (M_PI - 3.0*t) * pow(SLOPE,3.0) / (sqr(u) * t2) +
27.0/2.0 * sqr(2.0*x - 2.5) * (M_PI - 3.0*t) * pow(SLOPE,3.0) / (sqr(u) * t2) -
9.0/4.0 * sqr(2.0*y + 0.5) * SLOPE / (u * pow(t,3.0)) -
9.0/4.0 * sqr(2.0*x - 2.5) * SLOPE / (u * pow(t,3.0)) +
18.0 * SLOPE / (u * t);
}
template<typename Real, typename Scalar>
Scalar linear_form(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
return -int_F_v<Real, Scalar>(n, wt, rhs, v, e);
}
Solution:

Final mesh (h-FEM with linear elements):

Final mesh (h-FEM with quadratic elements):

Final mesh (hp-FEM):

DOF convergence graphs:

CPU time convergence graphs:

Git reference: Benchmark layer-2.
This example is a singularly perturbed problem with known exact solution that exhibits a thin boundary layer, that the reader can use to perform various experiments with adaptivity for problems with boundary layers. The sample numerical results presented below imply that:
Equation solved: Poisson equation
(12)
Domain of interest: Square
.
Exact solution:
where
is the exact solution of the 1D singularly-perturbed problem
in
with zero Dirichlet boundary conditions. This solution has the form
Right-hand side: Calculated by inserting the exact solution into the equation. Here is the code snippet with both the exact solution and the right-hand side:
// Solution to the 1D problem -u'' + K*K*u = K*K in (-1,1) with zero Dirichlet BC.
double uhat(double x) {
return 1. - (exp(K*x) + exp(-K*x)) / (exp(K) + exp(-K));
}
double duhat_dx(double x) {
return -K * (exp(K*x) - exp(-K*x)) / (exp(K) + exp(-K));
}
double dduhat_dxx(double x) {
return -K*K * (exp(K*x) + exp(-K*x)) / (exp(K) + exp(-K));
}
// Exact solution u(x,y) to the 2D problem is defined as the
// Cartesian product of the 1D solutions.
static double sol_exact(double x, double y, double& dx, double& dy)
{
dx = duhat_dx(x) * uhat(y);
dy = uhat(x) * duhat_dx(y);
return uhat(x) * uhat(y);
}
// Right-hand side.
double rhs(double x, double y) {
return -(dduhat_dxx(x)*uhat(y) + uhat(x)*dduhat_dxx(y)) + K*K*uhat(x)*uhat(y);
}
The weak forms are very simple and they are defined as follows. The only thing worth mentioning here is that we integrate the non-polynomial right-hand side with a very hign order for accuracy:
// Weak forms.
template<typename Real, typename Scalar>
Scalar bilinear_form(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *u, Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
return int_grad_u_grad_v<Real, Scalar>(n, wt, u, v) + K*K * int_u_v<Real, Scalar>(n, wt, u, v);
}
template<typename Real, typename Scalar>
Scalar linear_form(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
return int_F_v<Real, Scalar>(n, wt, rhs, v, e);;
}
// Integration order for linear_form_0.
Ord linear_form_ord(int n, double *wt, Func<Ord> *u_ext[], Func<Ord> *v, Geom<Ord> *e, ExtData<Ord> *ext)
{
return 24;
}
The numerical results follow:
Solution:

Below we present a series of convergence comparisons. Note that the error plotted is the true approximate error calculated wrt. the exact solution given above.
Let us first compare the performance of h-FEM (p=1), h-FEM (p=2) and hp-FEM with isotropic refinements:
Final mesh (h-FEM, p=1, isotropic refinements):

Final mesh (h-FEM, p=2, isotropic refinements):

Final mesh (hp-FEM, isotropic refinements):

DOF convergence graphs:

CPU convergence graphs:

Next we compare the performance of h-FEM (p=1), h-FEM (p=2) and hp-FEM with anisotropic refinements:
Final mesh (h-FEM, p=1, anisotropic refinements):

Final mesh (h-FEM, p=2, anisotropic refinements):

Final mesh (hp-FEM, anisotropic refinements):

DOF convergence graphs:

CPU convergence graphs:

DOF convergence graphs:

CPU convergence graphs:

DOF convergence graphs:

CPU convergence graphs:

In the hp-FEM one has two kinds of anisotropy – spatial and polynomial. In the following, “iso” means isotropy both in h and p, “aniso h” means anisotropy in h only, and “aniso hp” means anisotropy in both h and p.
DOF convergence graphs (hp-FEM):

CPU convergence graphs (hp-FEM):

The reader can see that enabling polynomially anisotropic refinements in the hp-FEM is equally important as allowing spatially anisotropic ones.
Git reference: Benchmark line-singularity.
The is another example with anisotropic solution that is suitable for testing anisotropic element refinements.
Equation solved: Poisson equation
(13)
Domain of interest: Square
.
Boundary conditions: Zero Neumann on left edge, Dirichlet given by the exact solution on the rest of the boundary.
Exact solution:
where
and
are real constants.
Right-hand side: Obtained by inserting the exact solution into the equation. The corresponding code snippet is shown below:
scalar rhs(scalar x, scalar y)
{
if (x < 0) return fn(x, y)*K*K;
else return fn(x, y)*K*K-ALPHA*(ALPHA-1)*pow(x, ALPHA - 2.) - K*K*pow(x, ALPHA);
}
Solution for
and
:

Final mesh (h-FEM, p=1, anisotropic refinements):

Final mesh (h-FEM, p=2, anisotropic refinements):

Final mesh (hp-FEM, h-anisotropic refinements):

DOF convergence graphs:

CPU convergence graphs:

Final mesh (hp-FEM, isotropic refinements):

Final mesh (hp-FEM, h-anisotropic refinements):

Final mesh (hp-FEM, hp-anisotropic refinements):

The following convergence comparisons still correspond to an older version of Hermes2D when we did not have anisotropic hp-refinements. With those, the convergence of adaptive hp-FEM has improved a lot. These results will be updated soon.
DOF convergence graphs:

CPU convergence graphs:

Git reference: Benchmark kellogg.
The solution to this elliptic problems contains a severe singularity that poses a challenge to adaptive methods.
Equation solved:
where the parameter
is piecewise-constant,
in the first and third quadrants and
in the remaining two quadrants.
Domain of interest: Square
.
Right-hand side:
.
Boundary conditions: Dirichlet given by exact solution.
Exact solution: Quite complicated, see the code below.
// Problem constants.
const double R = 161.4476387975881; // Equation parameter.
const double TAU = 0.1; // Equation parameter.
const double RHO = M_PI/4.; // Equation parameter
const double SIGMA = -14.92256510455152; // Equation parameter
// Exact solution.
static double fn(double x, double y)
{
double theta = atan2(y,x);
if (theta < 0) theta = theta + 2.*M_PI;
double r = sqrt(x*x + y*y);
double mu;
if (theta <= M_PI/2.) {
mu = cos((M_PI/2. - SIGMA)*TAU) * cos((theta - M_PI/2. + RHO)*TAU);
}
else {
if (theta <= M_PI) {
mu = cos(RHO*TAU) * cos((theta - M_PI + SIGMA)*TAU);
}
else {
if (theta <= 3.*M_PI/2.) {
mu = cos(SIGMA*TAU) * cos((theta - M_PI - RHO)*TAU);
}
else {
mu = cos((M_PI/2. - RHO)*TAU) * cos((theta - 3.*M_PI/2. - SIGMA)*TAU);
}
}
}
return pow(r, TAU) * mu;
}
The weak forms are as follows:
// Weak forms
template<typename Real, typename Scalar>
Scalar bilinear_form_I_III(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *u, Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
return R*int_grad_u_grad_v<Real, Scalar>(n, wt, u, v);
}
template<typename Real, typename Scalar>
Scalar bilinear_form_II_IV(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *u, Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
return 1.*int_grad_u_grad_v<Real, Scalar>(n, wt, u, v);
}
Solution:

Final mesh (h-FEM with linear elements):

Final mesh (h-FEM with quadratic elements):

Final mesh (hp-FEM):

DOF convergence graphs:

CPU time convergence graphs:

Git reference: Benchmark bessel.
This example solves time-harmonic Maxwell’s equations in an L-shaped domain and it describes the diffraction of an electromagnetic wave from a re-entrant corner. It comes with an exact solution that contains singularity.
Equation solved: Time-harmonic Maxwell’s equations
(14)
Domain of interest is the square
missing the quarter lying in the
fourth quadrant. It is filled with air:

Boundary conditions: Combined essential and natural, see the main.cpp file.
Exact solution:
(15)
where
is the Bessel function of the first kind,
the polar coordinates and
. In
computer code, this reads:
void exact_sol(double x, double y, scalar& e0, scalar& e1)
{
double t1 = x*x;
double t2 = y*y;
double t4 = sqrt(t1+t2);
double t5 = jv(-1.0/3.0,t4);
double t6 = 1/t4;
double t7 = jv(2.0/3.0,t4);
double t11 = (t5-2.0/3.0*t6*t7)*t6;
double t12 = atan2(y,x);
if (t12 < 0) t12 += 2.0*M_PI;
double t13 = 2.0/3.0*t12;
double t14 = cos(t13);
double t17 = sin(t13);
double t18 = t7*t17;
double t20 = 1/t1;
double t23 = 1/(1.0+t2*t20);
e0 = t11*y*t14-2.0/3.0*t18/x*t23;
e1 = -t11*x*t14-2.0/3.0*t18*y*t20*t23;
}
Here jv() is the Bessel function
. For its source code see the
forms.cpp file.
Code for the weak forms:
template<typename Real, typename Scalar>
Scalar bilinear_form(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *u, Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
return 1.0/mu_r * int_curl_e_curl_f<Real, Scalar>(n, wt, u, v) -
sqr(kappa) * int_e_f<Real, Scalar>(n, wt, u, v);
}
template<typename Real, typename Scalar>
Scalar bilinear_form_surf(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *u, Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
cplx ii = cplx(0.0, 1.0);
return ii * (-kappa) * int_e_tau_f_tau<Real, Scalar>(n, wt, u, v, e);
}
scalar linear_form_surf(int n, double *wt, Func<scalar> *u_ext[], Func<double> *v, Geom<double> *e, ExtData<scalar> *ext)
{
scalar result = 0;
for (int i = 0; i < n; i++)
{
double r = sqrt(e->x[i] * e->x[i] + e->y[i] * e->y[i]);
double theta = atan2(e->y[i], e->x[i]);
if (theta < 0) theta += 2.0*M_PI;
double j13 = jv(-1.0/3.0, r), j23 = jv(+2.0/3.0, r);
double cost = cos(theta), sint = sin(theta);
double cos23t = cos(2.0/3.0*theta), sin23t = sin(2.0/3.0*theta);
double Etau = e->tx[i] * (cos23t*sint*j13 - 2.0/(3.0*r)*j23*(cos23t*sint + sin23t*cost)) +
e->ty[i] * (-cos23t*cost*j13 + 2.0/(3.0*r)*j23*(cos23t*cost - sin23t*sint));
result += wt[i] * cplx(cos23t*j23, -Etau) * ((v->val0[i] * e->tx[i] + v->val1[i] * e->ty[i]));
}
return result;
}
// Maximal polynomial order to integrate surface linear form.
Ord linear_form_surf_ord(int n, double *wt, Func<Ord> *u_ext[], Func<Ord> *v, Geom<Ord> *e, ExtData<Ord> *ext)
{ return Ord(v->val[0].get_max_order()); }
Solution:

Final mesh (h-FEM with linear elements):

Note that the polynomial order indicated corresponds to the tangential components of approximation on element interfaces, not to polynomial degrees inside the elements (those are one higher).
Final mesh (h-FEM with quadratic elements):

Final mesh (hp-FEM):

DOF convergence graphs:

CPU time convergence graphs:

Git reference: Benchmark screen.
This example solves time-harmonic Maxwell’s equations. It describes an electromagnetic wave that hits a thin screen under the angle of 45 degrees, causing a singularity at the tip of the screen. The strength of the singularity makes this example rather difficult.
Equation solved: Time-harmonic Maxwell’s equations
(16)
Domain of interest is the square
missing the edge that connects the center with
the midpoint of the left side. It is filled with air:

Boundary conditions: Tangential component of solution taken from known exact solution (essential BC). See the main.cpp file.
Exact solution: This is rather complicated in this case - please look into the corresponding file exact_sol.cpp.
Code for the weak forms:
template<typename Real, typename Scalar>
Scalar bilinear_form(int n, double *wt, Func<Scalar> *u_ext[], Func<Real> *u, Func<Real> *v, Geom<Real> *e, ExtData<Scalar> *ext)
{
return int_curl_e_curl_f<Real, Scalar>(n, wt, u, v) - int_e_f<Real, Scalar>(n, wt, u, v);
}
Solution (real part of
):

Solution (real part of
):

Solution (imaginary part of
):

Solution (imaginary part of
):

Final mesh (h-FEM with linear elements):

Note that the polynomial order indicated corresponds to the tangential components of approximation on element interfaces, not to polynomial degrees inside the elements (those are one higher).
Final mesh (h-FEM with quadratic elements):

Final mesh (hp-FEM):

DOF convergence graphs:

CPU time convergence graphs:
