Mismatch between submesh `dx` and parent `ds` integrals with `entity_maps` in DOLFINx

dolfinx
1

Hello,

I have noticed unexpected behavior when integrating a UFL expression involving spatial derivatives over a parent mesh boundary using entity_maps. The integrand is defined on a 1D boundary submesh of a 2D parent mesh. Integrating with a dx measure defined on the submesh gives the correct result, but integrating with a ds measure defined on the parent mesh using entity_maps yields an incorrect value.

Minimal reproducible example:

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

comm = MPI.COMM_WORLD

msh2d = mesh.create_unit_square(
    comm, 8, 8, mesh.CellType.quadrilateral)

fdim = 1
top_facets = mesh.locate_entities_boundary(msh2d, fdim, lambda x: np.isclose(x[1], 1.0))
mt = mesh.meshtags(msh2d, fdim, top_facets, 1.0)

msh1d, entity_map = mesh.create_submesh(msh2d, fdim, top_facets)[0:2]

entity_maps = [entity_map]

dx1 = ufl.Measure("dx", domain=msh1d)
ds2 = ufl.Measure("ds", domain=msh2d, subdomain_id=1, subdomain_data=mt)

V = fem.functionspace(msh1d, ("Lagrange", 1))
u = fem.Function(V)
u.interpolate(lambda x: x[0])

w = ufl.Dx(u, 0) + 1.0

v1 = comm.allreduce(assemble_scalar(fem.form(w * dx1)), op=MPI.SUM)
v2 = comm.allreduce(
    assemble_scalar(fem.form(w * ds2, entity_maps=entity_maps)), op=MPI.SUM
)

if comm.rank == 0:
    print(f"v1:   {v1}")
    print(f"v2:   {v2}")
    print(f"diff: {abs(v1 - v2)}")

Output:

v1:   2.0
v2:   0.9999999999999999
diff: 1.0

Both v1 and v2 should equal 2.0. The integrand w = Dx(u, 0) + 1.0 with u(x) = x[0] evaluates to 2 everywhere on the top edge (length 1), so the integral should be 2.0. Here, the constant term integrates correctly, but the derivative term Dx(u, 0) appears to contribute 0.0 when using ds with entity_maps.

Environment:

  • DOLFINx 0.10.0.post2 (docker v0.10.0-r1)
  • Python 3.12.3
  • Ubuntu 24.04.3 LTS

Is this a known issue or am I making a mistake with entity_maps / measures? Thank you.

2

Hm, something is definitely wrong here.
Especially because if we change the derivative to be in “y”- direction, the result changes:

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

comm = MPI.COMM_WORLD

msh2d = mesh.create_unit_square(comm, 1, 1, mesh.CellType.quadrilateral)

fdim = 1
top_facets = mesh.locate_entities_boundary(msh2d, fdim, lambda x: np.isclose(x[1], 1.0))
mt = mesh.meshtags(msh2d, fdim, top_facets, 1.0)

msh1d, entity_map = mesh.create_submesh(msh2d, fdim, top_facets)[0:2]

entity_maps = [entity_map]

dx1 = ufl.Measure("dx", domain=msh1d)
ds2 = ufl.Measure("ds", domain=msh2d, subdomain_id=1, subdomain_data=mt)

V = fem.functionspace(msh1d, ("Lagrange", 1))

u = fem.Function(V)
u.interpolate(lambda x: x[0])
w = ufl.Dx(u, 1)
f1 = fem.form(w * dx1)
f2 = fem.form(w * ds2, entity_maps=entity_maps)
v1 = comm.allreduce(assemble_scalar(f1), op=MPI.SUM)
v2 = comm.allreduce(assemble_scalar(f2), op=MPI.SUM)


if comm.rank == 0:
    print(f"v1:   {v1}")
    print(f"v2:   {v2}")
    print(f"diff: {abs(v1 - v2)}")

Issue made at:

1 Like