Given a closed curve $\mathcal C$ in three dimensions, is it possible to use Mathematica's built-in functionality to compute a minimal surface whose boundary is $\mathcal C$? For simplicity, let us assume the surface to be a topological disk.

We could choose a domain $U\subset\mathbb R^2$, say the unit disk or the square $[0,1]\times[0,1]$, and take the unknown surface $\mathcal S$ and the given curve $\mathcal C$ to be parametrized by $U$ and its boundary $\partial U$ respectively. That is, we specify $\mathcal C$ as the image of a function $g:\partial U\to\mathbb R^3$, and seek a function $f:U\to\mathbb R^3$ that satisfies the boundary condition $f=g$ on $\partial U$, and such that the surface $\mathcal S=f(U)$ has zero mean curvature everywhere.

This seems a lot like some of the problems that the new FEM functionality in NDSolve can handle. But it's highly nonlinear, so maybe not.

Here's what I've tried so far; maybe it can help you get started. We'll use J.M.'s implementation of mean curvature, and try to recover Scherk's first surface $\exp z=\cos x/\cos y$ in the region $-1\le x\le1$, $-1\le y\le1$.

region = Rectangle[{-1, -1}, {1, 1}];
f[u_, v_] := Through@{x, y, z}[u, v];
g[u_, v_] := {u, v, Log@Cos@u - Log@Cos@v};

meanCurvature[f_?VectorQ, {u_, v_}] := 
  Simplify[(Det[{D[f, {u, 2}], D[f, u], D[f, v]}] D[f, v].D[f, v] - 
      2 Det[{D[f, u, v], D[f, u], D[f, v]}] D[f, u].D[f, v] + 
      Det[{D[f, {v, 2}], D[f, u], D[f, v]}] D[f, u].D[f, 
         u])/(2 PowerExpand[
       Simplify[(D[f, u].D[f, u]*
            D[f, v].D[f, v] - (D[f, u].D[f, v])^2)]^(3/2)])];
eq = meanCurvature[f[u, v], {u, v}] == 0;
bcs = Flatten@{Thread[f[-1, v] == g[-1, v]], Thread[f[1, v] == g[1, v]],
   Thread[f[u, -1] == g[u, -1]], Thread[f[u, 1] == g[u, 1]]};

NDSolve[{eq}~Join~bcs, f[u, v], {u, v} ∈ region]

Of course, this doesn't work, because

NDSolve::underdet: There are more dependent variables, {x[u, v], y[u, v], z[u, v]}, than equations, so the system is underdetermined.

The problem is that we can "slide around" the parametrization along the surface and it doesn't change the geometry. Formally, for any smooth bijection $\phi$ from $U$ to itself, $f$ and $f\circ\phi$ represent the same surface. Even if I introduce additional conditions to fix a unique solution (which I don't know how to do), I expect I'll just end up with

NDSolve::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve.

Is there a better way to do this?

There are two related questions already on this site. "4 circular arcs, how plot the minimal surface?" is a special case with no posted answer. In "How can I create a minimal surface with trefoil knot as inner edge and circle as outer edge?", the desired minimal surface is not a topological disk (i.e. not simply connected), but using rotational symmetry one can divide it into six identical simply-connected pieces.

Other resources dealing with minimal surfaces in Mathematica are O. Michael Melko's article "Visualizing Minimal Surfaces" and the Mathematica code provided by the Minimal Surface Archive, but at first glance they both seem to be concerned with visualizing and manipulating minimal surfaces whose parametrization is already known.

  • 1
    $\begingroup$ Very nice question. Thanks for taking the time ask it and write this up so carefully. +1 $\endgroup$
    – halirutan
    Commented Jan 21, 2015 at 17:39
  • 1
    $\begingroup$ "we can slide around the parametrization" It is diffeomorphism :-) $\endgroup$
    – ybeltukov
    Commented Jan 22, 2015 at 2:18
  • $\begingroup$ My two cents.The catenoid/helicoid isometric morphing preserves Gauss curvature K and zero mean curvature H as a special minimal surface case. As one workaround instead of starting with given closed boundary of disc and attempting to find the minimal surface spanned in it in direct computation, it would be perhaps insightful to take arbitrary closed loops written on the catenoid, find minimal surface using a Mathematica FEM algorithm and to directly verify with known solution. contd $\endgroup$
    – Narasimham
    Commented Jan 22, 2015 at 19:31
  • $\begingroup$ This way it allows one to discover or formulate a certain relationship among differentials with functional relationships that allow generalization to advantage into other cases also. Using physical soap films spanning inside loops and using holographic optical methods is another easy experimental verification method. $\endgroup$
    – Narasimham
    Commented Jan 22, 2015 at 19:38
  • $\begingroup$ @Narasimham: Sounds like a good idea. And with ybeltukov's code, you can now try it! Let me know if you find anything interesting :) $\endgroup$
    – user484
    Commented Jan 22, 2015 at 19:56

4 Answers 4


Here is a method that utilizes $H^1$-gradient flows. This is far quicker than the $L^2$-gradient flow (a.k.a. mean curvature flow) or using FindMinimum and friends, in particular when dealing with finely discretized surfaces.


For those who are interested: A major reason for numerical slowness of $L^2$-gradient flow is the Courant–Friedrichs Lewy condition, which enforces the time step size in explicit integration schemes for parabolic PDEs to be proportial to the maximal cell diameter of the mesh. This leads to the need for many time iterations for fine meshes. Another problem is that the Hessian of the surface area with respect to the surface positions is highly ill-conditioned (both in the continuous as well as in the discretized setting.)

In order to compute $H^1$-gradients, we need the Laplace-Beltrami operator of an immersed surface $\varSigma$, or rather its associated bilinear form

$$ a_\varSigma(u,v) = \int_\varSigma \langle\operatorname{d} u, \operatorname{d} v \rangle \, \operatorname{vol}, \quad u,\,v\in H^1(\varSigma;\mathbb{R}^3).$$

The $H^1$-gradient $\operatorname{grad}^{H^1}_\varSigma(F) \in H^1_0(\varSigma;\mathbb{R}^3)$ of the area functional $F(\varSigma)$ solves the Poisson problem

$$a_\varSigma(\operatorname{grad}^{H^1}_\varSigma(F),v) = DF(\varSigma) \, v \quad \text{for all $v\in H^1_0(\varSigma;\mathbb{R}^3)$}.$$

When the gradient at the surface configuration $\varSigma$ is known, we simply translate $\varSigma$ by $- \delta t \, \operatorname{grad}^{H^1}_\varSigma(F)$ with some step size $\delta t>0$. By the way, this leads to the same algorithm as in Pinkal, Polthier - Computing discrete minimal surfaces and their conjugates (although the authors motivate the method in quite a different way). Surprisingly, the Fréchet derivative $DF(\varSigma)$ is given by

$$ DF(\varSigma) \, v = \int_\varSigma \langle\operatorname{d} \operatorname{id}_\varSigma, \operatorname{d} v \rangle \, \operatorname{vol},$$

so, we can also use the discretized Laplace-Beltrami operator to compute it.


Unfortunately, Mathematica cannot deal with finite elements on surfaces (yet). Therefore, I provide some code to assemble the Laplace-Beltrami operator of a triangular mesh.

   getLaplacian = Quiet[Block[{xx, x, PP, P, UU, U, f, Df, u, Du, g, integrand, quadraturepoints, quadratureweights}, 
    xx = Table[Part[x, i], {i, 1, 2}];
    PP = Table[Compile`GetElement[P, i, j], {i, 1, 3}, {j, 1, 3}];
    UU = Table[Part[U, i], {i, 1, 3}];

    (*local affine parameterization of the surface with respect to the "standard triangle"*)
    f[x_] = PP[[1]] + x[[1]] (PP[[2]] - PP[[1]]) + x[[2]] (PP[[3]] - PP[[1]]);
    Df[x_] = D[f[xx], {xx}];
    (*the Riemannian pullback metric with respect to f*)        
    g[x_] = Df[xx]\[Transpose].Df[xx];
    (*an affine function u and its derivative*)
    u[x_] = UU[[1]] + x[[1]] (UU[[2]] - UU[[1]]) + x[[2]] (UU[[3]] - UU[[1]]);
    Du[x_] = D[u[xx], {xx}];
    integrand[x_] = 1/2 D[Du[xx].Inverse[g[xx]].Du[xx] Sqrt[Abs[Det[g[xx]]]], {{UU}, 2}];
    (*since the integrand is constant over each triangle, we use a one- point Gauss quadrature rule (for the standard triangle)*)        
    quadraturepoints = {{1/3, 1/3}};
    quadratureweights = {1/2};
    With[{code = N[quadratureweights.Map[integrand, quadraturepoints]]},
     Compile[{{P, _Real, 2}},
      CompilationTarget -> "C",
      RuntimeAttributes -> {Listable},
      Parallelization -> True,
      RuntimeOptions -> "Speed"

getLaplacianCombinatorics = Quiet[Module[{ff},
      code = Flatten[Table[Table[{Compile`GetElement[ff, i], Compile`GetElement[ff, j]}, {i, 1, 3}], {j, 1, 3}], 1]
     Compile[{{ff, _Integer, 1}},
      CompilationTarget -> "C",
      RuntimeAttributes -> {Listable},
      Parallelization -> True,
      RuntimeOptions -> "Speed"

LaplaceBeltrami[pts_, flist_, pat_] := With[{
    spopt = SystemOptions["SparseArrayOptions"],
    vals = Flatten[getLaplacian[Partition[pts[[flist]], 3]]]
    SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
    SparseArray[Rule[pat, vals], {Length[pts], Length[pts]}, 0.],

Now we can minimize: We utilize that the differential of area with respect to vertex positions pts equals LaplaceBeltrami[pts, flist, pat].pts. I use constant step size dt = 1; this works surprisingly well. Of course, one may add a line search method of one's choice.

areaGradientDescent[R_MeshRegion, stepsize_: 1., steps_: 10, 
   reassemble_: False] := 
  Module[{method, faces, bndedges, bndvertices, pts, intvertices, pat,
     flist, A, S, solver}, Print["Initial area = ", Area[R]];
   method = If[reassemble, "Pardiso", "Multifrontal"];
   pts = MeshCoordinates[R];
   faces = MeshCells[R, 2, "Multicells" -> True][[1, 1]];
   bndedges = Developer`ToPackedArray[Region`InternalBoundaryEdges[R][[All, 1]]];
   bndvertices = Union @@ bndedges;
   intvertices = Complement[Range[Length[pts]], bndvertices];
   pat = Flatten[getLaplacianCombinatorics[faces], 1];
   flist = Flatten[faces];
   Do[A = LaplaceBeltrami[pts, flist, pat];
    If[reassemble || i == 1, 
     solver = LinearSolve[A[[intvertices, intvertices]], Method -> method]];
    pts[[intvertices]] -= stepsize solver[(A.pts)[[intvertices]]];, {i, 1, steps}];
   S = MeshRegion[pts, MeshCells[R, 2], PlotTheme -> "LargeMesh"];
   Print["Final area = ", Area[S]];

Example 1

We have to create some geometry. Any MeshRegion with triangular faces and nonempty boundary will do (although it is not guaranteed that an area minimizer exists).

h = 0.9;
R = DiscretizeRegion[
  ImplicitRegion[{x^2 + y^2 + z^2 == 1}, {{x, -h, h}, {y, -h, h}, {z, -h, h}}],
  MaxCellMeasure -> 0.00001, 
  PlotTheme -> "LargeMesh"

enter image description here

And this is all we have to do for minimization:

areaGradientDescent[R, 1., 20., False]

Initial area = 8.79696

Final area = 7.59329

enter image description here

Example 2

Since creating interesting boundary data along with suitable initial surfaces is a bit involved and since I cannot upload MeshRegions here, I decided to compress the initial surface for this example into these two images:

enter image description here

enter image description here

The surface can now be obtained with

R = MeshRegion[
  Polygon[Round[#/Min[#]] &@ Transpose[ ImageData[Import["https://i.sstatic.net/WfjOL.png"]]]]

enter image description here

With the function LoopSubdivide from this post, we can successively refine and minimize with

SList = NestList[areaGradientDescent@*LoopSubdivide, R, 4]

enter image description here

Here is the final minimizer in more detail:

enter image description here

Final Remarks

If huge deformations are expected during the gradient flow, it helps a lot to set reassemble = True. This uses always the Laplacian of the current surface for the gradient computation. However, this is considerably slower since the Laplacian has to be refactorized in order to solve the linear equations for the gradient. Using "Pardiso" as Method helps a bit.

Of course, the best we can hope to obtain this way is a local minimizer.


The package "PardisoLink`" makes reassembly more comfortable. This is possible due to the fact that the Pardiso solver can reuse its symbolic factorization and because I included the contents of this post into the package. Here is the new optimization routine that can be used as alternative to areaGradientDescent above.



Options[areaGradientDescent2] = {
   "StepSize" -> 1.,
   "MaxIterations" -> 20,
   "Tolerance" -> 10^-6,
   "Reassemble" -> True

areaGradientDescent2[R_MeshRegion, OptionsPattern[]] := 
  Module[{faces, flist, bndedges, bndvertices, pts, intvertices, pat, 
    A, S, solver, assembler, TOL, maxiter, reassemble, stepsize, b, u, res, iter
    }, Print["Initial area = ", Area[R]];
   TOL = OptionValue["Tolerance"];
   maxiter = OptionValue["MaxIterations"];
   reassemble = OptionValue["Reassemble"];
   stepsize = OptionValue["StepSize"];
   pts = MeshCoordinates[R];
   faces = MeshCells[R, 2, "Multicells" -> True][[1, 1]];
   bndedges = 
    Developer`ToPackedArray[Region`InternalBoundaryEdges[R][[All, 1]]];
   bndvertices = Union @@ bndedges;
   intvertices = Complement[Range[Length[pts]], bndvertices];
   pat = Flatten[getLaplacianCombinatorics[faces], 1];
   flist = Flatten[faces];
   faces =.;
   assembler = Assembler[pat, {Length[pts], Length[pts]}];
   A = assembler[getLaplacian[Partition[pts[[flist]], 3]]];
   solver = 
    Pardiso[A[[intvertices, intvertices]], "MatrixType" -> 2];
   b = (A.pts)[[intvertices]];
   u = solver[b];
   res = Sqrt[Flatten[u].Flatten[b]];
   iter = 0;
   While[res > TOL && iter < maxiter,
    pts[[intvertices]] -= stepsize u;
    A = assembler[getLaplacian[Partition[pts[[flist]], 3]]];
    b = (A.pts)[[intvertices]];
     Quiet[solver["Update"[A[[intvertices, intvertices]]]]]];
    u = solver[b];
    res = Sqrt[Flatten[u].Flatten[b]];

   S = MeshRegion[pts, MeshCells[R, 2], PlotTheme -> "LargeMesh"];
   Print["Final area = ", Area[S], ". Iterations used = ", iter, 

Applied to example 1 above, it can be used as follows (note the different argument pattern).

h = 0.9;
R = DiscretizeRegion[
   ImplicitRegion[{x^2 + y^2 + z^2 == 1}, {{x, -h, h}, {y, -h, 
      h}, {z, -h, h}}], MaxCellMeasure -> 0.00001, 
   PlotTheme -> "LargeMesh"];

areaGradientDescent[R, 1., 10, False]; // AbsoluteTiming // First
areaGradientDescent[R, 1., 10, True]; // AbsoluteTiming // First
areaGradientDescent2[R, "MaxIterations" -> 10, "Reassemble" -> True]; // AbsoluteTiming // First




Even with reassembly activated, it is now faster than areaGradientDescent with deactivated assembly and more then twice as fast as areaGradientDescent with activated assembly.

  • $\begingroup$ Is it valid to use a degenerate bilinear form in gradient descent? After all, the Laplace-Beltrami operator has a nontrivial null space containing all constant functions, I believe. Or does the restriction to interior vertices take care of that? $\endgroup$
    – user484
    Commented May 19, 2018 at 4:52
  • $\begingroup$ The bilinear form is positive definite (and coercive w.r.t the $H^1$-norm) on the space of functions that vanish at the boundary. Or in terms of the Laplacian: Each function in its null space is locally constant; since we are dealing here with surfaces whose connected components need to have nontrivial boundary, only the zero-function lies in the null space. So you are right: the restriction to the interior vertices takes care of that. $\endgroup$ Commented May 19, 2018 at 7:53
  • $\begingroup$ Rahul, I feel honored. $\endgroup$ Commented May 19, 2018 at 12:11

Edit: added Gradient -> grad[vars] option. Without this small option the code was several orders of magnitude slower.

Yes, it can! Unfortunately, not automatically.

There are different algorithms to do it (see special literature, e.g. Dziuk, Gerhard, and John E. Hutchinson. A finite element method for the computation of parametric minimal surfaces. Equadiff 8, 49 (1994) [pdf] and references therein). However I'm going to implement the simplest method as possible. Just split a trial initial surface to triangles and minimize the total area of triangles.

boundary = HoldPattern[{_, _, z_} /; Abs[z] > 0.0001 && Abs[z - 1] > 0.0001];
g = ParametricPlot3D[{Cos[u] (1 + 0.3 Sin[5 u + π v]), 
   Sin[u] (1 + 0.3 Sin[5 u + π v]), v}, {u, 0, 2 π}, {v, 0, 
   1}, PlotPoints -> {100, 15}, MaxRecursion -> 0, Mesh -> None, 
  NormalsFunction -> None]

enter image description here

It is far from ideal. Let's convert it to MeshRegion.

R = DiscretizeGraphics@Normal@g;
vc = MeshCoordinates@R;
cells = MeshCells[R, 2];
{t0, t1, t2} = Transpose@cells[[All, 1]];
pts = Flatten@Position[vc, boundary];
P = SparseArray[Transpose@{Join[t0, t1, t2], Range[3 Length@t0]} -> 
    ConstantArray[1, 3 Length@t0]];
Ppts = P[[pts]];

Here P is an auxiliary matrix which converts a triangle number to a vertex number. pts is a list of numbers of vertices which did't lie on boundaries (the current implementation contains explicit conditions; it is ugly, but I don't know how to do it better).

The total area and the corresponding gradient

area[v_List] := Module[{vc = vc, u1, u2},
   vc[[pts]] = v;
   u1 = vc[[t1]] - vc[[t0]];
   u2 = vc[[t2]] - vc[[t0]];
   Total@Sqrt[(u1[[All, 1]] u2[[All, 2]] - u1[[All, 2]] u2[[All, 1]])^2 +
       (u1[[All, 2]] u2[[All, 3]] - u1[[All, 3]] u2[[All, 2]])^2 +
       (u1[[All, 3]] u2[[All, 1]] - u1[[All, 1]] u2[[All, 3]])^2]/2];

grad[v_List] := Flatten@Module[{vc = vc, u1, u2, a, g1, g2},
    vc[[pts]] = v;
    u1 = vc[[t1]] - vc[[t0]];
    u2 = vc[[t2]] - vc[[t0]];
    a = Sqrt[(u1[[All, 1]] u2[[All, 2]] - u1[[All, 2]] u2[[All, 1]])^2 +
        (u1[[All, 2]] u2[[All, 3]] - u1[[All, 3]] u2[[All, 2]])^2 +
        (u1[[All, 3]] u2[[All, 1]] - u1[[All, 1]] u2[[All, 3]])^2]/2;
    g1 = (u1 Total[u2^2, {2}] - u2 Total[u1 u2, {2}])/(4 a);
    g2 = (u2 Total[u1^2, {2}] - u1 Total[u1 u2, {2}])/(4 a);
    Ppts.Join[-g1 - g2, g1, g2]];

In other words, grad is finite-difference form of the mean curvature flow. Such exact calculation of grad considerably increases the speed of the evaluation.

vars = Table[Unique[], {Length@pts}];
v = vc;
v[[pts]] = First@FindArgMin[area[vars], {vars, vc[[pts]]}, Gradient -> grad[vars],
     MaxIterations -> 10000, Method -> "ConjugateGradient"];

Graphics3D[{EdgeForm[None], GraphicsComplex[v, cells]}]

enter image description here

The result is fine! However the visualization will be better with VertexNormal option and different colors

normals[v_List] := Module[{u1, u2},
  u1 = v[[t1]] - v[[t0]];
  u2 = v[[t2]] - v[[t0]];
  P.Join[#, #, #] &@
   Transpose@{u1[[All, 2]] u2[[All, 3]] - u1[[All, 3]] u2[[All, 2]],
     u1[[All, 3]] u2[[All, 1]] - u1[[All, 1]] u2[[All, 3]],
     u1[[All, 1]] u2[[All, 2]] - u1[[All, 2]] u2[[All, 1]]}]

Graphics3D[{EdgeForm[None], FaceForm[Red, Blue], 
  GraphicsComplex[v, cells, VertexNormals -> normals[v]]}]

enter image description here

Costa Minimal Surface

Let's try something interesting, e.g. Costa-like minimal surface. The main problem is the initial surface with a proper topology. We can do it with "knife and glue".

Pieces of surfaces (central connector, middle sheet, top&bottom sheet):

r1 = 10.;
r2 = 6.;
h = 5.0;
n = 60;
m = 50;
hole0 = Table[{Cos@φ, Sin@φ} (2 - Abs@Sin[2 φ]), {φ, 2 π/n, 2 π, 2 π/n}];
hole1 = Table[{Cos@φ, Sin@φ} (2 + Abs@Sin[2 φ]), {φ, 2 π/n, 2 π, 2 π/n}];
hole2 = Table[{Cos@φ, Sin@φ} (2 + Sin[2 φ]), {φ, 2 π/n, 2 π, 2 π/n}];
circle = Table[{Cos@φ, Sin@φ}, {φ, 2 π/m, 2 π, 2 π/m}];
bm0 = ToBoundaryMesh["Coordinates" -> hole0, 
   "BoundaryElements" -> {LineElement@Partition[Range@n, 2, 1, 1]}];
{bm1, bm2} = ToBoundaryMesh["Coordinates" -> Join[#, #2 circle], 
     "BoundaryElements" -> {LineElement@
        Join[Partition[Range@n, 2, 1, 1], 
         n + Partition[Range@m, 2, 1, 1]]}] & @@@ {{hole1, 
     r1}, {hole2, r2}};
{em0, em1, em2} = ToElementMesh[#, "SteinerPoints" -> False, "MeshOrder" -> 1, 
     "RegionHoles" -> #2, MaxCellMeasure -> 0.4] & @@@ {{bm0, 
     None}, {bm1, {{0, 0}}}, {bm2, {0, 0}}};
MeshRegion /@ {em0, em1, em2}

enter image description here

The option "SteinerPoints" -> False holds boundary points for further gluing. The option "MeshOrder" -> 1 forbids unnecessary second-order mid-side nodes. A final glued surface

boundary = HoldPattern[{x_, y_, z_} /; 
    Not[x^2 + y^2 == r1^2 && z == 0 || x^2 + y^2 == r2^2 && Abs@z == h]];
g = Graphics3D[{FaceForm[Red, Blue], 
   GraphicsComplex[em0["Coordinates"] /. {x_, y_} :> {-x, y, 0.}, 
    Polygon@em0["MeshElements"][[1, 1]]], 
   GraphicsComplex[em1["Coordinates"] /. {x_, y_} :> {x, y, 0}, 
    Polygon@em1["MeshElements"][[1, 1]]], 
   GraphicsComplex[em2["Coordinates"] /. {x_, y_} :> {-x, y, 
       h Sqrt@Rescale[Sqrt[
          x^2 + y^2], {2 + (2 x y)/(x^2 + y^2), r2}]}, 
    Polygon@em2["MeshElements"][[1, 1]]], 
   GraphicsComplex[em2["Coordinates"] /. {x_, y_} :> {y, 
       x, -h Sqrt@Rescale[Sqrt[x^2 + y^2], {2 + (2 x y)/(x^2 + y^2), r2}]}, 
    Polygon@em2["MeshElements"][[1, 1]]]}]

enter image description here

After the minimization code above we get

enter image description here

  • 1
    $\begingroup$ In your answer you have an initial guess for the surface which is the parametric plot. When I understand the question correctly, then the starting point is only a curve. Do you know a solution when we have for instance this curve, let's assume not even analytically but as coordinate list. Do you have an idea for this? Additionally, it would be awesome if you could provide some links to algorithms/literature. $\endgroup$
    – halirutan
    Commented Jan 22, 2015 at 2:54
  • 1
    $\begingroup$ Fantastic! I for one am perfectly happy to provide an initial surface. @halirutan: For your curve one could simply form a "cone" by connecting all the points to the origin. I think that works for arbitrary curves, but I don't know if self-intersections will cause the result to get stuck in local minima. $\endgroup$
    – user484
    Commented Jan 22, 2015 at 4:53
  • $\begingroup$ @halirutan: This surface is like the monkey saddle where there are now four humps/falls instead of three.The surface is Re or Im part of (x + I y)^4. $\endgroup$
    – Narasimham
    Commented Jan 23, 2015 at 0:58
  • 1
    $\begingroup$ With respect to generating normals: have you already seen this? $\endgroup$ Commented May 26, 2015 at 1:35
  • 1
    $\begingroup$ I know, that's why I linked you to it; it seems Max's weighted average normal might help here. $\endgroup$ Commented May 26, 2015 at 23:32

I've wrapped up @ybeltukov's code into a function that works for an arbitrary MeshRegion surface.

First we need to find the boundary vertices, which will remain fixed. If the MeshRegion represents a 2-dimensional manifold with boundary, then every internal vertex has as many edges as it has faces, but every boundary vertex has one extra edge.

boundaryVertices[mr_] := Module[{edges, faces},
  edges = First /@ MeshCells[mr, 1];
  faces = First /@ MeshCells[mr, 2];
    Merge[{Counts[Flatten@edges], Counts[Flatten@faces]}, 
     Greater @@ # &], TrueQ]]

Then computing the minimal surface is a near-verbatim copy of @ybeltukov's code:

findMinimalSurface[mr_] := 
 Module[{vc, cells, t0, t1, t2, bc, pts, P, area, grad, vars, v},
  vc = MeshCoordinates@mr;
  cells = MeshCells[mr, 2];
  {t0, t1, t2} = Transpose@cells[[All, 1]];
  pts = Complement[Range[Length@vc], boundaryVertices[mr]];
  P = SparseArray[
     Transpose@{Join[t0, t1, t2], Range[3 Length@t0]} -> 
      ConstantArray[1, 3 Length@t0]][[pts]];
  area[v_List] := Module[{vc = vc, u1, u2}, vc[[pts]] = v;
    u1 = vc[[t1]] - vc[[t0]];
    u2 = vc[[t2]] - vc[[t0]];
      Sqrt[(u1[[All, 1]] u2[[All, 2]] - 
           u1[[All, 2]] u2[[All, 1]])^2 + (u1[[All, 2]] u2[[All, 3]] -
            u1[[All, 3]] u2[[All, 2]])^2 + (u1[[All, 3]] u2[[All, 
             1]] - u1[[All, 1]] u2[[All, 3]])^2]/2];
  grad[v_List] := 
   Flatten@Module[{vc = vc, u1, u2, a, g1, g2}, vc[[pts]] = v;
     u1 = vc[[t1]] - vc[[t0]];
     u2 = vc[[t2]] - vc[[t0]];
     a = Sqrt[(u1[[All, 1]] u2[[All, 2]] - 
            u1[[All, 2]] u2[[All, 1]])^2 + (u1[[All, 2]] u2[[All, 
              3]] - 
            u1[[All, 3]] u2[[All, 2]])^2 + (u1[[All, 3]] u2[[All, 
              1]] - u1[[All, 1]] u2[[All, 3]])^2]/2;
     g1 = (u1 Total[u2^2, {2}] - u2 Total[u1 u2, {2}])/(4 a);
     g2 = (u2 Total[u1^2, {2}] - u1 Total[u1 u2, {2}])/(4 a);
     P.Join[-g1 - g2, g1, g2]];
  vars = Table[Unique[], {Length@pts}];
  v = vc;
  v[[pts]] = 
   First@FindArgMin[area[vars], {vars, vc[[pts]]}, 
     Gradient -> grad[vars], MaxIterations -> 1000];
  MeshRegion[v, cells]]

If all you have is the boundary curve, you can create an initial surface as a "cone" that connects every point on the curve to the origin. I do this by defining a Disk-shaped region and moving its vertices to lie on the cone. I prefer ToElementMesh because it lets you choose a finer resolution at the boundary than in the interior, allowing us to closely follow the input curve without wasting too many triangles on the smooth interior of the surface.

createInitialSurface[g_, {t_, t0_, t1_}] := 
 With[{mesh = 
     ToElementMesh[Disk[], "MeshOrder" -> 1, MaxCellMeasure -> 0.01, 
      "MaxBoundaryCellMeasure" -> 0.05]}, 
   With[{r = Norm@#, θ = ArcTan @@ #}, 
      r (g /. t -> Rescale[θ, {-π, π}, {t0, t1}])] & /@
     MeshCoordinates[mesh], MeshCells[mesh, 2]]]

By the way, here's some code to draw a prettier rendering of a MeshRegion surface.

showRegion[mr_] := 
   GraphicsComplex[MeshCoordinates[mr], MeshCells[mr, 2]]}]

Now we can solve @halirutan's example:

s = createInitialSurface[{Cos[t], Sin[t], Cos[4 t]/2}, {t, 0, 2 π}];
m = findMinimalSurface[s];

enter image description here

enter image description here

It's similar to the plot of $z=\operatorname{Re}\bigl((x+iy)^4\bigr)$, but if you draw both surfaces together you find that they're not identical.

We can also solve the previous question, "4 circular arcs, how plot the minimal surface?":

g[t_] := Piecewise[{{{1, -Cos@t, Sin@t}, 0 <= t <= π},
                    {{-Cos@t, 1, Sin@t}, π <= t <= 2 π},
                    {{-1, Cos@t, Sin@t}, 2 π <= t <= 3 π},
                    {{Cos@t, -1, Sin@t}, 3 π <= t <= 4 π}}];
ParametricPlot3D[g[t], {t, 0, 4 π}]

enter image description here

showRegion@findMinimalSurface@createInitialSurface[g[t], {t, 0, 4 π}]

enter image description here

There are a few magic numbers in the code that you can change to adjust the quality of the results. In findMinimalSurface, there's MaxIterations -> 1000 (which I reduced from @ybeltukov's 10000 because I didn't want to wait that long). You could also try a different Method such as "ConjugateGradient" in the same FindArgMax call. In createInitialSurface, there's MaxCellMeasure -> 0.01 and "MaxBoundaryCellMeasure" -> 0.05 which I picked to look OK and not be too slow on the presented examples. Also, if your boundary curve is only piecewise smooth, such as the square-shaped example I gave in the question, you may want to replace the Disk by a different 2D region that is closer to the shape of the expected surface.

  • $\begingroup$ I am giving code for closed contour on a cone with $ \alpha = \pi/6 $ tangential radius 0.4 $\endgroup$
    – Narasimham
    Commented Mar 14, 2015 at 9:40
  • $\begingroup$ Fine, it works so well. One or two things: 1) Can they be exported to be used as input for structural analysis? 2) Can you use it to find the film surface for a half of a Viviani Curve?(as the coordinates are readily available). 3) If you care to, I am also giving Mma program for a " Ring on Cone" a circular patch wrapped on a cone? It would be like the connected 4 semi-circle arcs but without curvature discontinuities at connecting points , for the boundary. $\endgroup$
    – Narasimham
    Commented Mar 14, 2015 at 9:42
  • $\begingroup$ It didn't work so well for computing surfaces bounded by Enneper's wire. (These are not Enneper's surface.) Try the following: <code> Enneper[R_] := ( en = createInitialSurface[{R Cos[t] - (1/3) R^3 Cos[3 t], -R Sin[t] - (1/3)R^3 Sin[3 t], R^2 Cos[2t]}, {t, 0,2 Pi}]; </code> enneper2 = findMinimalSurface[en] ) Manipulate[showRegion[Enneper[R]],{R,1,1.7,0.1}] $\endgroup$ Commented Mar 3, 2016 at 0:22
  • $\begingroup$ Sorry, I didn't know "comments can only be edited for five minutes", and I didn't manage to get the code displayed right in that time. $\endgroup$ Commented Mar 3, 2016 at 0:32
  • $\begingroup$ Also, I would like to draw the Gauss map of these surfaces. There is mention on the web of a GaussMap package, but I can't find it. Is it available? There's also a "vertexNormals" option to display a MeshRegion, but you don't see them when you ask Information[myMeshRegion]. $\endgroup$ Commented Mar 3, 2016 at 0:55

No answer here but only further forward suggestions with my thoughts on the topic.

We can start with any contour C but more conveniently consider a loop with known closed form parametrization. Supposing we start with an "ellipse" contour C written on a unit sphere ( defined by achille hui in SE Math in reply to my question or any Monkey saddle variant) with constant mean curvature H, implement the code and make this constant $ H = 0 $ in Manipulation of the minimal surface spanning across the ellipse.

For implementation it is a physical realization of soap-film on a cut-out contour C on the sphere where we apply pressure $ p/N = 2 H $ physically from inside the spherical container.

Equilibrium of forces equation across the minimal film is taken from membrane structural mechanics theory:

$$ \frac{N_1}{R_1} + \frac{N_2}{R_2} = p $$

where $N$ is surface tension or force per unit length, $R$ are principal radii of normal curvature, not along asymptotes. It reduces by notation

$$ k_1 = 1/ R_1, k_2 = 1/ R_2 ,N_1 = N_2 = N $$


$$ k_1 + k_2 = p/N = T = 2H, const. $$

which ODE describes const. H surface reducible to $ H=0$ minimal surfaces.

( For H = constant surfaces of revolution we have Delaunay Unduloids, and for $ p = 0, H = 0 $, the catenoid r = c \cosh z/c)

Integration is performed for surface on an orthogonal curvilinear net of asymptotic lines which is associated with all minimal films/patches. We can imagine a net placed on a soldier's helmet, but where there is a hole on the helmet.

So all we need to implement in the algorithm is only on the boundary interfacing with an asymptotic net.

Assume an initial $ \alpha_{start} $ at an arbitrary starting point. Go around the loop and find $ \alpha_{end} $. Iterate in a DO loop around the contour until satisfaction of accuracy of $\alpha $ error/difference.


Sine-Gordon Equation exists for constant K < 0 asymptotic lines in a more general situation, but for minimal surfaces no ode or pde is availble afaik. Its formulation appears to be an immediately needed good area for surface theory math research.

A bonanza for future is... adopting the same program numerically later on to find a fishnet, as an asymptotic net with constant $ K< 0 $. There is none available for non-surfaces of revolution as of now, attractive for Mathematica FEM.


Consideration of asymptotic lines of an orthogonal asymptotic net I guess would be very rewarding.


For the orthogonal net, asymptotic lines ($k_n=0$) are placed at $ \pm \pi/4 $ to the principal curvature directions.

Simple but significant curvature relations are shown on Mohr's circle of curvature tensors:

$$ k_n = \tau_g = \pm\; c $$

These arise from second fundamental form $$ L=0, N=0, M \neq 0, $$

$$ K = -(M/H)^2 , \tau_g = \pm (M/H), H = - FM/H^2. $$

For consideration of a circular patch isometrically draped on a cone:

(* Constant Geodesic Curvature on Cone semi vert angle \[Alpha]  " \
ChapathiChalupa.nb " *)
ri = 0.6 ; Rg = 0.4; smax = 3; \[Alpha] = Pi/6; zi = ri Cot[\[Alpha]] ;
Chapathi = {SI'[s] == 1/Rg - Sin[PH[s]] Sin[SI[s]]/R[s],  
   SI[0] == Pi/2, PH'[s] == 0., PH[0] == \[Alpha], 
   R'[s] == Sin[PH[s]] Cos[SI[s]], Z'[s] == Cos[PH[s]] Cos[SI[s]], 
   TH'[s] == Sin[SI[s]]/R[s], R[0] == ri, TH[0] == 0, Z[0] == zi};
NDSolve[Chapathi, {SI, PH, R, TH, Z}, {s, 0, smax}];
{si[t_], ph[t_], r[t_], th[t_], 
   z[t_]} = {SI[t], PH[t], R[t], TH[t], Z[t]} /. First[%];
Plot[{r[s], th[z], z[s]}, {s, 0, smax}];
ParametricPlot[{{z[s], r[s]}, {z[s], -r[s]}, {z[s], 
   r[s] Cos[th[s]]}, {z[s], r[s] Sin[th[s]]}}, {s, .0, smax}, 
 PlotLabel -> PROJNS, GridLines -> Automatic, 
 AspectRatio -> Automatic, PlotStyle -> {Red, Thick}]
Cone = ParametricPlot3D[{r[s] Cos[th[s] + v], r[s] Sin[th[s] + v], 
    z[s]}, {s, .0, smax}, {v, 0, 2 Pi}, PlotLabel -> ChapathiChalupa, 
   PlotStyle -> {Yellow, Opacity[0.05]}, Boxed -> False, Axes -> None];
Boundary = 
  ParametricPlot3D[{r[s] Cos[th[s]], r[s] Sin[th[s]], z[s]}, {s, .0, 
    smax}, PlotLabel -> 3 _D Projn, PlotStyle -> {Red, Thick}];
Show[Cone, Boundary]
Table[ {Z, z[s], R, r[s], ThDeg, th[s] 180/Pi, s} , { s, 0, smax, 
   smax/20.}] // TableForm



Posting after a long pause! This is not a solution , a hyperbolic geodesic boundary suggestion I defined this earlier, that can be used now to demarcate a simple boundary on the familiar catenoid (of minimum radius 1) below. The boundary bifurcates area along 3-axes symmetry. Hope it could be useful as a standard model to accompany any new computation minimal surface spanned with triangulation and meshing using Mathematica. If found useful for our FEA work here shall give its parameterization.



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.