# Projecting an expression onto a boundary submesh

**URL:** <https://fenicsproject.discourse.group/t/projecting-an-expression-onto-a-boundary-submesh/19772>\
**Category:** dolfinx\
**Created:** [July 14, 2026, 8:01pm UTC](https://fenicsproject.discourse.group/t/projecting-an-expression-onto-a-boundary-submesh/19772 "2026-07-14T20:01:38Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![rmsc](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/rmsc/32/9586_2.png) [@rmsc](https://fenicsproject.discourse.group/u/rmsc)\
**Post date:** [July 14, 2026, 8:01pm UTC](https://fenicsproject.discourse.group/t/projecting-an-expression-onto-a-boundary-submesh/19772/1 "2026-07-14T20:01:38Z")

</div>

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:

```python
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:

```bash
[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!

---

<div class="post-metadata">

**Author:** ![andreaantonelli1](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/andreaantonelli1/32/9291_2.png) [@andreaantonelli1](https://fenicsproject.discourse.group/u/andreaantonelli1)\
**Post date:** [July 15, 2026, 6:49am UTC](https://fenicsproject.discourse.group/t/projecting-an-expression-onto-a-boundary-submesh/19772/2 "2026-07-15T06:49:08Z")

</div>

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

```auto
# 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](https://fenicsproject.discourse.group/t/interpolate-expression-on-submesh-with-parent-mesh-function/17467)

---

<div class="post-metadata">

**Author:** ![rmsc](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/rmsc/32/9586_2.png) [@rmsc](https://fenicsproject.discourse.group/u/rmsc)\
**Post date:** [July 15, 2026, 10:52am UTC](https://fenicsproject.discourse.group/t/projecting-an-expression-onto-a-boundary-submesh/19772/3 "2026-07-15T10:52:07Z")

</div>

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:

```python
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)

```

```bash
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:

```python
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)

```

```bash
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

```

---

<div class="post-metadata">

**Author:** ![rmsc](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/rmsc/32/9586_2.png) [@rmsc](https://fenicsproject.discourse.group/u/rmsc)\
**Post date:** [July 15, 2026, 11:17am UTC](https://fenicsproject.discourse.group/t/projecting-an-expression-onto-a-boundary-submesh/19772/4 "2026-07-15T11:17:44Z")

</div>

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

```python
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!!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [July 18, 2026, 11:21am UTC](https://fenicsproject.discourse.group/t/projecting-an-expression-onto-a-boundary-submesh/19772/5 "2026-07-18T11:21:41Z")

</div>

> [@rmsc](#):
>
> ```auto
> # 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)
> 
> ```

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.

> <https://github.com/scientificcomputing/scifem/blob/7c0b2cb69144e37cae5a0a8e78dc75c84284a5d9/src/scifem/interpolation.py#L332>

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

---

<div class="post-metadata">

**Author:** ![rmsc](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/rmsc/32/9586_2.png) [@rmsc](https://fenicsproject.discourse.group/u/rmsc)\
**Post date:** [July 19, 2026, 10:14am UTC](https://fenicsproject.discourse.group/t/projecting-an-expression-onto-a-boundary-submesh/19772/6 "2026-07-19T10:14:16Z")

</div>

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`..

> [@dokken](#):
>
> This works for general `ufl.core.expr.Expr` since v0.11.0

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 😅

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [July 19, 2026, 3:28pm UTC](https://fenicsproject.discourse.group/t/projecting-an-expression-onto-a-boundary-submesh/19772/7 "2026-07-19T15:28:09Z")

</div>

> [@rmsc](#):
>
> Thanks, but to be very honest, I just couldn’t figure out how to use `interpolate_to_surface_submesh`

An example can be found at:

> <https://github.com/scientificcomputing/scifem/pull/226#issue-4586778760>
>
> Following the updates in DOLFINx (https://github.com/FEniCS/dolfinx/pull/4140) w…e can now evaluate expressions (not backwards compatible) in interpolate\_to\_surface\_submesh.
> Additional support for sub-meshes through entity maps.
> 
> \`\`\`python
> 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
> \`\`\`bash
> 1.4704992620052726 \[0.823 -0.2 \]
> \`\`\`
> \<img width="2272" height="2218" alt="image" src="https://github.com/user-attachments/assets/7fd38ea9-11ba-4acb-840b-fdccd93c3ac5" /\>

> [@rmsc](#):
>
> You mean in `u_volume`? The function signature expects a `fem.Function`..

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.

---

<div class="post-metadata">

**Author:** ![rmsc](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/rmsc/32/9586_2.png) [@rmsc](https://fenicsproject.discourse.group/u/rmsc)\
**Post date:** [July 20, 2026, 11:42am UTC](https://fenicsproject.discourse.group/t/projecting-an-expression-onto-a-boundary-submesh/19772/8 "2026-07-20T11:42:26Z")

</div>

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

Is the goal to eventually just accept an EntityMap?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [July 20, 2026, 12:09pm UTC](https://fenicsproject.discourse.group/t/projecting-an-expression-onto-a-boundary-submesh/19772/9 "2026-07-20T12:09:51Z")

</div>

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.