Neumann boundary conditions are not satisfied

Hi folks,

I have a quick question about Neumann BCs. I am trying to solve a nonlinear PDE with only Neumann Bcs at boundaries. Instead of Dirichlet BC, I am trying to impose an integral constraint. One may interpret this problem as a nonlinear version of demo_neumann-poisson.py.
My problem is that the Neumann BCs are not satisfied, which results in an incorrect solution. I would be grateful if someone could take a look at my code and see why the Neumann BCs are not satisfied.

Thank you for your help.
Mo

A simplified version of my code:

from dolfin import *

# Inputs
epsa = 0.485
epsfa = 0.0326
temp = 298.0
FF = 96487.0
RR = 8.314
iex = 1.65090528134
sa = 723600.0
iapp = 30
c0 = 1000
sig = 100.0
kappa = (4.1253E-2 + 5.007E-4*c0 - 4.7212E-7*(c0**2) + 1.5094E-10*(c0**3) - 1.6018E-14*(c0**4))
kapae = (epsa**4.0)*kappa
sige = sig*(1.0 - epsa - epsfa)
sigha = (sige*kapae)/(sige + kapae)
coef = FF/(2.0*RR*temp)

# Create mesh and define function space
nel = 64
Lx = 88.0E-6
mesh = Interval(nel, 0, Lx)
V = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R

# Define variational problem
v_a = TestFunction(W)
dwa = TrialFunction(W)
wa = Function(W)

v, d = split(v_a)
u, c = split(wa)

# Initialize the Newton solver (0.5*iapp/kapae)
g = Expression("(x[0]/Lx)*(0.5*iapp/kapae) - ((Lx - x[0])/Lx)*(iapp/sige)", Lx=Lx, iapp=iapp, kapae=kapae, sige=sige)
j = iex*sa*(exp(coef*u) - exp(-coef*u))
#j = iex*sa*u/(RR*temp)
a = -inner(sigha*grad(u), grad(v))*dx + c*v*dx - j*v*dx + j*d*dx + g*v*ds - iapp*d*dx
solve(a == 0, wa, solver_parameters={"newton_solver":{"relative_tolerance": 1e-14, "maximum_iterations": 200}})

# Solve the problem
gg = Expression("(x[0]/Lx)*(iapp/kapae) - ((Lx - x[0])/Lx)*(iapp/sige)", Lx=Lx, iapp=iapp, kapae=kapae, sige=sige)
aa = -inner(sigha*grad(u), grad(v))*dx + c*v*dx - j*v*dx + j*d*dx + gg*v*ds - iapp*d*dx
solve(aa == 0, wa, solver_parameters={"newton_solver":{"relative_tolerance": 1e-14, "maximum_iterations": 200}})

(u, c) = wa.split()
uu = project(u, V)

# Gradient of the solution
g_u = project(grad(uu), V)

print "Slope at x=0:::>", g_u.vector().array()[0], " Should be equal to -iapp/sige:::>", -iapp/sige
print "Slope at x=Lx::>:", g_u.vector().array()[nel], " Should be equal to iapp/kapae:::>", iapp/kapae

ja = Function(V)
ja = iex*sa*(exp(coef*uu) - exp(-coef*uu))
print "\int j_vol=", assemble(ja*dx), "Should be equal to iapp=", iapp

# Plot solution
plot(ja)
plot(u)
interactive()