BUG for 2D triangle boundary evaluation in dolfinx post5 release

dolfinx
1

Consider the following MWE:

#!/usr/bin/env python3
from dolfinx import mesh, fem, io
import dolfinx.fem.petsc
import ufl
from mpi4py import MPI
import numpy as np

msh = mesh.create_unit_square(MPI.COMM_WORLD, 1, 1, mesh.CellType.triangle)   # quadrilateral, triangle

n0 = ufl.FacetNormal(msh)

tdim = msh.topology.dim
fdim = tdim - 1

facets = mesh.locate_entities(msh, fdim,  lambda x: np.logical_or(np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)), np.logical_or(np.isclose(x[1], 1.0), np.isclose(x[1], 0.0))))

facet_tags = mesh.meshtags(msh, fdim, facets, np.full_like(facets, 1))

ds_ = ufl.Measure("ds", domain=msh, subdomain_data=facet_tags, metadata={"quadrature_degree": 5})

V = fem.functionspace(msh, ("Lagrange", 2, (2,)))
u = fem.Function(V)
varu = ufl.TestFunction(V)

wform = ufl.dot(n0, varu) * ds_(1)
res = fem.form(wform)

r = fem.petsc.assemble_vector(res)

print(r.norm())

which evaluates a constant on the boundary of the unit 2D domain. It should give 1.4142… regardless of triangles or quadrilaterals. It produces this value for both element types in an older dolfinx release (confirmed for Aug 6, 2025 install). However, for the post5 release version, the value deviates to 1.7320… for triangles (while it stays correct for quadrilaterals).

2

This relates to