Brakke’s samples: soap films and foams

Soap films spanning wires, a bubble, and foams in a periodic box, from the sample datafiles that come with Surface Evolver (pyse.examples). The helpers are in docs/brakke_helpers.py.

Each sample is run twice from its datafile: once with the recipe Brakke wrote into it (its gogo command, or the usage its comments describe), once with a default ev.relax() given as many refinements as the recipe does (levels). The table shows the energy each ends with, the residual (how far from equilibrium, see Relaxing a surface), the facet count, the time, and what relax() reports: converged, stable, and any health issues. A lower energy on the same mesh is not always the better answer: a coarse mesh can lower its energy in ways the smooth surface can’t, so where an exact value is known, the error against it is the measure. Times are from the documentation build and vary.

This gallery assumes Getting started and Relaxing a surface; its pages can be read in any order.

import numpy as np
import pyvista as pv
import pysurfaceevolver as pyse
import brakke_helpers as bh

pv.set_jupyter_backend("static")

cube.fe: a cube becoming a sphere

A unit cube of fixed volume 1, its surface area the energy: it relaxes toward a sphere (area \((36\pi)^{1/3} \approx 4.836\)). The recipe refines twice with gradient and Newton steps in between.

print(bh.recipe("cube.fe"))
runs = bh.compare("cube.fe")
bh.table(runs, reference=(36*np.pi)**(1/3))
bh.show(runs)

g 5; r; g 5; hessian; r; g 5; hessian;

	energy	residual	facets	seconds	converged	stable	notes
recipe: g 5; r; g 5; hessian; r; g 5; hessian;	4.854459815	5.6e-04	384	0.00
relax(levels=2)	4.852344416	2.8e-09	384	0.02	yes	yes

	energy	error	residual	facets	seconds	converged	stable	notes
recipe: g 5; r; g 5; hessian; r; g 5; hessian;	4.854459815	+3.82e-03	5.6e-04	384	0.00
relax(levels=2)	4.852344416	+3.38e-03	2.8e-09	384	0.02	yes	yes

Both meshes are polyhedra inscribed in the sphere, with a little more area than it; relax() takes its residual to about 1e-9 where the recipe stops sooner, at nearly the same energy.

catman.fe: the catenoid of the manual’s tutorial

A film between two rings of radius 1 at height \(\pm 0.55\); the stable catenoid has area 6.511674. The recipe’s comments say what happens on this coarse mesh: after 50 steps “at this point have nearly a saddle point”, then “triangulation twists around to lower energy”. Twisting lowers the area of the coarse mesh, not its error.

print(bh.recipe("catman.fe"))
runs = bh.compare("catman.fe")
bh.table(runs, reference=6.511674)
bh.show(runs)

g; u; r; g 50; g 200;

	energy	residual	facets	seconds	converged	stable	notes
recipe: g; u; r; g 50; g 200;	6.433703167	4.9e-02	96	0.01
relax(levels=1)	6.458665438	2.0e-10	96	0.01	yes	yes

	energy	error	residual	facets	seconds	converged	stable	notes
recipe: g; u; r; g 50; g 200;	6.433703167	-1.20e-02	4.9e-02	96	0.01
relax(levels=1)	6.458665438	-8.14e-03	2.0e-10	96	0.01	yes	yes

Both are below the smooth catenoid (the rings are coarse polygons); the twisted mesh of the recipe is further below. relax() stops at the untwisted equilibrium, and leaves sliding along the wires (where the twist happens) to its gradient rounds.

cat.fe: a catenoid past its limit

Here the rings are as far apart as a catenoid can span (RMAX set at the fold). The neck collapses, and the recipe deletes the tiny edges there (t .05) and pops the pinched vertex (o): the film becomes two flat discs. relax() can’t converge the film, so as a last resort it takes gradient steps without tidying, lets the neck collapse, deletes the collapsed edges and pops (collapsed, popped).

print(bh.recipe("cat.fe"))
runs = bh.compare("cat.fe")
R = 1.5088795                       # RMAX: the rings' radius
discs = 2*3*R**2                    # two 12-gons inscribed in the rings (after one refinement)
bh.table(runs, reference=discs)
bh.show(runs)

r; u; g 120; t .05; o; g 5;

	energy	residual	facets	seconds	converged	stable	notes
recipe: r; u; g 120; t .05; o; g 5;	13.71287913	5.7e-01	72	0.01
relax(levels=1), 1 popped, 9 collapsed edges	13.66030407	0.0e+00	78	0.05	yes	yes

	energy	error	residual	facets	seconds	converged	stable	notes
recipe: r; u; g 120; t .05; o; g 5;	13.71287913	+3.85e-03	5.7e-01	72	0.01
relax(levels=1), 1 popped, 9 collapsed edges	13.66030407	-4.34e-10	0.0e+00	78	0.05	yes	yes

The recipe stops after g 5 on the discs; relax() finishes them: two regular 12-gons, the polygons the rings’ vertices make after one refinement.

catbody.fe: a catenoid enclosing a fixed volume

The same catenoid shape, now closed by plane ends and holding a fixed volume, so it has a pressure. The file adjusts that volume to find zero pressure (the unstable catenoid) in a separate script; here only its gogo.

print(bh.recipe("catbody.fe"))
runs = bh.compare("catbody.fe")
bh.show(runs)

u; r; g 5; r; g 5; r; g 5; hessian; hessian

	energy	residual	facets	seconds	converged	stable	notes
recipe: u; r; g 5; r; g 5; r; g 5; hessian; hessian	6.991913875	5.8e-04	1536	0.01
relax(levels=3)	6.99440838	1.1e-10	1536	0.27	yes	yes

octa.fe: a film in an octahedral frame

A symmetric film from the edges of an octahedral wire frame to its centre: six edges there carry four films each. Plateau’s laws (soap films meet three at a time along an edge, at 120 degrees, and such edges four at a time at a point) say that is unstable: splitting the junction lowers the area. The recipe pops junctions (o), deletes small edges (t .1) and iterates, several times; relax() pops after converging, relaxes again, and deletes the films that shrink away.

print(bh.recipe("octa.fe"))
runs = bh.compare("octa.fe")
bh.show(runs)

o; U; g 20; t .1; g 10; o; g 20; u; V; t .1; o; g 12; o; u; V; u; V; g 10; t .1; o; V; g 10; t .1; u; V; g 10; t .1; g 10;

	energy	residual	facets	seconds	converged	stable	notes
recipe: o; U; g 20; t .1; g 10; o; g 20; u; V; t .1; o; g 12; o; u; V; u; V; g 10; t .1; o; V; g 10; t .1; u; V; g 10; t .1; g 10;	22.627417	1.5e-06	90	0.00
relax(levels=0), 11 popped, 4 collapsed edges	22.627417	6.0e-14	90	0.04	yes	yes

Both reach the same film, with a central tetrahedral point.

Kelvin against Weaire-Phelan: twointor.fe and phelanc.fe

Lord Kelvin asked in 1887 how to divide space into cells of equal volume with the least interface area, and proposed the cell of twointor.fe (two cells in a periodic box). In 1993 Weaire and Phelan found a structure with two kinds of cells (phelanc.fe, eight cells) with about 0.3% less area. Per cell, the area over the cell volume to the power 2/3 compares them independently of the box.

kelvin = bh.compare("twointor.fe")
wp = bh.compare("phelanc.fe")

def per_cell(name, levels):
    ev = pyse.Evolver(name)
    ev.relax(levels=levels)
    n = len(ev.bodies().volume)
    V = float(np.sum(ev.bodies().volume))
    return ev.total_energy/n/(V/n)**(2/3)

k, w = per_cell("twointor.fe", 2), per_cell("phelanc.fe", 2)
print(f"area per cell / volume^(2/3): Kelvin {k:.5f}, Weaire-Phelan {w:.5f}, "
      f"{1 - w/k:.2%} less")

	energy	residual	facets	seconds	converged	stable	notes
recipe: g 5; r; g 5; hessian; r; g 5; hessian; hessian;	3.343153075	3.8e-10	1152	0.03
relax(levels=2)	3.343152078	3.6e-11	1152	0.13	yes	yes

	energy	residual	facets	seconds	converged	stable	notes
recipe: g 5; r; g 5; hessian; r; g 5; hessian;	21.15396187	1.8e-05	4416	0.11
relax(levels=2)	21.15392109	3.9e-13	4416	0.58	yes	yes

area per cell / volume^(2/3): Kelvin 2.65346, Weaire-Phelan 2.64424, 0.35% less

bh.show(wp)

symtest.fe: the Kelvin foam in a quotient space

The same two Kelvin cells, but in Evolver’s quotient-space representation (a symmetry group instead of the periodic box). The recipe is twointor.fe’s.

runs = bh.compare("symtest.fe", bh.recipe("twointor.fe"))

	energy	residual	facets	seconds	converged	stable	notes
recipe: g 5; r; g 5; hessian; r; g 5; hessian; hessian;	3.343153075	3.8e-10	1152	0.01
relax(levels=2)	3.343152077	3.6e-11	1152	0.09	yes	no	unstable: 1 direction(s) lower the energy

The energy is that of twointor.fe, but relax() reports the surface unstable, with one negative Hessian eigenvalue (as Brakke’s Evolver 2.70a does), and twointor.fe stable. The instability is real: the volumes in this file’s bodies section are comments (// volume 0.500), so its two cells have no volume constraint, and one can grow at the expense of the other. Moving along that mode (Evolver’s saddle) lowers the energy, and relaxing from there shrinks one cell away.

ev = pyse.Evolver("symtest.fe")
print(ev.bodies().fixed)                    # no volume constraints
ev.relax(levels=1)
e0 = ev.total_energy
ev.saddle()                                 # a step along the negative mode
print(f"energy {e0:.6f} -> {ev.total_energy:.6f}")

[False False]
energy 3.344298 -> 3.330655

/opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/pysurfaceevolver/_relax.py:235: UnstableEquilibriumWarning: relax() converged to an unstable equilibrium: 1 direction(s) lower the energy (lowest eigenvalue -2.65)
  return self._relax(tol=tol, max_iter=max_iter, levels=levels, newton=newton,

Triply periodic minimal surfaces: pcell.fe, dcell.fe, IWP.fe

From Brakke’s collection of periodic minimal surfaces online (included with pySE): Schwarz’s P and D surfaces and Schoen’s I-WP, each a small piece of the surface in a quotient space (a symmetry group generates the rest).

for name in ("pcell.fe", "dcell.fe", "IWP.fe"):
    runs = bh.compare(name)

	energy	residual	facets	seconds	converged	stable	notes
recipe: g 5; r; g 5; r; g 5; r; g 10; u; refine edge where on_constraint 2; refine edge where on_constraint 4; hessian; hessian;	0.1955095825	1.1e-06	272	0.01
relax(levels=3)	0.1955119036	1.4e-10	256	0.03	yes	yes

	energy	residual	facets	seconds	converged	stable	notes
recipe: g 5; r; g 5; r; g 5; r; g 10; r; g 10; hessian; hessian;	0.3199260248	1.7e-10	256	0.03
relax(levels=4)	0.3199242245	7.4e-12	256	0.06	yes	yes

	energy	residual	facets	seconds	converged	stable	notes
recipe: r; u; g 5; refine edge where valence==1; u; V; g 5; r; g 12; u; V; convert_to_quantities; hessian; hessian; r; hessian; hessian;	0.2889451228	3.0e-05	384	0.01
relax(levels=3)	0.2891108142	2.1e-11	256	0.05	yes	yes

relax() agrees with the recipes of P and D to about 1e-5. On I-WP the recipe refines the edges along the boundary once more and stops unconverged, so the two end on different meshes.

Next

Drops against walls, gravity and solder: Brakke’s samples: drops and capillarity. The tutorial builds a double bubble from arrays instead of a datafile.