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.
Background
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.
Implementation
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}},
code,
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
]
]
]
];
getLaplacianCombinatorics = Quiet[Module[{ff},
With[{
code = Flatten[Table[Table[{Compile`GetElement[ff, i], Compile`GetElement[ff, j]}, {i, 1, 3}], {j, 1, 3}], 1]
},
Compile[{{ff, _Integer, 1}},
code,
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
]
]]];
LaplaceBeltrami[pts_, flist_, pat_] := With[{
spopt = SystemOptions["SparseArrayOptions"],
vals = Flatten[getLaplacian[Partition[pts[[flist]], 3]]]
},
Internal`WithLocalSettings[
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
SparseArray[Rule[pat, vals], {Length[pts], Length[pts]}, 0.],
SetSystemOptions[spopt]]
];
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]];
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"
]
And this is all we have to do for minimization:
areaGradientDescent[R, 1., 20., False]
Initial area = 8.79696
Final area = 7.59329
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:
The surface can now be obtained with
R = MeshRegion[
Transpose[ImageData[Import["https://i.sstatic.net/aaJPM.png"]]],
Polygon[Round[#/Min[#]] &@ Transpose[ ImageData[Import["https://i.sstatic.net/WfjOL.png"]]]]
]
With the function LoopSubdivide from this post, we can successively refine and minimize with
SList = NestList[areaGradientDescent@*LoopSubdivide, R, 4]
Here is the final minimizer in more detail:
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.
Update
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.
Needs["PardisoLink`"];
ClearAll[areaGradientDescent2];
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,
++iter;
pts[[intvertices]] -= stepsize u;
A = assembler[getLaplacian[Partition[pts[[flist]], 3]]];
b = (A.pts)[[intvertices]];
If[reassemble,
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,
"."];
S
];
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
0.875273
1.82184
0.848296
Even with reassembly activated, it is now faster than areaGradientDescent with deactivated assembly and more then twice as fast as areaGradientDescent with activated assembly.
