Integral Constraints on PDEs with FEniCSx - Real Finite Elements

In a problem I stumbled upon I have to incorporate some sort of non trivial integral constrain condition.

I tried to start from a simple example, as shown here fenics code. It is a calculation in the legacy FEniCS. However, I find it hard to repeat it in FEniCSx. Seems like the problem is in the Real Finite Element. Some minimal not working code follows:

from mpi4py import MPI
import dolfinx
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType

domain = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10)

element1 = ufl.FiniteElement("CG", domain.ufl_cell(), 1)
element2 = ufl.FiniteElement("Real", domain.ufl_cell(), 0)
mel = ufl.MixedElement([element1, element2])
V = dolfinx.fem.FunctionSpace(domain, mel)
V0, submap = V.sub(0).collapse()

# Test functions
v, r = ufl.TestFunctions(V)
# Fields to solve for
u = dolfinx.fem.Function(V)
u1, u2 = u.split()
u.name = "u"

def left_boundary(x):
    return np.isclose(x[0], 0)
def right_boundary(x):
    return np.isclose(x[0], 1)

# Determine boundary DOFs
boundary_dofs = dolfinx.fem.locate_dofs_geometrical((V.sub(0),V0), left_boundary)
bc_l = dolfinx.fem.dirichletbc(ScalarType(0), boundary_dofs[0], V.sub(0))
boundary_dofs = dolfinx.fem.locate_dofs_geometrical((V.sub(0),V0), right_boundary)
bc_r = dolfinx.fem.dirichletbc(ScalarType(0), boundary_dofs[0], V.sub(0))

# Functions
C0 = dolfinx.fem.Constant(domain, ScalarType(1))
f = dolfinx.fem.Constant(domain, ScalarType(0))

a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx + (u - C0) * r * ufl.dx
L = f * v * ufl.dx
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc_r])
uh = problem.solve()
print(u1.x.array)
print(u2.x.array)