# Nedelec DOFs give suboptimal interpolation operator

Issue #25 new
Jan Blechta created an issue

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.