It appears to me that passing in boundary conditions in the solve command when using pre-assembled matrices is not working correctly.  Consider the following problem:

**************************************************

mesh = UnitSquareMesh(10,10)

V = FunctionSpace(mesh,'CG',1)
x,y = SpatialCoordinate(mesh)
f_ex = 10.*exp(-(x-0.5)**2/0.02 - (y-0.5)**2/0.02)
f = interpolate(f_ex,V)
g_ex = sin(5*x)
g = interpolate(g_ex,V)
u, v = TrialFunction(V), TestFunction(V)
bc = DirichletBC(V,0.0,(1,2))

A = assemble(inner(grad(u),grad(v))*dx,bcs=bc)
b = assemble(f*v*dx + g*v*ds((3,4)),bcs=bc)
U_1 = Function(V)
solve(A,U_1,b)

A = assemble(inner(grad(u),grad(v))*dx)
b = assemble(f*v*dx + g*v*ds((3,4)))
U_2 = Function(V)
solve(A,U_2,b,bcs=bc)

A = inner(grad(u),grad(v))*dx
b = f*v*dx + g*v*ds((3,4))
U_3 = Function(V)
solve(A==b,U_3,bcs=bc)

****************************************************

The solve commands for U_1 and U_3 both converge, and the solutions agree to within machine precision.  However, Firedrake fails to solve for U_2, with the error message:

ConvergenceError: ('LinearSolver failed to converge after %d iterations with reason: %s', 1320, 'DIVERGED_DTOL'

I suspect that boundary conditions are not being applied in the second solve command to the pre-assembled system.  But according to the documentation on firedrakeproject.org,
this should work as coded above.

Am I doing something incorrectly?  I'm also wondering how this would affect LinearVariationalSolver if there is indeed a problem with the boundary conditions.

Thanks,

Peter Sentz