Brakke’s samples: quadratic and Lagrange elements

Curved elements: the quadratic model (edges with midpoints) and Lagrange elements of higher order, a short final stage on a mesh settled in the linear model. relax() in those models first does what the manual recommends, five gradient steps and then Newton steps.

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")

quad.fe: a skew quadrilateral

A film on a skew quadrilateral wire in the quadratic model.

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

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

	energy	residual	facets	seconds	converged	stable	notes
recipe: r; g 5; r; g 5; r; g 5; hessian; hessian;	4.318202947	2.6e-11	256	0.00
relax(levels=3)	4.318322217	4.6e-11	256	0.01	yes	yes

qmound.fe: the drop on a table, quadratic

mound.fe in the quadratic model, its volume a named quantity.

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

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

	energy	residual	facets	seconds	converged	stable	notes
recipe: re; g 5; r; g 5; r; g 5; hessian; hessian;	4.12755008	2.8e-05	448	0.09
relax(levels=2)	4.127333992	8.9e-09	384	0.41	yes	yes

The Lagrange ladders of cube.fe and catman.fe

cube.fe’s gogo2 refines twice in the linear model, then switches to Lagrange elements of order 2, 4 and 6 with g 5; hessian at each. The same with relax(): relax with two refinements, then for each order ev.set_model("lagrange", n) and relax.

import time

def ladder(name, orders, levels=2, setup=None):
    ev = pyse.Evolver(name)
    if setup:
        setup(ev)
    t = time.perf_counter()
    ev.relax(levels=levels)
    for n in orders:
        ev.set_model("lagrange", n)
        r = ev.relax()
    return bh.Run(f"relax ladder, orders {orders}", ev.total_energy, r.residual,
                  ev.counts["facets"], time.perf_counter() - t, r.converged, r.stable, [],
                  ev.mesh())

print(bh.recipe("cube.fe", "gogo2"))
cube = [bh.run_recipe("cube.fe", "gogo2", label="recipe: gogo2"), ladder("cube.fe", (2, 4, 6))]
bh.table(cube, reference=(36*np.pi)**(1/3))

g 5; r; g 5; hessian; r; g 5; hessian; lagrange 2; g 5; hessian; lagrange 4; g 5; hessian; lagrange 6; g 5; hessian; ideal_rad := (3*body[1].volume/4/pi)^(1/3); printf "Area error: %g\n",total_area - 4*pi*ideal_rad^2; printf "Vertex radius spread: %g\n", max(vertex,sqrt((x-.5)^2+(y-.5)^2+(z-.5)^2)) - min(vertex,sqrt((x-.5)^2+(y-.5)^2+(z-.5)^2));

	energy	error	residual	facets	seconds	converged	stable	notes
recipe: gogo2	4.835975862	+4.44e-16	5.3e-11	384	0.50
relax ladder, orders (2, 4, 6)	4.835975862	+2.22e-16	1.9e-11	384	0.67	yes	yes

At order 6 both match the sphere’s area to rounding error (the recipe’s own check prints an area error of about 1e-15): sixth-order elements on this mesh represent the sphere essentially exactly.

catman.fe’s gogo3 goes up to order 8. It sets rmax := cosh(height) and compares with the area of a catenoid through rings of radius \(\cosh(0.55)\), but rmax is a parameter of cat.fe: in catman.fe the rings keep radius 1, so its printed “Difference from true area” compares with the wrong catenoid. The exact catenoid for rings of radius 1 at \(\pm 0.55\) has area 6.5116735 (from \(a \cosh(0.55/a) = 1\)).

lo, hi = 0.6, 1.0                   # the larger root of a cosh(0.55/a) = 1: the stable catenoid
for _ in range(100):
    a = (lo + hi)/2
    lo, hi = (lo, a) if a*np.cosh(0.55/a) > 1 else (a, hi)
exact = np.pi*a*(2*0.55 + a*np.sinh(2*0.55/a))
print(f"exact catenoid area {exact:.10f}")
catman = [bh.run_recipe("catman.fe", "gogo3", label="recipe: gogo3"),
          ladder("catman.fe", (2, 4, 6, 8))]
bh.table(catman, reference=exact)
bh.show(catman)

exact catenoid area 6.5116735251

	energy	error	residual	facets	seconds	converged	stable	notes
recipe: gogo3	6.511673526	+1.58e-10	3.7e-11	384	1.60
relax ladder, orders (2, 4, 6, 8)	6.511673526	+1.58e-10	3.6e-11	384	1.73	yes	yes

Next

The tutorial uses the same Lagrange ladder as a final stage for a double bubble built from arrays, and Performance times it on large surfaces. The API reference documents set_model and relax.