[BUG]: Interpolation on manifold broken

[jorgensd]

Summarize the issue

Reported at:

How to reproduce the bug

Run following test on a one triangle mesh

Minimal Example (Python)

from mpi4py import MPI

import numpy as np
import basix.ufl
import ufl
import dolfinx as dfx
degree = 2
vertices = np.array(
    [(0.0, 0.0, 1.0), (1.0, 1.0, 1.0), (1.0, 0.0, 0.0)],
    dtype=dfx.default_real_type,
)
cells = [(0, 1, 2)]
domain = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(3,), dtype=dfx.default_real_type))
msh = dfx.mesh.create_mesh(MPI.COMM_WORLD, cells, vertices, domain)


N1E = basix.ufl.element(basix.ElementFamily.N1E,basix.CellType.triangle,degree = degree)
BDM = basix.ufl.element(basix.ElementFamily.BDM,basix.CellType.triangle, degree = degree) 
P = basix.ufl.element(basix.ElementFamily.P,basix.CellType.triangle, shape = (3,), degree = degree, discontinuous=True)

N1E = dfx.fem.functionspace(msh, N1E)
BDM = dfx.fem.functionspace(msh, BDM)
P = dfx.fem.functionspace(msh, P)

n1e = dfx.fem.Function(N1E)
bdm = dfx.fem.Function(BDM)
p = dfx.fem.Function(P)

def f(x):
    return 0*x[0], 0*x[1], 1*x[2]



n1e.interpolate(f)
bdm.interpolate(f)
p.interpolate(f)

def error(ph):
    x_ = ufl.SpatialCoordinate(msh)
    f_ex = ufl.as_vector(f(x_))
    diff = ufl.inner(ph - f_ex, ph-f_ex) * ufl.dx
    compiled = dfx.fem.form(diff)
    local_err = dfx.fem.assemble_scalar(compiled)
    return np.sqrt(MPI.COMM_WORLD.allreduce(local_err, op=MPI.SUM))

print(f"{error(p)=}, {error(bdm)=}, {error(n1e)=}")

p_n1e = dfx.fem.Function(P)
p_n1e.interpolate(n1e)
p_bdm = dfx.fem.Function(P)
p_bdm.interpolate(bdm)
print(f"{error(p_n1e)=}, {error(p_bdm)=}")

Output (Python)

error(p)=9.124861240584765e-17, error(bdm)=0.7891307953048737, error(n1e)=0.5645633547617197
error(p_n1e)=0.5645633547617197, error(p_bdm)=0.7891307953048737

Version

main branch

DOLFINx git commit

e64178c

Installation

Docker, ghcr.io/fenics/dolfinx/lab:nightly

Additional information

No response

[jorgensd]

Works on cartesian-aligned manifold:

vertices = np.array(
    [(1.0, 0.0, 0.0), (1.0, 1.0, 0.0), (0.0, 0.0, 0.0)],
    dtype=dfx.default_real_type,
)