Adding Dirichlet into NonLinearProblem

Dear all,

I try to use a Newton linearization as a class in python AND implement nonhomogeneous Dirichlet conditions, where am I doing wrong?

Here is the model problem (cleaned so much that it is not a nonlinear problem anymore but anyhow) the boundary in the first iteration is correct but then it is overwritten with zero (for Newton perturbation this is fine) and I want to put it back in every single iteration. I tried in the def F and def J to do it, but no success! Any better ideas? Or should I copy the whole NonlinearProblem class every time :(

#--------------------------------------------------------------------------
from dolfin import *
parameters["form_compiler"]["cpp_optimize"] = True
class iterate(NonlinearProblem):
    def __init__(self, a, L, bc, exter_B):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
 self.bc = bc
 self.exter = exter_B
    def F(self, b, x):
        assemble(self.L, tensor=b, exterior_facet_domains=self.exter)
 for condition in self.bc: condition.apply(b, x)
    def J(self, A, x):
        assemble(self.a, tensor=A, exterior_facet_domains=self.exter)
 for condition in self.bc: condition.apply(A, x)

mesh = Box(0, 0, 0, 1000.0, 1000.0, 1000.0, 20, 50, 20)
V = VectorFunctionSpace(mesh, 'CG', 1)
D = mesh.topology().dim()
bottom = compile_subdomains('x[1] == 0.0')
top = compile_subdomains('x[1] == 1000.0')
# dummy object into the abstract solver
neumann_domains = MeshFunction("uint", mesh, D-1)
neumann_domains.set_all(0)
shear = Expression('5.0')
bc1 = DirichletBC(V.sub(0), shear, top)
bc2 = DirichletBC(V, (0.0,0.0,0.0), bottom)
bc=[bc1,bc2]
dv = TrialFunction(V) # increment
w = TestFunction(V) # test
v = Function(V) # actual velocity from previous iteration
i, j, k, l = indices(4)
G = 80000.0
lmbd = 120000.0
def vij(v):
 return as_matrix(1.0/2.0*(v[i].dx(j) + v[j].dx(i)), [i,j])

S = as_matrix(lmbd*vij(v)[k,k]*Identity(3)[i,j] + 2.0*G*vij(v)[i,j], [i,j])
form = S[i,j]*w[j].dx(i)*dx
gain = derivative(form,v,dv)
problem = iterate(gain, form, bc, neumann_domains)
solver = NewtonSolver('gmres') # for nonsymmetric case
solver.parameters['convergence_criterion'] = 'incremental'
solver.parameters['relative_tolerance'] = 1.0e-1
solver.solve(problem, v.vector())
#----------------------------------------