Minimal surface solution in Python

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?

• All you have is this boundary? May 19 '13 at 12:56
• @J.M. Yes, all I have is this set of points defining the boundary. May 19 '13 at 12:57
• Is there a reason you want the minimal surface equation solution as opposed to some other one? There are infinitely many surfaces that you can create that have this curve as their boundary. May 19 '13 at 13:50
• I would add that to the description of the question. May 19 '13 at 14:47
• Do you have your flow field? Is your flow incompressible? Depending on the answer to these, you might not need the minimal surface to estimate the flow. any surface will do. Dec 19 '19 at 16:01

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.

• Bill thanks for all the resources! I haven't heard of FEniCS, and doing some search I found this link: fenicsproject.org/documentation/tutorial/…. Would this be the approach for solving the PDF in python? May 19 '13 at 15:35
• That seems like one approach, but I'm not a FEniCS expert. May 19 '13 at 15:49

You can use FEniCS:

from fenics import (
UnitSquareMesh,
FunctionSpace,
Expression,
interpolate,
assemble,
sqrt,
inner,
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)*x", a=0.2, degree=5)
u = interpolate(u0, V)
print(
)

# Define the linearized weak formulation for the Newton iteration
du = TrialFunction(V)
v = TestFunction(V)
a = (
)
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(
)  # check increment size as convergence test
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)
• This code uses the undamped Newton method, which will converge if the initial guess is sufficiently close to the solution, but which can also diverge if the initial guess is not close enough. To work around this, you can use a line search strategy. I've shown how to do this here using Firedrake instead of FEniCS, but they're quite similar. Dec 19 '19 at 19:00