{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c00",
   "metadata": {},
   "source": "# Tutorial: a double bubble, from arrays to FEM mesh\n\nTwo bubbles of volume 1 and 2 share a film. Surface tension pulls the films\ntoward minimal total area while the volumes stay fixed: the result is the\n*standard double bubble*, three spherical caps meeting at 120° along a circle.\nThis notebook builds that surface from arrays, relaxes it, refines the answer\nwith high-order (Lagrange) elements, plots it, and exports watertight surfaces\nper bubble for a finite-element mesher. It assumes the drop on a table and the\nconventions of [Getting started](getting_started.md).\n\nIn Jupyter, use `pv.set_jupyter_backend(\"trame\")` (installed with\n`pip install \"pysurfaceevolver[jupyter]\"`) for interactive 3D views; these\ndocs draw static images."
  },
  {
   "cell_type": "code",
   "id": "c01",
   "metadata": {},
   "source": "import numpy as np\nimport pyvista as pv\nimport pysurfaceevolver as pyse\nfrom pysurfaceevolver import Body, Evolver\n\npv.set_jupyter_backend(\"static\")\npv.global_theme.window_size = [600, 400]",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "c02",
   "metadata": {},
   "source": "## 1. Build the starting surface\n\nStart from two unit cubes side by side: body 1 is $x \\in [0, 1]$, body 2 is\n$x \\in [1, 2]$, and the square at $x = 1$ is their shared film. Faces are\nvertex loops, in any order around the square; each body lists its faces. When a\nbody's faces close a surface, as here, `load_arrays` orients them outward\nitself (otherwise give `orientation`: +1 or -1 per face, whether the face's\nnormal by its vertex order points out of the body)."
  },
  {
   "cell_type": "code",
   "id": "c03",
   "metadata": {},
   "source": "vertices = np.array([(x, y, z) for x in (0, 1, 2) for y in (0, 1) for z in (0, 1)], float)\n\ndef square(axis, value, xmin=0, xmax=2):\n    # the face {coordinate `axis` = value} within xmin <= x <= xmax, as a loop\n    rows = [i for i, v in enumerate(vertices) if v[axis] == value and xmin <= v[0] <= xmax]\n    a, b = [k for k in range(3) if k != axis]\n    center = vertices[rows].mean(axis=0)\n    angle = np.arctan2(vertices[rows, b] - center[b], vertices[rows, a] - center[a])\n    return [rows[i] for i in np.argsort(angle)]\n\nfaces = [square(0, 0), square(0, 1), square(0, 2)]          # the film between: faces[1]\nfor lo, hi in ((0, 1), (1, 2)):\n    faces += [square(axis, value, lo, hi) for axis in (1, 2) for value in (0, 1)]\n\nbodies = [Body(faces=[0, 1, 3, 4, 5, 6], volume=1.0),\n          Body(faces=[1, 2, 7, 8, 9, 10], volume=2.0)]\nev = Evolver()\nev.load_arrays(vertices, faces, bodies=bodies)\nev",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "c04",
   "metadata": {},
   "source": "Evolver triangulated the squares (each got a center vertex). Body 2 starts with\nvolume 1 but its target is 2, so the first iterations inflate it."
  },
  {
   "cell_type": "code",
   "id": "c05",
   "metadata": {},
   "source": "ev.plot(color=\"lightblue\")",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "c06",
   "metadata": {},
   "source": "## 2. Relax\n\n`iterate` takes gradient-descent steps. `relax` brings the surface to an\nequilibrium: gradient steps that keep the mesh in shape, then Newton steps, until\nthe remaining force imbalance (`ev.residual()`) is small. It also checks that the\nequilibrium is stable, and warns if it isn't. `levels=3` refines and relaxes again\nthree times: refining in steps, relaxing after each, keeps the surface close to\nequilibrium. [Relaxing a surface](relaxing.md) says what `relax()` does on the way.\nWithout walls, `ev.energy` is simply the total area."
  },
  {
   "cell_type": "code",
   "id": "c07",
   "metadata": {},
   "source": "history = [ev.iterate(30)]\nresult = ev.relax(levels=3)\nhistory.append(result)\nprint(f\"{ev.counts['facets']} facets, energy {ev.energy:.10f}, converged {result.converged} \"\n      f\"(residual {result.residual:.1e}), {result.newton_steps} Newton steps, stable {result.stable}\")",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "c08",
   "metadata": {},
   "source": "import matplotlib.pyplot as plt\n\nenergy = np.concatenate([h.energy for h in history])\nfig, ax = plt.subplots(figsize=(6, 3))\nax.plot(energy)\nax.set_xlabel(\"iteration\")\nax.set_ylabel(\"energy (total area)\")\noffset = len(history[0].energy)\nfor s in offset + np.flatnonzero(np.diff(result.level)) + 1:\n    ax.axvline(s, color=\"gray\", lw=0.5)    # refinements\nax.set_ylim(energy[-1] - 0.05, energy[-1] + 0.6);",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "c09",
   "metadata": {},
   "source": "The volume constraints have Lagrange multipliers: the bubbles' pressures. By\nthe Young–Laplace law the pressure jump across a film is twice its mean\ncurvature (the energy is area, so the tension is 1): the outer films have mean\ncurvature $p_i/2$, and the shared film $(p_1 - p_2)/2$. The smaller bubble has the\nhigher pressure, so the shared film bulges into the larger one."
  },
  {
   "cell_type": "code",
   "id": "c10",
   "metadata": {},
   "source": "ev.bodies()",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "c11",
   "metadata": {},
   "source": "mesh = ev.mesh()\nshared = (mesh.facet_bodies > 0).all(axis=1)          # films with a body on both sides\nH = ev.values(\"vertex\", \"mean_curvature\")\np1, p2 = ev.bodies().pressure\nprint(f\"pressures {p1:.4f}, {p2:.4f}: mean curvatures p1/2 = {p1/2:.4f}, \"\n      f\"p2/2 = {p2/2:.4f}, (p1 - p2)/2 = {(p1 - p2)/2:.4f}\")\n# the shared film sits inside: draw the outer films translucent\npoly = mesh.to_pyvista(point_values={\"mean curvature\": H})\nplotter = pv.Plotter()\nstyle = dict(scalars=\"mean curvature\", clim=[0, 1.6], cmap=\"viridis\")\nplotter.add_mesh(poly.extract_cells(np.flatnonzero(~shared)), opacity=0.5, **style)\nplotter.add_mesh(poly.extract_cells(np.flatnonzero(shared)), show_edges=True, **style)\nplotter.show()",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "c12",
   "metadata": {},
   "source": "The colors show the three caps of constant mean curvature; the shared film,\ndrawn opaque, is the flattest. (Vertices on the triple junction, where three films\nmeet, have no meaningful mean curvature.) A live view redraws while a run\ngoes on; in Jupyter with the trame backend it updates in place:\n\n```python\nview = ev.live_view(\"area\")\nev.iterate(200, callback=view.update, every=10)\n```\n\n## 3. High-order final stage\n\nThe linear model approximates each cap by flat triangles, so the area\nconverges only as the mesh is refined. Once the triangulation is settled,\nLagrange elements of order $n$ give much more accurate energies on the same\nmesh. In those models `relax()` does what Brakke's manual recommends, five gradient\nsteps and then Newton steps. Compare with one more linear refinement:"
  },
  {
   "cell_type": "code",
   "id": "c13",
   "metadata": {},
   "source": "settled = ev.save()                       # the linear equilibrium\nlinear = {ev.counts[\"facets\"]: ev.energy}\nev.refine(); ev.relax()\nlinear[ev.counts[\"facets\"]] = ev.energy\n\nev.restore(settled)\nlagrange = {}\nfor order in (2, 3, 4):\n    ev.set_model(\"lagrange\", order)\n    ev.relax()\n    lagrange[order] = ev.energy\n\ncoarse = min(linear)\nfor facets, e in linear.items():\n    print(f\"linear      {facets:6d} facets  energy {e:.12f}\")\nfor order, e in lagrange.items():\n    print(f\"Lagrange {order}  {coarse:6d} facets  energy {e:.12f}\")",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "c14",
   "metadata": {},
   "source": "Lagrange 3 on the coarser mesh already agrees with Lagrange 4 to about\na few $10^{-9}$, far better than the linear model on four times as many facets.\nPlots and exports sample the curved facets; the facet sizes still show where\nthe cube's corners were:"
  },
  {
   "cell_type": "code",
   "id": "c15",
   "metadata": {},
   "source": "ev.plot(\"area\", n=4, show_edges=False, cmap=\"plasma\",\n        scalar_bar_args={\"title\": \"facet area\"})",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "c16",
   "metadata": {},
   "source": "## 4. Export for finite elements\n\n`body_surfaces` gives each bubble's closed surface with outward normals (the\nshared film belongs to both). `write_bodies` writes one file per body, either\nas flat triangles sampled on the curved facets (`curved=\"tessellate\"`) or as\nnative high-order triangles (`curved=\"native\"`; Gmsh format, here quadratic\nsince Gmsh meshers take order-2 surfaces)."
  },
  {
   "cell_type": "code",
   "id": "c17",
   "metadata": {},
   "source": "for s in ev.body_surfaces(n=4).values():\n    print(f\"body {s.body}: {len(s.cells)} triangles, watertight {s.watertight}\")",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "c18",
   "metadata": {},
   "source": "import os, tempfile\nimport meshio\n\nout = tempfile.mkdtemp()\nev.write_bodies(os.path.join(out, \"bubble_{id}.stl\"), n=4)\n\nev.set_model(\"lagrange\", 2)\nev.relax()\nev.write_bodies(os.path.join(out, \"bubble_{id}.msh\"), curved=\"native\")\nprint(sorted(os.listdir(out)))\n\ndef enclosed_volume(points, triangles):\n    a, b, c = (points[triangles[:, i]] for i in range(3))\n    return np.einsum(\"ij,ij->i\", a, np.cross(b, c)).sum() / 6\n\nimport warnings\nwarnings.filterwarnings(\"ignore\", \"overflow\", RuntimeWarning)   # meshio's STL reader\n\nfor body in (1, 2):\n    stl = meshio.read(os.path.join(out, f\"bubble_{body}.stl\"))\n    msh = meshio.read(os.path.join(out, f\"bubble_{body}.msh\"))\n    tris = stl.cells_dict[\"triangle\"]\n    print(f\"body {body}: STL volume {enclosed_volume(stl.points, tris):.6f}; \"\n          f\".msh cells {[(c.type, len(c.data)) for c in msh.cells]}\")",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "c19",
   "metadata": {},
   "source": "The flat triangles enclose the target volumes to within 0.05% (the curved\nsurface is sampled at their corners); the\n`.msh` files hold quadratic triangles (`triangle6`) with the facet ids as\nGmsh physical and geometrical tags, ready for a volume mesher.\n\n## Next\n\n* [Relaxing a surface](relaxing.md) explains what `relax()` did here, how to read its\n  result, and where its defaults stop working.\n* The [liquid bridge](liquid_bridge.ipynb) adds walls: a drop between two beads, with\n  contact angles, forces and a volume mesh.\n* `ev.command(...)` runs any Evolver command; the [manual](https://kenbrakke.com/evolver/html/evolver.htm)\n  describes them all. `pyse.map` runs parameter sweeps in worker processes, one engine\n  each. The [API reference](api.md) documents every public class and function."
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}