"""
Douglas N. Arnold, 2011-12-01, 2014-12-03

Solve the Dirichlet problem for the minimal surface equation.

-div( q grad u) = 0,

where q = q(grad u) = (1 + |grad u|^2)^{-1/2}

using Newton's method
"""

from fenics import *

# Create mesh and define function space
mesh = UnitSquareMesh(100, 100)
V = FunctionSpace(mesh, 'Lagrange', 2)

# initial guess (its boundary values specify the Dirichlet boundary conditions)
# (larger coefficient in front of the sin term makes the problem "more nonlinear")
u0 = Expression('a*sin(2.5*pi*x[1])*x[0]', a=.2, degree=5)
u = interpolate(u0, V)
plot(u)
print 'initial surface area:', assemble(sqrt(1+inner(grad(u), grad(u)))*dx)

# Define the linearized weak formulation for the Newton iteration
du = TrialFunction(V)
v = TestFunction(V)
q = (1 + inner(grad(u), grad(u)))**(-.5)
a =  q * inner(grad(du), grad(v)) * dx - q**3 * inner(grad(u), grad(du)) * inner(grad(u), grad(v)) * dx
L = -q*inner(grad(u), grad(v))*dx
du = Function(V)

# Newton iteration
tol = 1.0E-5
maxiter = 30
for iter in range(maxiter):
    # compute the Newton increment by solving the linearized problem;
    # note that the increment has *homogeneous* Dirichlet boundary conditions
    solve(a == L, du, DirichletBC(V, 0.0, DomainBoundary()))
    u.vector()[:] += du.vector() # update the solution
    eps = sqrt(abs(assemble(inner(grad(du), grad(du))*dx))) # check increment size as convergence test
    area = assemble(sqrt(1+inner(grad(u), grad(u)))*dx)
    print('iteration{:3d}  H1 seminorm of delta: {:10.2e}  area: {:13.5e}'.format(iter+1, eps, area))
    if eps < tol:
      break

if eps > tol:
  print('no convergence after {} Newton iterations'.format(iter+1))
else:
  print('convergence after {} Newton iterations'.format(iter+1))

plot(u, interactive=True)