Minimal surface solution in Python

Question

Note: this question was also posted in StackOverflow and math.stackexchange.

I have a set of 3D points defining a 3D contour, as shown below. The points in this contour lie in their best-fit plane and I want to obtain a 3D triangular mesh representation of the surface inside this contour.

Doing some research I found that this is basically a minimal surface problem and its solution is related with the Biharmonic Equation. I also found that the Thin-plate spline is the fundamental solution to this equation.

So I think the approach would be to try to fit this sparse representation of the surface (given by the 3D contour of points) using thin-plate splines. I found this example in scipy.interpolate where scattered data (x,y,z format) is interpolated using thin-plate splines to obtain the ZI coordinates on a uniform grid (XI,YI).

Two questions arise:

1. Would thin-plate spline interpolation be the correct approach for the problem of computing the surface from the set of 3D contour points?

2. If so, how to perform thin-plate interpolation on scipy with a non-uniform grid?

Answer 1

I don't think you can just interpolate. You may know the Fundamental Solution to the governing equations, but you have no data about the interior points of the problem, and I suspect that the solution will be poor far from the boundary.

Also, it's been awhile, but I usually associate the following equation with minimal surfaces (for small displacements, at least): $$ (1+u_x^2)u_{yy} - 2u_x u_y u_{xy} + (1+u_y^2)u_{xx} = 0 $$ where $u(x,y)$ is the small displacement of the surface in question.

You could solve this by projecting your boundary points down into the $(x,y)$-plane and using their $z$ value as $u(x,y)$ on the boundary. Using the $(x,y)$ coordinates of the boundary points as the boundary of your domain, you could create a triangular mesh on the interior and then solve the above equation with above boundary conditions.

deal.II's example #15 does this, though I think you'd need a quadrilateral mesh rather than triangles. I can't tell from the figure, but if your boundary is a nominally a circle, then generating a quad mesh for it should be straightforward. If it's arbitrary, then you'll have to do some more work.

There might also be a FEniCS example for this if you prefer Python to C++, but I couldn't tell from a quick search around the web. There's also a quick MATLAB example using the PDE Toolbox if you have access to that. It will even do the mesh generation for you.

Answer 2

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.)