Nonlocal (boundary) terms in interior form due to rescaling

I’m trying to solve a continuity equation

\frac{\partial u}{\partial t} = - \frac{\partial \mathbf{j}_u}{\partial x}

in 1D on a growing domain \Omega from 0 to x_l(t). The growth speed of the domain depends on the value of u at the boundary x_l,

\frac{\mathrm{d} x_l}{\mathrm{d}t} = f(x|_{x_l}).

To avoid remeshing after every timestep, I want to scale the domain as \xi = \frac{x}{x_l(t)}. This results in the following equation:

\frac{\partial u}{\partial t} = - \frac{1}{x_l} \frac{\partial \mathbf{j}_u}{\partial \xi} + \frac{\xi}{x_l} \frac{\mathrm{d} x_l}{\mathrm{d}t} \frac{\partial u}{\partial \xi}.

If I insert the equation for \frac{\mathrm{d} x_l}{\mathrm{d} t}, I will end up with a nonlocal boundary term in the interior form,

\int v \frac{\xi}{x_l} f(u|_{\xi = 1}) \frac{\partial u}{\partial \xi} \mathrm{d} \Omega,

which I don’t know how to represent with the UFL.

What I currently do is represent the boundary value f(u|_{x_l}) as a separate constant, find the boundary DoF and update the constant manually before every timestep.

u = Function(V)
f = Constant(mesh, 0.0)
boundary_dof = locate_dofs_geometrical(...)
f.value = my_function(u.x.array[boundary_dof])

But this will obviously not give me the correct bilinear form when I use the automatic functional differentiation of UFL.

What is the best way to handle these sort of terms in FEniCSx?