Nedelec H(curl) interpolation

General
1

Hi,

I was wondering how the interpolation coefficients for H(curl) order 2 functions are calculated in e.g. 2 dimensions?

Lets take e.g. one triangle. As far as I understand for order 1 the 3 coefficients c_i are calculated as follows

\vec{F}(P_i)\cdot \vec{e}_i=\left(\vec{N_i}(P_i)\cdot\vec{e}_i\right)\,c_i

where the P_i are the midpoints of the edges, the \vec{e}_i are the edge vectors, the \vec{N}_i are the base functions and \vec{F} is the given vector field.

For order 2 the vector field \vec{F} is calculated at 9 positions during interpolation (see below). How are the 8 coefficients determined finally?
Is there a way to get these positions from the dofmap?

import numpy as np
import dolfinx as dox
import ufl
from mpi4py import MPI

# points and elements to dolfin mesh
def pym2fm(p,t):

  gdim = 2
  shape = "triangle"
  degree = 1

  cell = ufl.Cell(shape, geometric_dimension=gdim)
  domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))

  cells = np.array(t, dtype=np.int32)
  mesh = dox.mesh.create_mesh(MPI.COMM_WORLD, cells, p, domain)

  return mesh,mesh.geometry.x[:,[0,1]],mesh.geometry.dofmap.array.reshape(-1,3) 

gxy=0
def VectorField(x):
  global gxy  
  gxy = x.T
  values = np.empty( (2, x.shape[1]) )
  values[0] = x[0]
  values[1] = 3*x[1]
  return values

# Definition mesh
pois = [[0,0],[1,0.2],[0.5,1]]
tris = [[0,1,2]]
triangle,pt,tt=pym2fm(pois,tris)

# space and interpolation 
order=2
VV = dox.fem.FunctionSpace(triangle, ("Nedelec 1st kind H(curl)", order))
uF = dox.fem.Function(VV)
uF.interpolate( VectorField )

print(gxy)
print(uF.x.array)
2

For Nedelec elements, the functionals are facet integrals, see for instance:
https://defelement.com/elements/examples/triangle-N1curl-2.html
These integrals are approximated with the relevant quadrature.

To get the interpolation points on the reference element, you can call

print(VV.element.interpolation_points())

which returns

[[0.78867513 0.21132487]
 [0.21132487 0.78867513]
 [0.         0.21132487]
 [0.         0.78867513]
 [0.21132487 0.        ]
 [0.78867513 0.        ]
 [0.16666667 0.16666667]
 [0.16666667 0.66666667]
 [0.66666667 0.16666667]]

THe basis functions are evaluated at these points, then the evaluated functions are combined as defined by:

VV.element.basix_element.interpolation_matrix
3

Thanks dokken,

just to see if I got that right for first order.

For first order Nedelec the c_i are actually line integrals along the edges of \vec{F}\cdot\vec{t}_i, as given in the defelement link for first order. However the line integrals are calculated using quadrature and the number of points used for this quadrature is only one, the midpoint. So if I would like to increase the accuracy of this quadrature integration more interpolation points could be used and the c_i would differ. Not sure if that is technically possible?

The interpolation is finally done with the interpolation_matrix (which is probably dependent on the quadrature points).
Switching to order=1 in my code above and using VV.element.basix_element.interpolation_matrix

FF = VectorField(gxy.T)
cci = np.dot(VV.element.basix_element.interpolation_matrix,FF.flatten())
print( "Matrix*vector :",cci )
print( "Interpolation values :", uF.x.array )

results in
cci = [1.05, 1.5 , 0.5 ])

but the interpolation coefficients are
[1.065, 1.625, 0.56] ?

4

This is correct. As one knows the order of the polynomials along the edges, you do not need to increase the order, as the quadrature rule used is accurate up to the degree of the polynomial.

You are missing the pull-back from the physical space (where the function VectorField is defined) and the reference element. This happens in: https://github.com/FEniCS/dolfinx/blob/main/cpp/dolfinx/fem/interpolate.h#L974
i.e. your basis function is dependent of the jacobian of your cell.
i.e. changing it to:

pois = [[0, 0], [1, 0.0], [0.0, 1]]

you get the same result.
The mapping for N1curl is covariant Piola (see Mapping at: DefElement: Nédélec (first kind)) and “Mapping of finite elements” at: DefElement: How to define a finite element

5

thanks

This is correct. As one knows the order of the polynomials along the edges, you do not need to increase the order, as the quadrature rule used is accurate up to the degree of the polynomial.

Right, I had a general non linear vector field \vec{F} in mind. Is there an option to set the number of integration points? Thanks for your patient.

6

Im not sure why you want to change it.

It is hardcoded in Basix:

and

However, I guess you could change this by making a custom element, as explained in:
https://docs.fenicsproject.org/basix/v0.6.0/python/demo/demo_custom_element.py.html#degree-1-ravairt-thomas-element
and
https://docs.fenicsproject.org/basix/v0.6.0/python/demo/demo_custom_element_conforming_cr.py.html

You can use them in DOLFINx as shown in
https://docs.fenicsproject.org/dolfinx/v0.6.0/python/demos/demo_tnt-elements.html