### Custom Newton solver for nonlinear problem and AttributeError: 'list' object has no attribute 'homogenize'

**First,**there is an error with applying boundary condition to the incremental du:

`AttributeError: 'list' object has no attribute 'homogenize'`

This error is related with the code:

`bcs_du = bcs.homogenize()`

where bcs is a list of DirichletBC.

**Second,**if I manually specify the BCs of du, instead of using

`bcs_du = bcs.homogenize()`

, i.e.,`bcs_du = [ DirichletBC(V.sub(0), 0, left),`

`DirichletBC(V.sub(1), 0, bot),`

`DirichletBC(V.sub(1), 0, top) ]`

The error disappears.

*However, the result from the the custom Newton solver is different from the default NonlinearVariationalSolver.*The minimum working code is as follows. Note that currently the code is dealing with a constant body force B(u) = [0 , 1], I am using a nonlinear solver to solve a linear elasticity problem. A nonlinear problem appears if the body force vector B(u) indeed depends on u, e.g., B(u) = [u.sub(0)^2, u.sub(1)^2 ]. The nonlinear case is actually what I want to do. But first I need to fix the problems mentioned before.

`from dolfin import *`

`import os, shutil, math, sys, sympy`

`import numpy as np`

`# result folder`

`savedir = "Result/Example"`

`if os.path.isdir(savedir):`

` shutil.rmtree(savedir) `

`# Create mesh and define function space`

`mesh = RectangleMesh(Point(0.0, 0.0), Point(1.0, 2.0), 25, 50, "right/left")`

`d = mesh.topology().dim()`

`e1 = [Constant([1.,0.]),Constant((1.,0.,0.))][d-2]`

`e2 = [Constant([0.,1.]),Constant((0.,1.,0.))][d-2]`

`# Mark boundary subdomains`

`class left(SubDomain): `

` def inside(self,x,on_boundary): `

` return near(x[0], 0.)`

`class bot(SubDomain):`

` def inside(self,x,on_boundary):`

` return near(x[1], 0.)`

`class top(SubDomain):`

` def inside(self,x,on_boundary):`

` return near(x[1], 2.)`

`left, bot, top = left(), bot(), top()`

`boundaries = FacetFunction("uint", mesh)`

`boundaries.set_all(0)`

`left.mark(boundaries, 1)`

`bot.mark(boundaries, 2)`

`top.mark(boundaries, 3)`

`ds = Measure("ds", subdomain_data=boundaries)`

`# Define functions`

`V = VectorFunctionSpace(mesh, "Lagrange", 1)`

`du = TrialFunction(V) # Incremental displacement`

`v = TestFunction(V) # Test function`

`u = Function(V) # Displacement from previous iteration`

`# Define Dirichlet boundary`

`bc_l = DirichletBC(V.sub(0), 0, left)`

`bc_b = DirichletBC(V.sub(1), 0, bot)`

`bc_t = DirichletBC(V.sub(1), 0.2, top)`

`bcs = [ bc_l, bc_b, bc_t ]`

`# Elasticity parameters`

`mu, lmbda = 1.0, 1.0`

`# Strain and stress`

`def eps(u):`

` return sym(grad(u))`

`def sigma(u):`

` return 2.0*mu*(eps(u)) + lmbda*tr(eps(u))*Identity(d)`

`def B(u):`

` uy = u.sub(1)`

` return as_vector( [ 0.0, 1 ] )`

`T = Constant((0.0, 0.0)) # Traction force on the boundary`

`file_u = File(savedir+"/u.pvd");`

`file_s = File(savedir+"/stress.pvd")`

`# Total potential energy`

`Pi = 0.5*inner(sigma(u),eps(u))*dx - dot(B(u), u)*dx - dot(T, u)*ds`

`# Compute first variation of Pi (directional derivative about u in the direction of v)`

`F = derivative(Pi, u, v)`

`# Compute Jacobian of F`

`Jac = derivative(F, u, du)`

`u_init = Expression(("0", "0"), degree=1)`

`u.interpolate(u_init)`

`abs_tol, rel_tol = 1e-10, 1e-8`

`abs_res, rel_res = 1.0, 1.0`

`max_iter = 100`

`# bcs_du = [ DirichletBC(V.sub(0), 0, left),`

`# DirichletBC(V.sub(1), 0, bot),`

`# DirichletBC(V.sub(1), 0, top) ]`

`bcs_du = bcs.homogenize()`

`u_inc = Function(V)`

`nIter = 0`

`flag = bool(int(sys.argv[1]))`

`if flag:`

print ("Default solver")
` problem = NonlinearVariationalProblem(F, u, bcs, Jac)`

` solver = NonlinearVariationalSolver(problem)`

` solver.solve()`

`else:`

print ("Custom solver")
` while abs_res > abs_tol and nIter < max_iter: # Newton iterations`

` nIter += 1`

` abs_res_old = abs_res`

` A, b = assemble_system(Jac, -F, bcs_du)`

` solve(A, u_inc.vector(), b) # Determine step direction`

` u_inc_norm = np.linalg.norm(u_inc.vector().array(), ord = 2)`

` a = assemble(F)`

` for bc in bcs_du:`

` bc.apply(a)`

` abs_res = b.norm('l2')`

` rel_res = abs_res/abs_res_old`

` alpha = 1.0 # step size, initially 1`

` u.vector()[:] += alpha*u_inc.vector() # New u vector`

` string = " Newton iteration %d: r (abs) = %.3e (tol = %.3e) " \`

` +"r (rel) = %.3e (tol = %.3e)"`

` print string % (nIter, abs_res, abs_tol, rel_res, rel_tol)`

`file_u << u`

### 1 Answer

```
bc0 = [DirichletBC(bc) for bc in bcs] # List of copies of your original bcs
[bc.homogenize() for bc in bc0] # homogenize every item (DirichletBC object)
```

The list bc0 contains the homogenized versions of your bcs.

An idea as to why your residuals do not match:

The relative residual in the FEniCS NewtonSolver method is the residual in the current iteration divided by the residual in the zeroth iteration, before anything is solved. So basically the residual, your initial guess produces. If I read your code correctly you define your relative residual as the ratio of the current residual and the residual from the previous iteration.

https://fenicsproject.org/qa/536/newton-method-programmed-manually

are all homogeneous, but it should have been included anyway because the use of homogenize really doesn't make any sense without it.

You have to apply the boundary condition to the unknown before entering newton, because the direction of decent is set to zero along Dirichlet boundaries, meaning that the value of the unknown

*will not change*as it is incremented down the gradient.

` bcs_du = [ DirichletBC(V.sub(0), 0, left),`

` DirichletBC(V.sub(1), 0, bot),`

` DirichletBC(V.sub(1), 0, top) ]`

and

`bcs_du = [DirichletBC(bc) for bc in bcs] # List of copies of your original bcs`

`[bc.homogenize() for bc in bcs_du] # homogenize every item`

.

Regarding to the second question, the criterion of relative residual is not used for the convergence check in my code, although I computed it ("incorrectly". Thank you for pointing out how the relative residual is computed in Fenics, I do not know that before). I apply the homogenized boundary conditions to du, as you indicated. However, the result of the Newton solver is still not correct. So I think the mismatach of the two solvers are not due to the relative residual convergence check.

I found the mismatch of the result from the two solver disappears if I delete the non-zero boundary condition bc_t on the top boundary. They match perfectly. But if I add the bc_t boundary, the result of the Newton solve is not correct. I have no idea why this happens. Could you help me with this?

First I would delete the lines

```
u_init = Expression(("0", "0"), degree=1)
u.interpolate(u_init)
```

since the initialization of u already writes zeros too every entry.

Then add the line

`[bc.apply(u.vector()) for bc in bcs]`

- where bcs are your original non-homogenous Dirichlet boundary conditions - somewhere before your iteration loop.

The only reason you have to set homogenous Dirichlet boundary conditions on your increment is that you suppose that your initial guess already satisfies all Dirichlet boundary conditions. When the initial guess gets incremented in a Newton iteration, you'll want to add a zero in every entry where Dirichlet boundary conditions are prescribed otherwise you would change the values that originally satisfied the BC - thus homogenous DirichletBC on the increment. But for this to work, your initial u vector has to satisfy the Dirichlet BC in the first place, hence the line above.

I change my initial guess to

`u_init = Expression(("0", "0.1*x[1]"), degree=1)`

`u.interpolate(u_init)`

to meet the boundary conditions.

```
bcs_du = [ DirichletBC(V.sub(0), 0, left),
DirichletBC(V.sub(1), 0, bot),
DirichletBC(V.sub(1), 0, top) ]
```

is an explicitly given special case, tailored to the non-homogenous BC

```
# Define Dirichlet boundary
bc_l = DirichletBC(V.sub(0), 0, left)
bc_b = DirichletBC(V.sub(1), 0, bot)
bc_t = DirichletBC(V.sub(1), 0.2, top)
bcs = [ bc_l, bc_b, bc_t ]
```

whereas the expression I suggested works for any case of Dirichlet BC. Thus I'd prefer the second variant so I can simply change the non-homogenous BC and the homogenized version adapts accordingly.