Double derivative of Facet normal broken for higher order meshes

[jorgensd]

First referenced:
https://fenicsproject.discourse.group/t/ulf-issue-for-higher-order-derivatives-on-quadratic-elements/17075/2

Minimal reproducible example

from ufl import (
    Mesh,
    tetrahedron,
    FacetNormal,
    ds,
)
from ufl.algorithms import compute_form_data
from ufl.finiteelement import FiniteElement
from ufl.pullback import identity_pullback
from ufl.sobolevspace import H1

domain = Mesh(FiniteElement("Lagrange", tetrahedron, 2, (3,), identity_pullback, H1))

# Define the facet normal
n = FacetNormal(domain)

# Define the boundary integral
c = n[0].dx(0).dx(0) * ds


A = compute_form_data(c)

[jorgensd]

It seems like for the expression at hand, one has to apply geometry lowering prior to computing derivatives, i.e.

form = ufl.algorithms.apply_algebra_lowering.apply_algebra_lowering(c)
form = ufl.algorithms.remove_complex_nodes.remove_complex_nodes(form)
form = ufl.algorithms.apply_geometry_lowering.apply_geometry_lowering(form)
form = ufl.algorithms.apply_derivatives.apply_derivatives(form)

Does anyone have an idea of how we could work around this?