Nodal Approximations for Triply Periodic Minimal Surfaces

Ulrich Schwarz, 29 October 2001

This data and the way it was generated were published in U. S. Schwarz and G. Gompper, Systematic approach to bicontinuous cubic phases in ternary amphiphilic systems, Phys. Rev. E 59: 5528 - 5541(1999).
Each structure corresponds to a field phi(x,y,z), and surfaces are constructed as isosurfaces to the value phi=0 of this field. Here we give the nodal approximations (which represents only
the space group information and not the minimal surface property) for different triply periodic minimal surfaces: P, G, D, I-WP, F-RD, S, C(P) and O,C-TO. [C(D) is missing since the formula given in PRE seems to be incorrect.]
 For the physically most relevant ones (P, D, G and I-WP) we also give the improved nodal approximations, which are approximations to the real minimal surfaces and for which the maximum in mean curvature over the surface is 
 reduced by one order of magnitude compared with the nodal approximation. We also give formulae for Gaussian curvature K and mean curvature H, which can be evaluated at the positions of the isosurface, and even used 
 to color the surface for its curvature properties (unfortunately, this does not look so nice in Mathematica, and one has to use many triangles anyway, there are better programs for doing this).

Needs["Graphics`ContourPlot3D`"]
Needs["Graphics`Graphics`"]

Ia[pf_, qf_, rf_, h_, k_, l_, x_, y_, z_]  :=  
      pf[h x] qf[k y] rf[l z]
Ea[pf_, qf_, rf_, h_, k_, l_, x_, y_, 
        z_]  :=   pf[h x] qf[k y] rf[l z] +   pf[h y] qf[
            k z] rf[l x] +  
        pf[h z] qf[k x] rf[l y]
Oa[pf_, qf_, rf_, h_, k_, l_, x_, y_, 
        z_]  :=   pf[h x] qf[k z] rf[l y] +   pf[h y] qf[
            k x] rf[l z] +  pf[h z] qf[k y] rf[l x]

Nodal Approximations with Space Group Information only

NodalP[x_,y_,z_] := Cos[x] + Cos[y] + Cos[z];
NodalG[x_,y_,z_] := Sin[x] Cos[y] + Sin[y] Cos[z] + Sin[z] Cos[x];
NodalD[x_,y_,z_] := Cos[x-y] Cos[z] + Sin[x+y] Sin[z];
NodalIWP[x_,y_,z_] := 2 (Cos[x] Cos[y] + Cos[y] Cos[z] + Cos[x] Cos[z]) - (Cos[2 x] + Cos[2 y] + Cos[2 z]);
NodalFRD[x_,y_,z_] := 4 (Cos[x] Cos[y] Cos[z]) - (Cos[2 x] Cos[2 y] + Cos[2 y] Cos[2 z] + 
          Cos[2z] Cos[2 x]);
NodalS[x_,y_,z_] := Cos[2 x] Sin[y] Cos[z] + Cos[x] Cos[2 y] Sin[z] + Sin[x] Cos[y] Cos[2 z];
NodalCP[x_,y_,z_] := Cos[x] + Cos[y] + Cos[z] + 4 Cos[x] Cos[y] Cos[z];
NodalOCTO[x_,y_,z_] := 4 (Cos[x] Cos[y] + Cos[y] Cos[z] + Cos[x] Cos[z]) - (Cos[2 x] + 
          Cos[2 y] + Cos[2 z]) + 8  Cos[x] Cos[y] Cos[z];

Improved Nodal Approximations from Minimization of Ginzburg-Landau Hamiltonian

FullP[x_, y_, z_]  :=  
      0.2260 (Ea[Cos, Cos, Cos, 1, 0, 0, x, y, z] + 
              Oa[Cos, Cos, Cos, 1, 0, 0, x, y, z])
 - 
        0.0516 (Ea[Cos, Cos, Cos, 1, 1, 1, x, y, z] + 
              Oa[Cos, Cos, Cos, 1, 1, 1, x, y, z])
 - 
        0.0196 (Ea[Cos, Cos, Cos, 2, 1, 0, x, y, z] + 
              Oa[Cos, Cos, Cos, 2, 1, 0, x, y, z])
 - 
        0.0027 (Ea[Cos, Cos, Cos, 3, 0, 0, x, y, z] + 
              Oa[Cos, Cos, Cos, 3, 0, 0, x, y, z])
FullIWP[x_, y_, z_]  :=  
      0.0652  +  
        0.5739 Ea[Cos, Cos, Cos, 1, 1, 0, x, y, z]
 - 
        0.1712 Ea[Cos, Cos, Cos, 2, 0, 0, x, y, z]
 - 
        0.1314 Ea[Cos, Cos, Cos, 2, 1, 1, x, y, z]
 + 
        0.0184 Ea[Cos, Cos, Cos, 2, 2, 0, x, y, z]
FullG[x_, y_, z_]  :=  
      0.3435 (Ea[Sin, Cos, Cos, 1, 1, 0, x, y, z] + 
              Oa[Cos, Sin, Cos, 1, 1, 0, x, y, z])
 -  
        0.0202 (Ea[Cos, Sin, Cos, 0, 3, 1, x, y, z] - 
              Oa[Cos, Cos, Sin, 0, 3, 1, x, y, z])
 -  
        0.0106 (Ea[Sin, Sin, Sin, 2, 2, 2, x, y, z] + 
              Oa[Sin, Sin, Sin, 2, 2, 2, x, y, z])
 +  
        0.0298 (Ea[Cos, Sin, Cos, 2, 1, 3, x, y, z] + 
              Oa[Cos, Cos, Sin, 2, 1, 3, x, y, z])
 +  
        0.0016 (Ea[Cos, Cos, Sin, 3, 0, 3, x, y, z] + 
              Oa[Sin, Cos, Cos, 3, 0, 3, x, y, z])
 -  
        0.0011 (Ea[Sin, Cos, Cos, 1, 1, 4, x, y, z] + 
              Oa[Cos, Sin, Cos, 1, 1, 4, x, y, z])
FullD[x_, y_, z_]  :=  
      0.1407 (Ea[Cos, Cos, Cos, 1, 1, 1, x, y, z] + 
              Ea[Cos, Sin, Sin, 1, 1, 1, x, y, z] + 
              Ea[Sin, Cos, Sin, 1, 1, 1, x, y, z] + 
              Ea[Sin, Sin, Cos, 1, 1, 1, x, y, z] + 
              Oa[Cos, Cos, Cos, 1, 1, 1, x, y, z] + 
              Oa[Cos, Sin, Sin, 1, 1, 1, x, y, z] + 
              Oa[Sin, Cos, Sin, 1, 1, 1, x, y, z] + 
              Oa[Sin, Sin, Cos, 1, 1, 1, x, y, z])
 + 
        0.0209 (Ea[Cos, Cos, Cos, 3, 1, 1, x, y, z] + 
              Ea[Cos, Sin, Sin, 3, 1, 1, x, y, z] - 
              Ea[Sin, Cos, Sin, 3, 1, 1, x, y, z] - 
              Ea[Sin, Sin, Cos, 3, 1, 1, x, y, z] + 
              Oa[Cos, Cos, Cos, 3, 1, 1, x, y, z] + 
              Oa[Cos, Sin, Sin, 3, 1, 1, x, y, z] - 
              Oa[Sin, Cos, Sin, 3, 1, 1, x, y, z] - 
              Oa[Sin, Sin, Cos, 3, 1, 1, x, y, z])
 - 
        0.0138 (Ea[Cos, Cos, Cos, 3, 1, 3, x, y, z] - 
              Ea[Cos, Sin, Sin, 3, 1, 3, x, y, z] + 
              Ea[Sin, Cos, Sin, 3, 1, 3, x, y, z] - 
              Ea[Sin, Sin, Cos, 3, 1, 3, x, y, z] + 
              Oa[Cos, Cos, Cos, 3, 1, 3, x, y, z] - 
              Oa[Cos, Sin, Sin, 3, 1, 3, x, y, z] + 
              Oa[Sin, Cos, Sin, 3, 1, 3, x, y, z] - 
              Oa[Sin, Sin, Cos, 3, 1, 3, x, y, z])
 - 
        0.0028 (Ea[Cos, Cos, Cos, 1, 1, 5, x, y, z] + 
              Ea[Cos, Sin, Sin, 1, 1, 5, x, y, z] + 
              Ea[Sin, Cos, Sin, 1, 1, 5, x, y, z] + 
              Ea[Sin, Sin, Cos, 1, 1, 5, x, y, z] + 
              Oa[Cos, Cos, Cos, 1, 1, 5, x, y, z] + 
              Oa[Cos, Sin, Sin, 1, 1, 5, x, y, z] + 
              Oa[Sin, Cos, Sin, 1, 1, 5, x, y, z] + 
              Oa[Sin, Sin, Cos, 1, 1, 5, x, y, z])
 - 
        0.0021 (Ea[Cos, Cos, Cos, 3, 3, 3, x, y, z] + 
              Ea[Cos, Sin, Sin, 3, 3, 3, x, y, z] + 
              Ea[Sin, Cos, Sin, 3, 3, 3, x, y, z] + 
              Ea[Sin, Sin, Cos, 3, 3, 3, x, y, z] + 
              Oa[Cos, Cos, Cos, 3, 3, 3, x, y, z] + 
              Oa[Cos, Sin, Sin, 3, 3, 3, x, y, z] + 
              Oa[Sin, Cos, Sin, 3, 3, 3, x, y, z] + 
              Oa[Sin, Sin, Cos, 3, 3, 3, x, y, z])
 + 
        0.0011 (Ea[Cos, Cos, Cos, 3, 1, 5, x, y, z] + 
              Ea[Cos, Sin, Sin, 3, 1, 5, x, y, z] - 
              Ea[Sin, Cos, Sin, 3, 1, 5, x, y, z] - 
              Ea[Sin, Sin, Cos, 3, 1, 5, x, y, z] + 
              Oa[Cos, Cos, Cos, 3, 1, 5, x, y, z] + 
              Oa[Cos, Sin, Sin, 3, 1, 5, x, y, z] - 
              Oa[Sin, Cos, Sin, 3, 1, 5, x, y, z] - 
              Oa[Sin, Sin, Cos, 3, 1, 5, x, y, z])

Curvature Properties

Mean curvature

H[phi_, x_, y_, z_] := -(Derivative[0, 0, 2][phi][x, y, z] ((Derivative[0, 1, 0][phi][x, y, z])^2)-2 Derivative[0, 0, 1][phi][x, y, z] Derivative[0, 1, 0][phi][x, y, z] Derivative[0, 1, 1][phi][x, y, z]+((Derivative[0, 0, 1][phi][x, y, z])^2) Derivative[0, 2, 0][phi][x, y, z]+Derivative[0, 0, 2][phi][x, y, z] ((Derivative[1, 0, 0][phi][x, y, z])^2)+Derivative[0, 2, 0][phi][x, y, z] ((Derivative[1, 0, 0][phi][x, y, z])^2)-2 Derivative[0, 0, 1][phi][x, y, z] Derivative[1, 0, 0][phi][x, y, z] Derivative[1, 0, 1][phi][x, y, z]-2 Derivative[0, 1, 0][phi][x, y, z] Derivative[1, 0, 0][phi][x, y, z] Derivative[1, 1, 0][phi][x, y, z]+((Derivative[0, 0, 1][phi][x, y, z])^2) Derivative[2, 0, 0][phi][x, y, z]+((Derivative[0, 1, 0][phi][x, y, z])^2) Derivative[2, 0, 0][phi][x, y, z])/(((Derivative[0, 0, 1][phi][x, y, z])^2)+((Derivative[0, 1, 0][phi][x, y, z])^2)+((Derivative[1, 0, 0][phi][x, y, z])^2))^(3/2)/2

Gaussian curvature

K[phi_, x_, y_, z_] := -(((Derivative[0, 1, 1][phi][x, y, z])^2) ((Derivative[1, 0, 0][phi][x, y, z])^2)-Derivative[0, 0, 2][phi][x, y, z] Derivative[0, 2, 0][phi][x, y, z] ((Derivative[1, 0, 0][phi][x, y, z])^2)-2 Derivative[0, 1, 0][phi][x, y, z] Derivative[0, 1, 1][phi][x, y, z] Derivative[1, 0, 0][phi][x, y, z] Derivative[1, 0, 1][phi][x, y, z]+2 Derivative[0, 0, 1][phi][x, y, z] Derivative[0, 2, 0][phi][x, y, z] Derivative[1, 0, 0][phi][x, y, z] Derivative[1, 0, 1][phi][x, y, z]+((Derivative[0, 1, 0][phi][x, y, z])^2) ((Derivative[1, 0, 1][phi][x, y, z])^2)+2 Derivative[0, 0, 2][phi][x, y, z] Derivative[0, 1, 0][phi][x, y, z] Derivative[1, 0, 0][phi][x, y, z] Derivative[1, 1, 0][phi][x, y, z]-2 Derivative[0, 0, 1][phi][x, y, z] Derivative[0, 1, 1][phi][x, y, z] Derivative[1, 0, 0][phi][x, y, z] Derivative[1, 1, 0][phi][x, y, z]-2 Derivative[0, 0, 1][phi][x, y, z] Derivative[0, 1, 0][phi][x, y, z] Derivative[1, 0, 1][phi][x, y, z] Derivative[1, 1, 0][phi][x, y, z]+((Derivative[0, 0, 1][phi][x, y, z])^2) ((Derivative[1, 1, 0][phi][x, y, z])^2)-Derivative[0, 0, 2][phi][x, y, z] ((Derivative[0, 1, 0][phi][x, y, z])^2) Derivative[2, 0, 0][phi][x, y, z]+2 Derivative[0, 0, 1][phi][x, y, z] Derivative[0, 1, 0][phi][x, y, z] Derivative[0, 1, 1][phi][x, y, z] Derivative[2, 0, 0][phi][x, y, z]-((Derivative[0, 0, 1][phi][x, y, z])^2) Derivative[0, 2, 0][phi][x, y, z] Derivative[2, 0, 0][phi][x, y, z])/(((Derivative[0, 0, 1][phi][x, y, z])^2)+((Derivative[0, 1, 0][phi][x, y, z])^2)+((Derivative[1, 0, 0][phi][x, y, z])^2))^2

Contour Plot

f  :=  NodalG

np  := 5;

pp  =  
      ContourPlot3D[
        f[x, y, z], {x, -Pi, Pi}, {y, -Pi, Pi}, {z, -Pi, Pi}, 
        ViewPoint -> {-1.160,  -2.875,  1.356}, 
        PlotPoints -> {np, np, np}]

Contour Plot with Curvature Properties in Color Code

mycolor[{x_, y_, z_}] := SurfaceColor[Hue[K[f, x, y, z]]]

bild  =  
      Show[Graphics3D[{EdgeForm[{}], 
            Cases[pp, _Polygon, Infinity] /. 
              Polygon[pnts_] :> {mycolor[(Plus @@ pnts)/
                      Length[pnts]], Polygon[pnts]}}], 
        ViewPoint -> {2.183,  0.938,  2.409}]

Area and Euler characteristic

Extract midpoints of the polygons

MidList  =  
        Cases[pp, _Polygon, Infinity] /. 
          Polygon[pnts_] :> (Plus @@ pnts)/
              Length[pnts];

List of Gaussian curvatures at the midpoints

KK[{x_, y_, z_}]  :=  K[f, x, y, z]

KList  =  Map[KK, MidList];

Kmin  =  Min[KList]
Kmax  =  Max[KList]

List of corresponding surface areas

polyarea[p_]  :=  
      Module[{n, mid, c}, n = Length[p]; mid  = (Plus  @@  p)/n; 
        Sum[
c  =  
            Cross[Part[p, Mod[i + 1, n, 1]] - mid, 
              Part[p, Mod[i, n, 1]] - mid]; 
          Sqrt[c . c]/2, {i, 1, n}]]

AreaList  =  
        Cases[pp, _Polygon, Infinity] /. 
          Polygon[pnts_] :>  polyarea[pnts];

Surface area

TotalArea  =  (Plus  @@  AreaList)/(2 Pi)^2

Euler characteristic

chi  =  (Sum[
            Part[KList, i] Part[AreaList, i], {i, 1, 
              Length[KList]}])/(2 Pi)