Evaluate nonlocal functions: how to build a field whose value is taken from the nearest point on the boundary

dolfinx
1

Hi everyone,

I am working on a variational problem in which I need a field defined by nearest-wall evaluation.

Given a domain \Omega and selected boundary facets \Gamma_{\mathrm{wall}} \subset \partial\Omega, I would like to construct a field u_{\mathrm{wall}} on the whole domain such that

u_{\mathrm{wall}}(x) = f(p(x)), \qquad p(x) = \operatorname*{arg\,min}_{y \in \Gamma_{\mathrm{wall}}} |x-y|.

Here p(x) is the nearest point on the wall boundary, and f is a known quantity, available as a UFL expression that I can interpolate into a dolfinx.fem.Function. Ideally, I would only need to define or evaluate f on the boundary, where the normal vector is available.

As far as I understand, this kind of operation is not expressible in pure UFL. At the moment, I approximate u_{\mathrm{wall}} by solving an auxiliary PDE (for example, transport in the wall-normal direction or a Poisson problem).

Is there a way in DOLFINx to obtain or approximate such a field with something like a geometric look-up? I imagine something similar to this.

Illustration of u_{\mathrm{wall}} and f

The following MWE does not reproduce an error. It only illustrates the type of field I want. In this example, the relevant walls are the top and bottom boundaries.

from mpi4py import MPI
import numpy
import pyvista

from dolfinx import mesh, fem, plot

# square: top&bottom = walls; rest doesn't matter
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = fem.functionspace(domain, ("Lagrange", 1))

# function to interpolate into f
def stress(x):
    return 1 + x[0] * 2 * x[1] ** 2

# function to interpolate into u
def wall_stress(x):
    """Return stress value at upper wall for x[1]>0.5 and at lower wall elsewhere."""
    x_mod = x.copy()
    x_mod[1] = numpy.where(x[1] > 0.5, 1.0, 0.0)

    return stress(x_mod)


f = fem.Function(V, name="f")
f.interpolate(stress)

# function that holds f(nearest_wall)
u_wall = fem.Function(V, name="u_wall")
u_wall.interpolate(wall_stress)

# plotting f and u
tdim = domain.topology.dim
domain.topology.create_connectivity(tdim, tdim)
topology, cell_types, geometry = plot.vtk_mesh(V)

for field in [f, u_wall]:
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
    grid.point_data[field.name] = field.x.array.real
    grid.set_active_scalars(field.name)

    plotter = pyvista.Plotter()
    plotter.add_mesh(grid, show_edges=True)
    plotter.view_xy()
    if not pyvista.OFF_SCREEN:
        plotter.show()

Environment

I am using:

  • dolfinx: 0.10.0
  • python: 3.13.9
  • Installation method: conda-forge
  • Platform: macOS 26.3.1, arm64

For background information, I need to evaluate wall shear stresses in a turbulence model I am implementing. Many thanks for your consideration!

Best,

Kai

2

I think you could do this with FEniCSx_ii, by making a modification to the mapping operator described in: