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.

enter image description here

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?

  • $\begingroup$ All you have is this boundary? $\endgroup$ – J. M. May 19 '13 at 12:56
  • $\begingroup$ @J.M. Yes, all I have is this set of points defining the boundary. $\endgroup$ – CodificandoBits May 19 '13 at 12:57
  • 1
    $\begingroup$ 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. $\endgroup$ – Bill Barth May 19 '13 at 13:50
  • 3
    $\begingroup$ I would add that to the description of the question. $\endgroup$ – Bill Barth May 19 '13 at 14:47
  • 1
    $\begingroup$ 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. $\endgroup$ – MPIchael 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.

| cite | improve this answer | |
  • $\begingroup$ 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? $\endgroup$ – CodificandoBits May 19 '13 at 15:35
  • $\begingroup$ That seems like one approach, but I'm not a FEniCS expert. $\endgroup$ – Bill Barth May 19 '13 at 15:49

You can use FEniCS:

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=0.2, degree=5)
u = interpolate(u0, V)
    "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)
        f"iteration{iter + 1:3d}  H1 seminorm of delta: {eps:10.2e}  area: {area:13.5e}"
    if eps < tol:

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

with XDMFFile(MPI.comm_world, "out.xdmf") as xdmf_file:

(Modified from http://www-users.math.umn.edu/~arnold/8445/programs/minimalsurf-newton.py.)

| cite | improve this answer | |
  • 2
    $\begingroup$ 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. $\endgroup$ – Daniel Shapero Dec 19 '19 at 19:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.