Imposing weak boundary condition

I have adapted the demo "dg-advection-diffusion.py" to solve the following steady-state advection equation with Dirichlet BC on the inflow boundary:

          div(uw) = f (x,y in [0,1]x[0,1]),

          w = 1 (x = 0)

where we can assume w = x**2 + 1 and u = (1,0) (velocity), which implies f = 2x.

The code can be found below. Before reading it, here are the problems I am getting:

1) If I don't specify any BC, the solution I get is w = x**2 and not w = x**2 + 1, so it differs by a constant. Moreover, the error (L2 norm) is 1.00000000136, which becomes 1.36E-09 if we subtract it by 1 (the difference of the solutions).

2) If I impose a strong BC on x = 0, I get w = x**2 + 1, however, the error is now 1.8E-04, which is far larger than the above one and it does not improve when the mesh is refined.

So, how can I solve this problem by specifying proper BC that will result in more precise solutions? What about specifying the Dirichlet BC weakly (i.e. as a flux in the inflow boundary)? How can I do that?

==========================================================
from dolfin import *

# Exact solution
w_e = Expression('x[0] * x[0] + 1')

# Load mesh
mesh = UnitSquare(40, 40)

# Define inflow boundary
class InflowBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0]) < DOLFIN_EPS and on_boundary

# Defining the function spaces
V_dg = FunctionSpace(mesh, "DG", 1)
V_cg = FunctionSpace(mesh, "CG", 1)
V_u = VectorFunctionSpace(mesh, "CG", 1)

# Create velocity Function
u = Constant((1,0))

# Test and trial functions
w = TrialFunction(V_dg)
v = TestFunction(V_dg)

# Source term
f = Expression('2 * x[0]')

# Mesh-related functions
n = FacetNormal(mesh)

# (dot(v, n) + |dot(v, n)|)/2.0
un = (dot(u, n) + abs(dot(u, n)))/2.0

# Bilinear form
a = -dot(grad(v), u*w)*dx + \
     dot(jump(v), un('+')*w('+') - un('-')*w('-') )*dS + \
     dot(v, un*w)*ds

# Linear form
L = v*f*dx

# Set up boundary condition (apply strong BCs)
bc = DirichletBC(V_dg, w_e, InflowBoundary(), "geometric")

#bc = None # Uncomment this line to get w = x**2

# Solution function
w_h = Function(V_dg)

# Assemble linear system
A, b = assemble_system(a, L, bc)

# Solve system
solve(A, w_h.vector(), b)

# Project solution to a continuous function space
w_p = project(w_h, V_cg)

# Print the error
print "Error (L2): ", errornorm(w_e, w_h)

# Plot solution
plot(w_p, interactive=True)