Support ufl expressions and submesh quantities in interpolate_to_surface_submesh

[jorgensd]

Following the updates in DOLFINx (FEniCS/dolfinx#4140) we can now evaluate expressions (not backwards compatible) in interpolate_to_surface_submesh.
Additional support for sub-meshes through entity maps.

from mpi4py import MPI
import dolfinx
import ufl
import numpy as np
import scifem


def cell_locator(x):
    return x[0] > 0.1 - 1e-10


mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [[-0.2, -0.2], [1.1, 1.1]], [11, 13])

degree = 3
Vh = dolfinx.fem.functionspace(mesh, ("Lagrange", degree, (mesh.geometry.dim,)))
uh = dolfinx.fem.Function(Vh)
uh.interpolate(
    lambda x: (1.3*(x[0]-0.823)**2-(x[1]-0.9012)**2, (x[0]-0.475)**3 - 0.2*(x[1]-0.852)**2),
    cells0=dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, cell_locator),
)

tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(fdim, tdim)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)

nh = ufl.FacetNormal(mesh)


# Define submesh of the facets you would like to compute the maximum over
facet_submesh, entity_map = dolfinx.mesh.create_submesh(mesh, fdim, exterior_facets)[0:2]
num_submesh_cells = facet_submesh.topology.index_map(fdim).size_local
parent_facets = entity_map.sub_topology_to_topology(
    np.arange(num_submesh_cells, dtype=np.int64), inverse=False
)
exterior_facet_as_integration_entities = dolfinx.fem.compute_integration_domains(
    dolfinx.fem.IntegralType.exterior_facet, mesh.topology, parent_facets
).reshape(-1, 2)


expr = ufl.as_vector((uh[1] * nh[0], uh[0] * nh[1]))

Qh = dolfinx.fem.functionspace(facet_submesh, ("DG", 2, (mesh.geometry.dim,)))
facet_expr = dolfinx.fem.Function(Qh)
scifem.interpolation.interpolate_to_surface_submesh(expr, facet_expr, np.arange(len(parent_facets), dtype=np.int32), exterior_facet_as_integration_entities)

with dolfinx.io.VTXWriter(facet_submesh.comm, "facet_expr.bp", [facet_expr]) as writer:
    writer.write(0.0)

with dolfinx.io.VTXWriter(facet_submesh.comm, "uh.bp", [uh]) as writer:
    writer.write(0.0)


# Use scifem to compute maximum


max_expression = ufl.dot(facet_expr, facet_expr)

value, point_of_extrema = scifem.compute_extrema(max_expression, kind=max)
print(value, point_of_extrema)

yielding

1.4704992620052726 [ 0.823 -0.2  ]

[finsberg]

LGTM