Add Delaunay demo

demo/delaunay.ipynb

@@ -0,0 +1,205 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "d2f53310",
+   "metadata": {},
+   "source": [
+    "Here we'll show how to compute Delaunay triangulations in 2D.\n",
+    "Under the hood, computing the Delaunay triangulation of a 2D point set is equivalent to computing the 3D convex hull of those points lifted onto a paraboloid in 3-space.\n",
+    "This means that if you understand how convex hulls work, you basically understand how Delaunay triangulations work -- all the moving parts are the same, down to the visibility graph.\n",
+    "First, we'll generate some random input data."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "961f0a95",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "rng = np.random.default_rng(seed=1729)\n",
+    "num_points = 40\n",
+    "X = rng.normal(size=(num_points, 2))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "485820e4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import matplotlib.pyplot as plt\n",
+    "fig, axes = plt.subplots()\n",
+    "axes.set_aspect(\"equal\")\n",
+    "axes.scatter(X[:, 0], X[:, 1]);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bc863ca5",
+   "metadata": {},
+   "source": [
+    "The plot below shows what these points look like when lifted to a 3D paraboloid."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "675880b3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from mpl_toolkits import mplot3d\n",
+    "W = np.sum(X**2, axis=1)\n",
+    "\n",
+    "fig, ax = plt.subplots(subplot_kw={\"projection\": \"3d\"})\n",
+    "ax.scatter(*np.column_stack((X, W)).T);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "632d6212",
+   "metadata": {},
+   "source": [
+    "Much like for convex hulls, we'll use a state machine object that we'll call `delaunay_machine` to keep track of the progress of the algorithm."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "282842f3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import zmsh\n",
+    "delaunay_machine = zmsh.DelaunayMachine(X)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0f974c60",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from copy import deepcopy\n",
+    "geometries = [deepcopy(delaunay_machine.geometry)]\n",
+    "\n",
+    "while not delaunay_machine.is_done():\n",
+    "    delaunay_machine.step()\n",
+    "    geometries.append(deepcopy(delaunay_machine.geometry))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bc9af4ed",
+   "metadata": {},
+   "source": [
+    "There is only one extra step for Delaunay triangulations.\n",
+    "If we repurpose an existing algorithm to compute the convex hull of the points lifted up to a parabola, we're going to get two \"sides\" -- a top and a bottom.\n",
+    "We're only interested in the facets on the bottom of the parabola, so to get the desired output we need to filter out anything on top.\n",
+    "The code below does the filtering for us."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "49e3f281",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def filter_bottom_facets(geometry):\n",
+    "    topology = geometry.topology\n",
+    "    dimension = topology.dimension\n",
+    "    cells = topology.cells(dimension)\n",
+    "    cell_ids_to_remove = []\n",
+    "    for cell_id in range(len(cells)):\n",
+    "        faces_ids, matrices = cells.closure(cell_id)\n",
+    "        if len(faces_ids[0]) > 0:\n",
+    "            orientation = zmsh.simplicial.orientation(matrices)\n",
+    "            x = geometry.points[faces_ids[0]]\n",
+    "            if orientation * zmsh.predicates.volume(*x) >= 0:\n",
+    "                cell_ids_to_remove.append(cell_id)\n",
+    "\n",
+    "    D = topology.boundary(dimension)\n",
+    "    for cell_id in cell_ids_to_remove:\n",
+    "        D[:, cell_id] = 0\n",
+    "        \n",
+    "    for k in range(dimension - 1, 0, -1):\n",
+    "        cocells = topology.cocells(k)\n",
+    "        cell_ids_to_remove = []\n",
+    "        for cell_id in range(len(cocells)):\n",
+    "            if len(cocells[cell_id][0]) == 0:\n",
+    "                cell_ids_to_remove.append(cell_id)\n",
+    "                \n",
+    "        D = topology.boundary(k)\n",
+    "        for cell_id in cell_ids_to_remove:\n",
+    "            D[:, cell_id] = 0"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2843e08d-7252-4d6c-b4e2-8c8186b577a5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for geometry in geometries:\n",
+    "    filter_bottom_facets(geometry)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d79146eb",
+   "metadata": {},
+   "source": [
+    "Now we can see the progress of the algorithm at each step.\n",
+    "Some of the steps are adding facets to the top of the hull of the paraboloid; we'll see those in the animation below as steps that don't appear to make any progress."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9174934e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from ipywidgets import interact\n",
+    "@interact(step=(0, len(geometries) - 1))\n",
+    "def f(step=0):\n",
+    "    geometry = geometries[step]\n",
+    "\n",
+    "    fig, ax = plt.subplots()\n",
+    "    ax.set_aspect(\"equal\")\n",
+    "\n",
+    "    zmsh.visualize(geometry, dimension=1, ax=ax)\n",
+    "    zmsh.visualize(geometry, dimension=0, ax=ax)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "zmsh",
+   "language": "python",
+   "name": "zmsh"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}