I have also had this problem and spent a lot of time on various forums and I have finally come up with a good solution. This solution will allow you to specify a boundary condition on any whole subdomain within your geometry, and each subdomain can have its own separate boundary condition. I am posting this solution here as the old official Fenics QA forum ...
You need to define your function space of your trial and test functions as the product space of the corresponding elements.
# CG element, degree 1
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
# function space for all fields (here two)
W = FunctionSpace(mesh, P1*P1)
# trial and test functions
u_coolant, u_fuel = TrialFunction(W)
p, q = ...