[BUG]: interpolation on CR space

[cmaurini]

Summarize the issue

I try to interpolate the jump of a DG function on a CR space and I get unexpected results, see the MWE below.

I use the fact the CR function has dofs on facet mid-points, corresponding to integration points on dS with integration degree 1 to obtain the correct reference result.

Am I doing something wrong or there is a bug in the interpolation function?

How to reproduce the bug

See the script below. I use version 0.8.0.

Minimal Example (Python)

# Description: This script shows a possible bug in the interpolation of the jump of a function on a CR space.

from dolfinx import fem, mesh
import ufl
import basix
import numpy as np
from mpi4py import MPI

comm = MPI.COMM_WORLD
msh = mesh.create_unit_square(
    comm,
    2,
    2,
)
CR_el = basix.ufl.element("CR", msh.basix_cell(), 1, discontinuous=False)
DG_el = basix.ufl.element("CR", msh.basix_cell(), 1, discontinuous=True)

# Take a linear function on DG
V_u = fem.functionspace(msh, DG_el)
u = fem.Function(V_u)
u.interpolate(lambda x: x[0])

# Interpolate on CR the jump (CR has DOFs on the facet mid-points)
# The jump should be zero for the given linear function

V_j = fem.functionspace(msh, CR_el)
j = fem.Function(V_j)
j.interpolate(fem.Expression(ufl.jump(u), V_j.element.interpolation_points()))
# The results is not zero
print(f"Interpolation gives result: {j.vector.array}")

# I use the trick below to interpolate the jump on the CR space
# This gives me what I expect

dS = ufl.Measure("dS", domain=msh, metadata={"quadrature_degree": 1})
jump_1 = fem.assemble_vector(
    fem.form(
        ufl.jump(u)
        * ufl.avg(ufl.TestFunction(V_j))
        / ufl.avg(ufl.FacetArea(msh))
        * dS
    )
)
print(f"Expected result: {jump_1.array}")

Output (Python)

Interpolation gives result: 

[ 2.5e-001  5.0e-001  2.5e-001  7.5e-001 -5.0e-001 -5.0e-001 -5.0e-001
  7.5e-001  5.0e-001  2.5e-001 -5.0e-001  2.5e-001  2.5e-001  2.5e-001
  2.5e-001 -2.0e-323]

Expected result: 

[ 0.00000000e+00 -1.11022302e-16  0.00000000e+00  0.00000000e+00
  0.00000000e+00 -2.77555756e-17  0.00000000e+00  0.00000000e+00
 -1.11022302e-16  0.00000000e+00  0.00000000e+00  0.00000000e+00
  0.00000000e+00 -2.77555756e-17  0.00000000e+00  0.00000000e+00]

Version

0.8.0

DOLFINx git commit

No response

Installation

No response

Additional information

No response

[jorgensd]

I am suprised that the following: j.interpolate(fem.Expression(ufl.jump(u), V_j.element.interpolation_points())) doesn’t throw a compilation error, as the jump operator is not defined at the general interpolation points of any element.

As far as I can recall, there is no specific code path in Expression that would be able to take in two cells (to define a jump). There should be some possibilities to extend this in a similar fashion as what I did with: https://github.com/FEniCS/dolfinx/pull/3062

[cmaurini]

Yes I was surprised too to see that line compiling.

Should it throw an error? What bother me and the reason why I posted is that it gives a wrong result without throwing any error/warning.

[jorgensd]

I think it should. Doing interpolation with jump operators are a bit strange in general, as interpolation is usually done from the point of view of one cell.

I am not sure what one would call that operator, and what it should do on exterior facets.