{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "brakke_lagrange-0",
   "metadata": {},
   "source": [
    "# Brakke's samples: quadratic and Lagrange elements\n",
    "\n",
    "\n",
    "Curved elements: the quadratic model (edges with midpoints) and Lagrange elements of\n",
    "higher order, a short final stage on a mesh settled in the linear model. `relax()` in\n",
    "those models first does what the manual recommends, five gradient steps and then\n",
    "Newton steps.\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_lagrange-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_lagrange-2",
   "metadata": {},
   "source": [
    "## `quad.fe`: a skew quadrilateral\n",
    "\n",
    "A film on a skew quadrilateral wire in the quadratic model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_lagrange-3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bh.recipe(\"quad.fe\"))\n",
    "runs = bh.compare(\"quad.fe\")\n",
    "bh.show(runs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_lagrange-4",
   "metadata": {},
   "source": [
    "## `qmound.fe`: the drop on a table, quadratic\n",
    "\n",
    "`mound.fe` in the quadratic model, its volume a named quantity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_lagrange-5",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bh.recipe(\"qmound.fe\"))\n",
    "runs = bh.compare(\"qmound.fe\")\n",
    "bh.show(runs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_lagrange-6",
   "metadata": {},
   "source": [
    "## The Lagrange ladders of `cube.fe` and `catman.fe`\n",
    "\n",
    "`cube.fe`'s `gogo2` refines twice in the linear model, then switches to Lagrange\n",
    "elements of order 2, 4 and 6 with `g 5; hessian` at each. The same with `relax()`:\n",
    "relax with two refinements, then for each order `ev.set_model(\"lagrange\", n)` and\n",
    "relax."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_lagrange-7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "def ladder(name, orders, levels=2, setup=None):\n",
    "    ev = pyse.Evolver(name)\n",
    "    if setup:\n",
    "        setup(ev)\n",
    "    t = time.perf_counter()\n",
    "    ev.relax(levels=levels)\n",
    "    for n in orders:\n",
    "        ev.set_model(\"lagrange\", n)\n",
    "        r = ev.relax()\n",
    "    return bh.Run(f\"relax ladder, orders {orders}\", ev.total_energy, r.residual,\n",
    "                  ev.counts[\"facets\"], time.perf_counter() - t, r.converged, r.stable, [],\n",
    "                  ev.mesh())\n",
    "\n",
    "print(bh.recipe(\"cube.fe\", \"gogo2\"))\n",
    "cube = [bh.run_recipe(\"cube.fe\", \"gogo2\", label=\"recipe: gogo2\"), ladder(\"cube.fe\", (2, 4, 6))]\n",
    "bh.table(cube, reference=(36*np.pi)**(1/3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_lagrange-8",
   "metadata": {},
   "source": [
    "At order 6 both match the sphere's area to rounding error (the recipe's own check\n",
    "prints an area error of about 1e-15): sixth-order elements on this mesh represent the\n",
    "sphere essentially exactly.\n",
    "\n",
    "`catman.fe`'s `gogo3` goes up to order 8. It sets `rmax := cosh(height)` and compares\n",
    "with the area of a catenoid through rings of radius $\\cosh(0.55)$, but `rmax` is a\n",
    "parameter of `cat.fe`: in `catman.fe` the rings keep radius 1, so its printed\n",
    "\"Difference from true area\" compares with the wrong catenoid. The exact catenoid for\n",
    "rings of radius 1 at $\\pm 0.55$ has area 6.5116735 (from $a \\cosh(0.55/a) = 1$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brakke_lagrange-9",
   "metadata": {},
   "outputs": [],
   "source": [
    "lo, hi = 0.6, 1.0                   # the larger root of a cosh(0.55/a) = 1: the stable catenoid\n",
    "for _ in range(100):\n",
    "    a = (lo + hi)/2\n",
    "    lo, hi = (lo, a) if a*np.cosh(0.55/a) > 1 else (a, hi)\n",
    "exact = np.pi*a*(2*0.55 + a*np.sinh(2*0.55/a))\n",
    "print(f\"exact catenoid area {exact:.10f}\")\n",
    "catman = [bh.run_recipe(\"catman.fe\", \"gogo3\", label=\"recipe: gogo3\"),\n",
    "          ladder(\"catman.fe\", (2, 4, 6, 8))]\n",
    "bh.table(catman, reference=exact)\n",
    "bh.show(catman)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brakke_lagrange-10",
   "metadata": {},
   "source": [
    "## Next\n",
    "\n",
    "The [tutorial](tutorial.ipynb) uses the same Lagrange ladder as a final stage for a\n",
    "double bubble built from arrays, and [Performance](performance.md) times it on large\n",
    "surfaces. The [API reference](api.md) documents `set_model` and `relax`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}