Projecting an expression onto a boundary submesh

dolfinx
1

I’m trying to do this projection as part of a post-processing step. I have a Function over the entire volume, and I want to plot an expression containing that Function on a few boundaries.

I think this should illustrate what I’m trying to do:

import numpy as np
from dolfinx import fem, mesh
from mpi4py import MPI
import pyvista4dolfinx as p4d
import ufl

domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, mesh.CellType.triangle)

# Simple demonstration vector function f(x, y) = (x, y)
V = fem.functionspace(domain, ("Lagrange", 1, (2,)))
u = fem.Function(V)
x = ufl.SpatialCoordinate(domain)
u.interpolate(lambda x: (x[0] + 0.1, x[1] + 0.1))

p4d.plot(u, show=True, factor=0.1)

# Create a subdomain over the top boundary
fdim = domain.topology.dim - 1
top_facets = mesh.locate_entities_boundary(
    domain, fdim, lambda x: np.isclose(x[1], 1.0)
)
subdomain, entity_map, _, _ = mesh.create_submesh(domain, fdim, top_facets)

V_boundary = fem.functionspace(subdomain, ("Lagrange", 1, (2,)))
u_sub = fem.Function(V_boundary)

expr = fem.Expression(
    u, V_boundary.element.interpolation_points, entity_maps=[entity_map]
)
u_sub.interpolate(expr)

plt = p4d.plot(u_sub, factor=0.1)
p4d.plot(domain, plotter=plt, show=True)

But it’s not working:

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0

I’ve searched the forum and found a few attempts at doing precisely what I’m trying to do, but they’re all fairly outdated. I also looked at several interpolation and expression tests, nothing got me quite there.

There’s an example with an expression over quadrature points, but the only option seems to be to eval, which is not really helpful for plotting.

The closest I got was by handrolling an L2 projector, but there’s always a single facet with a completely wrong value that messes up the visualization (see issue 4291 in github).

Some help is really appreciated! Thanks!

2

This should fix if you only need this projection once as part of postprocessing:

# Get the number of cells in the submesh
sub_cell_map = subdomain.topology.index_map(subdomain.topology.dim)
num_sub_cells = sub_cell_map.size_local + sub_cell_map.num_ghosts
# For parent -> submesh interpolation
interpolation_data = fem.create_interpolation_data(V_boundary, V, cells=np.arange(num_sub_cells))

# interpolate the function u from the parent mesh to the submesh
u_sub.interpolate_nonmatching(u, cells=np.arange(num_sub_cells), interpolation_data=interpolation_data)

However, this approach is more general for non-matching meshes, while in your case the submesh is coincident with the parent (entity_map provides cell-to-cell or facet-to-facet correspondance between parent and subsmesh).
Another alternative is the one here: Interpolate Expression on submesh with parent mesh function

3

Thanks! The first suggestion works for u, but unfortunately not for ufl.dot(u, u). I really need it to work on an Expression..

The alternative suggestion is what I had tried before (with some updates from test_submesh_codim_one()), but unfortunately I can’t seem to be able to plot the result directly:

import numpy as np
import basix
from dolfinx import fem, mesh
from dolfinx.fem import petsc
from mpi4py import MPI
import pyvista4dolfinx as p4d
import ufl

domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, mesh.CellType.triangle)

# Simple demonstration vector function f(x, y) = (x, y)
V = fem.functionspace(domain, ("Lagrange", 1, (2,)))
u = fem.Function(V)
x = ufl.SpatialCoordinate(domain)
u.interpolate(lambda x: (x[0] + 0.1, x[1] + 0.1))

p4d.plot(u, show=True, factor=0.1)

# Create a subdomain over the top boundary
fdim = domain.topology.dim - 1
top_facets = mesh.locate_entities_boundary(
    domain, fdim, lambda x: np.isclose(x[1], 1.0)
)
subdomain, entity_map, _, _ = mesh.create_submesh(domain, fdim, top_facets)

# Create function space on subdomain
q_el = basix.ufl.quadrature_element(
    subdomain.basix_cell(), (), scheme="default", degree=3
)
Q = fem.functionspace(subdomain, q_el)
q = fem.Function(Q)

# Create expression on boundary
n = ufl.FacetNormal(domain)
expr = fem.Expression(ufl.dot(u, n), Q.element.interpolation_points)

# Get integration entities and evaluate expression
entities = fem.compute_integration_domains(
    fem.IntegralType.exterior_facet, domain.topology, top_facets
)
q.x.array[:] = expr.eval(domain, entities.reshape((-1, 2))).flatten()

plt = p4d.plot(q, factor=0.1)
p4d.plot(domain, plotter=plt, show=True)

Traceback (most recent call last):
  File "x.py", line 43, in <module>
    plt = p4d.plot(q, factor=0.1)
  File "x/.pixi/envs/dev/lib/python3.14/functools.py", line 982, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
           ~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^
  File "x/.pixi/envs/dev/lib/python3.14/site-packages/pyvista4dolfinx/plot.py", line 129, in plot_function_
    return plot_function_scalar(plottable, *args, **kwargs)
  File "x/.pixi/envs/dev/lib/python3.14/site-packages/pyvista4dolfinx/plot.py", line 295, in plot_function_scalar
    u, warp = _upgrade_u_and_warp_to_geometry(u, warp, min_degree=1)
              ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^
  File "x/.pixi/envs/dev/lib/python3.14/site-packages/pyvista4dolfinx/plot.py", line 1145, in _upgrade_u_and_warp_to_geometry
    poldeg = V.element.basix_element.degree
             ^^^^^^^^^^^^^^^^^^^^^^^
  File "x/.pixi/envs/dev/lib/python3.14/site-packages/dolfinx/fem/element.py", line 209, in basix_element
    return self._cpp_object.basix_element
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: No Basix element available. Maybe this is a mixed element?

Tried interpolating into a new CG function over the subdomain, didn’t work:

new_element = basix.ufl.element("DG", subdomain.basix_cell(), 1)
V_new = fem.functionspace(subdomain, new_element)
new_function = fem.Function(V_new)
new_function.interpolate(q)

Traceback (most recent call last):
  File "x.py", line 47, in <module>
    new_function.interpolate(q)
    ~~~~~~~~~~~~~~~~~~~~~~~~^^^
  File "x/.pixi/envs/dev/lib/python3.14/site-packages/dolfinx/fem/function.py", line 537, in interpolate
    _interpolate(u0)
    ~~~~~~~~~~~~^^^^
  File "x/.pixi/envs/dev/lib/python3.14/functools.py", line 982, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
           ~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^
  File "x/.pixi/envs/dev/lib/python3.14/site-packages/dolfinx/fem/function.py", line 523, in _
    self._cpp_object.interpolate(u0._cpp_object, cells0, cells1)
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Missing a Basix element. Cannot check for equivalence
4

Got it working! I just had to project q onto a new function f:

V_boundary = fem.functionspace(subdomain, ("DG", 1))
f = fem.Function(V_boundary)
u = ufl.TrialFunction(V_boundary)
v = ufl.TestFunction(V_boundary)

dx_sub = ufl.Measure("dx", domain=subdomain)
a = ufl.inner(u, v) * dx_sub
L = ufl.inner(q, v) * dx_sub

problem = petsc.LinearProblem(
    a,
    L,
    u=f,
    petsc_options={"ksp_type": "preonly", "pc_type": "lu"},
    petsc_options_prefix="projection",
)
problem.solve()

Thank you for pointing me on the right direction!!

5

Instead of interpolating onto a quadrature space on the submesh, you can use the function: scifem.interpolation.interpolate_to_surface_submesh to interpolate onto any function space on the facet mesh.

This works for general ufl.core.expr.Expr since v0.11.0

6

Thanks, but to be very honest, I just couldn’t figure out how to use interpolate_to_surface_submesh. The API documentation is outdated, and I have no idea what it expects in submesh_facets and integration_entities..

You mean in u_volume? The function signature expects a fem.Function..

BTW I did a benchmark with a 2000x2000 domain (same petsc_options) between the approach in #4291 and the final approach above, and surprisingly interpolating to a quadrature space is visibly faster (~1950 ms vs. ~2050 ms single process). I would expect the opposite, given it takes more steps :sweat_smile:

7

An example can be found at:

This is because I can’t have version dependent signatures (scifem is backwards compatible whenever it can, for ufl expressions I couldn’t find a way to make it work for older versions).

If you use 0.11 or main then you could pass: ufl.dot(u, u) into the interpolate_to_surface_submesh without doing an interpolation into an intermediate space.

8

That worked, and it’s about as fast! Thanks!

Is the goal to eventually just accept an EntityMap?

9

That might be possible to achieve yes. It’s not something I’ve though much about, but if one restricts it to the problem where one wants to interpolate onto the full submesh, that should be doable.