Nedelec DOFs give suboptimal interpolation operator

[Jan Blechta]

Interpolation operator of Nedelec degree 2 gives interpolation operator which does not preserve second order convergence in H(curl) norm. Tentative fix in 94b3af8.

If using the resulting interpolation operator to interpolate Dirichlet datum, it screws up convergence rate. Not so minimal working example:

from dolfin import *
import ufl

import matplotlib.pyplot as plt
import numpy as np

#import dijitso
#dijitso.set_log_level("debug")

#parameters['form_compiler']['representation'] = 'tsfc'
fcp = parameters['form_compiler']
fcp.add('enable_preintegration', False)  # Workaround FFC issue #173

parameters['std_out_all_processes'] = False


class PlaneWave2D(UserExpression):
    def eval(self, values, x):
        u = cos(-1.5*x[0] + x[1])
        values[0] = u
        values[1] = 1.5*u
    def value_shape(self):
        return (2,)


class PlaneWave3D(UserExpression):
    def eval(self, values, x):
        u = cos(-1.5*x[0] + x[1])
        values[0] = u
        values[1] = 1.5*u
        values[2] = 0.5*u
    def value_shape(self):
        return (3,)


def solve_cavity(N, element):
    # Prepare mesh and function space
    dim = element.cell().topological_dimension()
    mesh = {2: UnitSquareMesh, 3: UnitCubeMesh}[dim](*dim*(N,))
    V = FunctionSpace(mesh, element)

    # Prepare exaxt plane wave solution and its wave number
    SolutionClass = {2: PlaneWave2D, 3: PlaneWave3D}[dim]
    u0 = SolutionClass(domain=mesh, element=V.ufl_element())
    omega2 = 1.5**2 + 1.0**2

    # Prepare variational problem
    u, v = TrialFunction(V), TestFunction(V)
    a = (inner(curl(u), curl(v)) - Constant(omega2)*inner(u, v))*dx
    L = inner(Constant(dim*(0,)), v)*dx

    # Solve
    u = Function(V)
    bc = DirichletBC(V, u0, lambda x, b: b)
    solve(a == L, u, bc)

    return u0, u


def test_convergence(dim, degree_list, N_list):
    results = {}

    # Loop over finite elements
    for family in ['Nedelec 1st kind H(curl)', 'Nedelec 2nd kind H(curl)']:
        for degree in degree_list:
            cell = ufl.cell.simplex(dim)
            element = FiniteElement(family, cell, degree)

            # Solve problem on series of meshes and compute errors
            for N in N_list:
                info("****Element: {}, N: {}".format(element, N))
                u_exact, u = solve_cavity(N, element)
                e_L2 = errornorm(u_exact, u, norm_type='L2', degree_rise=1)
                e_Hcurl = errornorm(u_exact, u, norm_type='Hcurl0', degree_rise=1)

                # Store the result
                results.setdefault(element, []).append(
                    (N, u.function_space().dim(), e_L2, e_Hcurl)
                )

    return results


def plot_results(results):
    for element, data in results.items():
        data = np.array(data)

        # Plot L2 error
        plt.subplot(2, 1, 1)
        plt.xlabel('Number DOFs')
        plt.ylabel('L^2 error')
        plt.loglog(data[:,1], data[:,2], '+-', label=element)

        # Plot Hcurl error
        plt.subplot(2, 1, 2)
        plt.xlabel('Number DOFs')
        plt.ylabel('H(curl) error')
        plt.loglog(data[:,1], data[:,3], '+-', label=element)

        # Report convergence rates
        tdim = element.cell().topological_dimension()
        rate_l2_fit = np.polyfit(np.log10(data[:,1]), tdim*np.log10(data[:,2]), 1)[0]
        rate_l2_last = tdim*np.log10(data[-1,2]/data[-2,2])/np.log10(data[-1,1]/data[-2,1])
        rate_hcurl_fit = np.polyfit(np.log10(data[:,1]), tdim*np.log10(data[:,3]), 1)[0]
        rate_hcurl_last = tdim*np.log10(data[-1,3]/data[-2,3])/np.log10(data[-1,1]/data[-2,1])
        info("Element {}:".format(element))
        info("    L2    last: {}, L2    fit: {}".format(rate_l2_last, rate_l2_fit))
        info("    Hcurl last: {}, Hcurl fit: {}".format(rate_hcurl_last, rate_hcurl_fit))

    plt.legend()
    plt.show()


if __name__ == '__main__':
    dim = 2
    degree_list = [1, 2, 3]
    N_list = [2**i for i in range(2, 5)]
    results = test_convergence(dim, degree_list, N_list)
    plot_results(results)

    dim = 3
    degree_list = [1, 2, 3]
    N_list = [2**i for i in range(0, 4)]
    results = test_convergence(dim, degree_list, N_list)
    plot_results(results)

indicates that N1curl2 on tetrahedron is suboptimal and N1curl3 possibly too. Nedelecs of second kind look fine and Nedelecs on triangles too. Similar issue might apply to Raviart-Thomas and/or Brezzi-Douglas-Marini.

Element <N1curl1 on a tetrahedron>:
    L2    last: -1.0769174064599485, L2    fit: -1.1776251217910474
    Hcurl last: -1.0846220310481236, Hcurl fit: -1.1550957717074564
Element <N1curl2 on a tetrahedron>:
    L2    last: -2.0985796089600335, L2    fit: -2.2065164649533564
    Hcurl last: -1.6697668811148338, Hcurl fit: -1.8039142940261197
Element <N1curl3 on a tetrahedron>:
    L2    last: -3.3103730362301333, L2    fit: -3.4812745262580163
    Hcurl last: -2.779362850067247, Hcurl fit: -3.0317195930680985
Element <N2curl1 on a tetrahedron>:
    L2    last: -2.144749286717133, L2    fit: -2.297116180917629
    Hcurl last: -1.0789328691971223, Hcurl fit: -1.1339023221546727
Element <N2curl2 on a tetrahedron>:
    L2    last: -3.178139314923811, L2    fit: -3.3419668206713973
    Hcurl last: -2.111184218331062, Hcurl fit: -2.246378747273589
Element <N2curl3 on a tetrahedron>:
    L2    last: -4.196824248955414, L2    fit: -4.4164526146157606
    Hcurl last: -3.1180590790915796, Hcurl fit: -3.2457133331800105

Note of decreased rate in N1curl2 and N1curl3.