Error: Expecting the left-hand side to be a bilinear form (not rank 1)

[Mohammad Babar]

Hi everyone, 

I am solving an approximate steady state Poisson-Nernst-Plank equation,

Here is my attempt on a simple square grid,
```

P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

# Define boundary condition
c_D = Expression('0 + 0*x[0]*x[0] + 0*x[1]*x[1]', degree=2)
c_L = Expression('1 + 0*x[0]*x[0] + 0*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
    return on_boundary

q_1, q_2 = TestFunctions(V)
c = Function(V)
c_1, c_2 = split(c)

# Constants
z_1 = 1 # Charge number
z_2 = -1
eps = 80
eps0 = 8.854187e-12 # Vacuum permittivity F/m
e = 1.60218e-19 # Coulombs
kbT = 0.0259 # eV, room temp
F = 96485.332 # C/mol
a = (1/kbT)
b = a*(F/(eps*eps0))
a, b = 1, 1

r = Expression('x[0]', degree=1)
nabla_phi = (z_1*c_1 + z_2*c_2)

F = (dot(grad(c_1), grad(q_1)))*r*dx() \
    + (dot(grad(c_2), grad(q_2)))*r*dx() \
    + ((z_1*c_1*q_1 + z_2*c_2*q_2)*b*nabla_phi)*r*dx()
   
g = Expression('1.0', degree=2)

tol = 1E-14
#L = g*q_1*ds + g*q_2*ds
L = 0
def boundary_z(x, on_boundary):
    return on_boundary and (near(x[1], 1, tol))

def boundary_r(x, on_boundary):
    return on_boundary and near(x[0], 1, tol)

bc_z= DirichletBC(V.sub(0), c_D, boundary_z)
bc_z= DirichletBC(V.sub(0), c_L, boundary_z)
bc_r= DirichletBC(V.sub(1), c_L, boundary_r)
bcs = [bc_r, bc_z]

# Compute solution
solve(F == L, c, bcs)
_c_1, _c_2 = c.split()

# Plot solution and mesh
plot(_c_1)
plot(mesh)

```

It works fine if L = 0, however if I want to impose a Neumann BC with 

L = g*q_1*ds + g*q_2*ds,

It throws me this error,

```

*** Error: Unable to define linear variational problem a(u, v) == L(v) for all v.

 *** Reason: Expecting the left-hand side to be a bilinear form (not rank 1).

 *** Where: This error was encountered inside LinearVariationalProblem.cpp. 

*** Process: 0 *** *** DOLFIN version: 2019.1.0 

*** Git changeset: 

```

Here is the full traceback,

Thanks

-Mohammad