# Approximation of triply periodic minimal surfaces with trigonometric level sets

Some triply periodic minimal surfaces are known to be approximated by trigonometric level sets very accurately.

To see this, let's sample a gyroid scaled to the bounding box $$[0, 1]^3$$ exactly through its Weierstrass-Enneper parametrization. We can examine the distance to the nearest point on the level set

$$\cos(2\pi x)\sin(2\pi y) - \cos(2\pi y)\sin(2\pi z) + \cos(2\pi z)\sin(2\pi x) = 0:$$

The median relative error is ~$$0.0004$$ and max relative error is ~$$0.0013$$, and so this is a very good approximation for most applications.

A list of level sets that approximate other TPMS can be found here.

## Question

How were these approximations derived or discovered? Perhaps through a series approximation to the Weierstrass-Enneper parametrization?

Are there ways to find more accurate approximations?

## Edit

I found the following excerpt here.

These functions are not found by chance; they come from choosing an appropriate lowest-order term from the Fourier series of an electrostatic potential function derived from charges whose distribution has the desired space-group symmetry, or from a symmetrization procedure using generators of the space group. But there is as yet no real explanation as to why the match is so good in some cases and not at all accurate in others.

It’s not clear to me what this means exactly, and unfortunately there is no source listed.

I believe these equations first appear in Ian Barnes' 1990 PhD thesis Microstructure of bicontinuous phases in surfactant systems.

In section 6.2.1 he derives the Fourier representation of periodic equipotential surfaces:

$$\varphi(x, y, z) = \frac{8}{\pi}\sum_{j=1}^N q_j \sum_{h,k,l=0}^\infty{\vphantom{\sum}}{\!\!\!'\,\,\,} \frac{\cos(2\pi h(x-x_j))\cos(2\pi k(y-y_j))\cos(2\pi l(z-z_j))}{h^2+k^2+l^2},$$

where the prime on the last sum indicates that for each index, terms are to be multiplied by $$1/2$$ if it is zero, and the $$h = k = l = 0$$ term is taken to be $$0$$.

Here there are $$N$$ charges in $$[0, 1]^3$$ with values $$q_j$$ at points $$(x_j, y_j, z_j)$$, and these charges are then periodically tiled.

Explicit examples are covered in section 6.3 by truncating the above infinite sum, as first done in a paper by Mackay. For the gyroid, in section 6.3.5, there are 16 charges:

• 8 points with charge $$+1$$ at $$\small \left(\frac{1}{8},\frac{1}{8},\frac{1}{8}\right), \left(\frac{5}{8},\frac{3}{8},\frac{7}{8}\right), \left(\frac{7}{8},\frac{5}{8},\frac{3}{8}\right), \left(\frac{3}{8},\frac{7}{8},\frac{5}{8}\right), \left(\frac{5}{8},\frac{5}{8},\frac{5}{8}\right), \left(\frac{1}{8},\frac{7}{8},\frac{3}{8}\right), \left(\frac{3}{8},\frac{1}{8},\frac{7}{8}\right), \left(\frac{7}{8},\frac{3}{8},\frac{1}{8}\right)$$

• 8 points with charge $$-1$$ at $$\small \left(\frac{7}{8},\frac{7}{8},\frac{7}{8}\right), \left(\frac{3}{8},\frac{5}{8},\frac{1}{8}\right), \left(\frac{1}{8},\frac{3}{8},\frac{5}{8}\right), \left(\frac{5}{8},\frac{1}{8},\frac{3}{8}\right), \left(\frac{3}{8},\frac{3}{8},\frac{3}{8}\right), \left(\frac{7}{8},\frac{1}{8},\frac{5}{8}\right), \left(\frac{5}{8},\frac{7}{8},\frac{1}{8}\right), \left(\frac{1}{8},\frac{5}{8},\frac{7}{8}\right)$$

We can verify the well known equation that approximates the gyroid by truncating at $$h,k,l\leq1$$ in Mathematica:

Do[q[j] = 1, {j, 8}]

Do[q[j + 8] = -1, {j, 8}]

pts = {
{1, 1, 1}, {5, 3, 7}, {7, 5, 3}, {3, 7, 5},
{5, 5, 5}, {1, 7, 3}, {3, 1, 7}, {7, 3, 1},
{7, 7, 7}, {3, 5, 1}, {1, 3, 5}, {5, 1, 3},
{3, 3, 3}, {7, 1, 5}, {5, 7, 1}, {1, 5, 7}
}/8;

Do[{xx[j], yy[j], zz[j]} = pts[[j]], {j, 16}]

gyrapprox[n_] := FullSimplify[
8/π Sum[q[j] (Sum[Piecewise[{
{((1/2)^(Boole[h == 0] + Boole[k == 0] + Boole[l == 0])*
Cos[2π h (x - xx[j])] Cos[2π k (y - yy[j])] Cos[2π l (z - zz[j])])/(h^2 + k^2 + l^2),
h != 0 || k != 0 || l != 0}}], {h, 0, n}, {k, 0, n}, {l, 0, n}]), {j, 16}],
ComplexityFunction -> (LeafCount[#] + 1000 Total[Cases[#, p_Plus :> Length[p], ∞]] &)]


The approximate formula (after dividing both sides by a constant):

g1 = Simplify[π/16 gyrapprox[1] /. {v : x | y | z -> v/(2π)}] == 0

Cos[y]Sin[x] + Cos[z]Sin[y] + Cos[x]Sin[z] == 0


We can take more terms now:

g2 = Simplify[π/16 gyrapprox[2] /. {v : x | y | z -> v/(2π)}] == 0

Cos[y]Sin[x] + Cos[z]Sin[y] + Cos[x]Sin[z] + 2Sin[2x]Sin[2y]Sin[2z]/3 == 0


Omitting the output:

g3 = Simplify[π/16 gyrapprox[3] /. {v : x | y | z -> v/(2π)}] == 0;


Unfortunately only g1, the well known formula, gives an accurate approximation to the gyroid:

opts = {PlotTheme -> "Minimal", Axes -> False, Boxed -> False};

ContourPlot3D[#, {x, 0, 2π}, {y, 0, 2π}, {z, 0, 2π}, opts] & /@ ({g1, g2, g3} /. y -> -y)