You can use FEniCS:
from fenics import (
UnitSquareMesh,
FunctionSpace,
Expression,
interpolate,
assemble,
sqrt,
inner,
grad,
dx,
TrialFunction,
TestFunction,
Function,
solve,
DirichletBC,
DomainBoundary,
MPI,
XDMFFile,
)
# 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=0.2, degree=5)
u = interpolate(u0, V)
print(
"initial surface area: {}".format(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))) ** (-0.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(
f"iteration{iter + 1:3d} H1 seminorm of delta: {eps:10.2e} area: {area:13.5e}"
)
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))
with XDMFFile(MPI.comm_world, "out.xdmf") as xdmf_file:
xdmf_file.write(u)
(Modified from http://www-users.math.umn.edu/~arnold/8445/programs/minimalsurf-newton.py.)