{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "brakke_films-0",
   "metadata": {},
   "source": [
    "# Brakke's samples: soap films and foams\n",
    "\n",
    "\n",
    "Soap films spanning wires, a bubble, and foams in a periodic box, from the sample\n",
    "datafiles that come with Surface Evolver (`pyse.examples`). The helpers are in\n",
    "[`docs/brakke_helpers.py`](https://github.com/gautelinga/pySurfaceEvolver/blob/main/docs/brakke_helpers.py).\n",
    "\n",
    "Each sample is run twice from its datafile: once with the recipe Brakke wrote into\n",
    "it (its `gogo` command, or the usage its comments describe), once with a default\n",
    "`ev.relax()` given as many refinements as the recipe does (`levels`). The table\n",
    "shows the energy each ends with, the residual (how far from equilibrium, see\n",
    "[Relaxing a surface](relaxing.md)), the facet count, the time, and what `relax()`\n",
    "reports: converged, stable, and any health issues. A lower energy on the same mesh\n",
    "is not always the better answer: a coarse mesh can lower its energy in ways the\n",
    "smooth surface can't, so where an exact value is known, the error against it is\n",
    "the measure. Times are from the documentation build and vary.\n",
    "\n",
    "This gallery assumes [Getting started](getting_started.md) and\n",
    "[Relaxing a surface](relaxing.md); its pages can be read in any order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_films-1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pyvista as pv\n",
    "import pysurfaceevolver as pyse\n",
    "import brakke_helpers as bh\n",
    "\n",
    "pv.set_jupyter_backend(\"static\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-2",
   "metadata": {},
   "source": [
    "## `cube.fe`: a cube becoming a sphere\n",
    "\n",
    "A unit cube of fixed volume 1, its surface area the energy: it relaxes toward a sphere\n",
    "(area $(36\\pi)^{1/3} \\approx 4.836$). The recipe refines twice with gradient and Newton\n",
    "steps in between."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_films-3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bh.recipe(\"cube.fe\"))\n",
    "runs = bh.compare(\"cube.fe\")\n",
    "bh.table(runs, reference=(36*np.pi)**(1/3))\n",
    "bh.show(runs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-4",
   "metadata": {},
   "source": [
    "Both meshes are polyhedra inscribed in the sphere, with a little more area than it;\n",
    "`relax()` takes its residual to about 1e-9 where the recipe stops sooner, at nearly the\n",
    "same energy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-5",
   "metadata": {},
   "source": [
    "## `catman.fe`: the catenoid of the manual's tutorial\n",
    "\n",
    "A film between two rings of radius 1 at height $\\pm 0.55$; the stable catenoid has\n",
    "area 6.511674. The recipe's comments say what happens on this coarse mesh: after 50\n",
    "steps \"at this point have nearly a saddle point\", then \"triangulation twists around to\n",
    "lower energy\". Twisting lowers the area of the coarse mesh, not its error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_films-6",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bh.recipe(\"catman.fe\"))\n",
    "runs = bh.compare(\"catman.fe\")\n",
    "bh.table(runs, reference=6.511674)\n",
    "bh.show(runs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-7",
   "metadata": {},
   "source": [
    "Both are below the smooth catenoid (the rings are coarse polygons); the twisted mesh\n",
    "of the recipe is further below. `relax()` stops at the untwisted equilibrium, and\n",
    "leaves sliding along the wires (where the twist happens) to its gradient rounds."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-8",
   "metadata": {},
   "source": [
    "## `cat.fe`: a catenoid past its limit\n",
    "\n",
    "Here the rings are as far apart as a catenoid can span (`RMAX` set at the fold). The\n",
    "neck collapses, and the recipe deletes the tiny edges there (`t .05`) and pops the\n",
    "pinched vertex (`o`): the film becomes two flat discs. `relax()` can't converge the\n",
    "film, so as a last resort it takes gradient steps without tidying, lets the neck\n",
    "collapse, deletes the collapsed edges and pops (`collapsed`, `popped`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_films-9",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bh.recipe(\"cat.fe\"))\n",
    "runs = bh.compare(\"cat.fe\")\n",
    "R = 1.5088795                       # RMAX: the rings' radius\n",
    "discs = 2*3*R**2                    # two 12-gons inscribed in the rings (after one refinement)\n",
    "bh.table(runs, reference=discs)\n",
    "bh.show(runs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-10",
   "metadata": {},
   "source": [
    "The recipe stops after `g 5` on the discs; `relax()` finishes them: two regular\n",
    "12-gons, the polygons the rings' vertices make after one refinement."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-11",
   "metadata": {},
   "source": [
    "## `catbody.fe`: a catenoid enclosing a fixed volume\n",
    "\n",
    "The same catenoid shape, now closed by plane ends and holding a fixed volume, so it\n",
    "has a pressure. The file adjusts that volume to find zero pressure (the unstable\n",
    "catenoid) in a separate script; here only its `gogo`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_films-12",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bh.recipe(\"catbody.fe\"))\n",
    "runs = bh.compare(\"catbody.fe\")\n",
    "bh.show(runs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-13",
   "metadata": {},
   "source": [
    "## `octa.fe`: a film in an octahedral frame\n",
    "\n",
    "A symmetric film from the edges of an octahedral wire frame to its centre: six edges\n",
    "there carry four films each. Plateau's laws (soap films meet three at a time along an\n",
    "edge, at 120 degrees, and such edges four at a time at a point) say that is unstable:\n",
    "splitting the junction lowers the area. The recipe pops junctions (`o`), deletes small\n",
    "edges (`t .1`) and iterates, several times; `relax()` pops after converging, relaxes\n",
    "again, and deletes the films that shrink away."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_films-14",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bh.recipe(\"octa.fe\"))\n",
    "runs = bh.compare(\"octa.fe\")\n",
    "bh.show(runs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-15",
   "metadata": {},
   "source": [
    "Both reach the same film, with a central tetrahedral point."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-16",
   "metadata": {},
   "source": [
    "## Kelvin against Weaire-Phelan: `twointor.fe` and `phelanc.fe`\n",
    "\n",
    "Lord Kelvin asked in 1887 how to divide space into cells of equal volume with the\n",
    "least interface area, and proposed the cell of `twointor.fe` (two cells in a\n",
    "periodic box). In 1993 Weaire and Phelan found a structure with two kinds of cells\n",
    "(`phelanc.fe`, eight cells) with about 0.3% less area. Per cell, the area over the\n",
    "cell volume to the power 2/3 compares them independently of the box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_films-17",
   "metadata": {},
   "outputs": [],
   "source": [
    "kelvin = bh.compare(\"twointor.fe\")\n",
    "wp = bh.compare(\"phelanc.fe\")\n",
    "\n",
    "def per_cell(name, levels):\n",
    "    ev = pyse.Evolver(name)\n",
    "    ev.relax(levels=levels)\n",
    "    n = len(ev.bodies().volume)\n",
    "    V = float(np.sum(ev.bodies().volume))\n",
    "    return ev.total_energy/n/(V/n)**(2/3)\n",
    "\n",
    "k, w = per_cell(\"twointor.fe\", 2), per_cell(\"phelanc.fe\", 2)\n",
    "print(f\"area per cell / volume^(2/3): Kelvin {k:.5f}, Weaire-Phelan {w:.5f}, \"\n",
    "      f\"{1 - w/k:.2%} less\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_films-18",
   "metadata": {},
   "outputs": [],
   "source": [
    "bh.show(wp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-19",
   "metadata": {},
   "source": [
    "## `symtest.fe`: the Kelvin foam in a quotient space\n",
    "\n",
    "The same two Kelvin cells, but in Evolver's quotient-space representation (a symmetry\n",
    "group instead of the periodic box). The recipe is `twointor.fe`'s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_films-20",
   "metadata": {},
   "outputs": [],
   "source": [
    "runs = bh.compare(\"symtest.fe\", bh.recipe(\"twointor.fe\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-21",
   "metadata": {},
   "source": [
    "The energy is that of `twointor.fe`, but `relax()` reports the surface unstable, with\n",
    "one negative Hessian eigenvalue (as Brakke's Evolver 2.70a does), and `twointor.fe`\n",
    "stable. The instability is real: the volumes in this file's `bodies` section are\n",
    "comments (`// volume 0.500`), so its two cells have no volume constraint, and one can\n",
    "grow at the expense of the other. Moving along that mode (Evolver's `saddle`) lowers\n",
    "the energy, and relaxing from there shrinks one cell away."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_films-22",
   "metadata": {},
   "outputs": [],
   "source": [
    "ev = pyse.Evolver(\"symtest.fe\")\n",
    "print(ev.bodies().fixed)                    # no volume constraints\n",
    "ev.relax(levels=1)\n",
    "e0 = ev.total_energy\n",
    "ev.saddle()                                 # a step along the negative mode\n",
    "print(f\"energy {e0:.6f} -> {ev.total_energy:.6f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-23",
   "metadata": {},
   "source": [
    "## Triply periodic minimal surfaces: `pcell.fe`, `dcell.fe`, `IWP.fe`\n",
    "\n",
    "From Brakke's collection of periodic minimal surfaces online (included with pySE):\n",
    "Schwarz's P and D surfaces and Schoen's I-WP, each a small piece of the surface in a\n",
    "quotient space (a symmetry group generates the rest)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_films-24",
   "metadata": {},
   "outputs": [],
   "source": [
    "for name in (\"pcell.fe\", \"dcell.fe\", \"IWP.fe\"):\n",
    "    runs = bh.compare(name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-25",
   "metadata": {},
   "source": [
    "`relax()` agrees with the recipes of P and D to about 1e-5. On I-WP the recipe\n",
    "refines the edges along the boundary once more and stops unconverged, so the two\n",
    "end on different meshes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_films-26",
   "metadata": {},
   "source": [
    "## Next\n",
    "\n",
    "Drops against walls, gravity and solder: [Brakke's samples: drops and\n",
    "capillarity](brakke_capillarity.ipynb). The [tutorial](tutorial.ipynb) builds a double\n",
    "bubble from arrays instead of a datafile."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}