BUG: Boundary conditions not applied by LinearVariationalSolver when the right-hand-side is a Cofunction or FormSum

[jrmaddison]

Describe the bug
Boundary conditions are not applied in a linear variational solve when the right-hand-side form is a Cofunction or FormSum.

Steps to Reproduce

from firedrake import *

mesh = UnitIntervalMesh(2)
space = FunctionSpace(mesh, "Lagrange", 1)
test, trial = TestFunction(space), TrialFunction(space)

# Form RHS

u = Function(space, name="u")
problem = LinearVariationalProblem(
    inner(trial, test) * dx, inner(Constant(1.0), test) * dx, u,
    DirichletBC(space, 0.0, "on_boundary"))
solver = LinearVariationalSolver(problem)
solver.solve()

print(f"{assemble(inner(u, u) * ds)=}")  # Zero

# Cofunction RHS

b = assemble(inner(Constant(1.0), test) * dx)
u = Function(space, name="u")
problem = LinearVariationalProblem(
    inner(trial, test) * dx, b, u,
    DirichletBC(space, 0.0, "on_boundary"))
solver = LinearVariationalSolver(problem)
solver.solve()

print(f"{assemble(inner(u, u) * ds)=}")  # Nonzero

# FormSum RHS

b = assemble(inner(Constant(0.5), test) * dx) + inner(Constant(0.5), test) * dx
u = Function(space, name="u")
problem = LinearVariationalProblem(
    inner(trial, test) * dx, b, u,
    DirichletBC(space, 0.0, "on_boundary"))
solver = LinearVariationalSolver(problem)
solver.solve()

print(f"{assemble(inner(u, u) * ds)=}")  # Nonzero

Expected behavior
Boundary conditions should be applied in the solve. In the Cofunction case dofs corresponding to boundary evaluation should be ignored.

Error message
No error.

Environment:
Ubuntu 22.04.

Additional Info
The Cofunction case appears e.g. when defining a Riesz map on zero trace spaces.