examples/ex1b-pacman.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 1.B: Pac-Man\n",
    "### Sam Reynolds, 2023\n",
    "\n",
    "This is a continuation of Example 1-A. We will compute $L^2$ inner products\n",
    "and $H^1$ inner products of implicitly-defined functions on a punctured\n",
    "Pac-Man domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "\n",
    "current_dir = os.getcwd()\n",
    "parent_dir = os.path.abspath(os.path.join(current_dir, \"..\"))\n",
    "sys.path.append(parent_dir)\n",
    "\n",
    "import puncturedfem as pf\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# define quadrature schemes\n",
    "quad_dict = pf.get_quad_dict(n=64)\n",
    "\n",
    "# define vertices\n",
    "verts = []\n",
    "verts.append(pf.Vert(x=0.0, y=0.0))\n",
    "verts.append(pf.Vert(x=(np.sqrt(3) / 2), y=0.5))\n",
    "verts.append(pf.Vert(x=(np.sqrt(3) / 2), y=-0.5))\n",
    "verts.append(pf.Vert(x=-0.1, y=0.5))\n",
    "\n",
    "# define edges\n",
    "edges = []\n",
    "edges.append(pf.Edge(verts[0], verts[1], pos_cell_idx=0))\n",
    "edges.append(\n",
    "    pf.Edge(\n",
    "        verts[1],\n",
    "        verts[2],\n",
    "        pos_cell_idx=0,\n",
    "        curve_type=\"circular_arc_deg\",\n",
    "        theta0=300,\n",
    "    )\n",
    ")\n",
    "edges.append(pf.Edge(verts[2], verts[0], pos_cell_idx=0))\n",
    "edges.append(\n",
    "    pf.Edge(\n",
    "        verts[3],\n",
    "        verts[3],\n",
    "        neg_cell_idx=0,\n",
    "        curve_type=\"circle\",\n",
    "        radius=0.25,\n",
    "        quad_type=\"trap\",\n",
    "    )\n",
    ")\n",
    "\n",
    "# define mesh cell\n",
    "K = pf.MeshCell(idx=0, edges=edges)\n",
    "\n",
    "# parameterize edges\n",
    "K.parameterize(quad_dict)\n",
    "\n",
    "# plot boundary\n",
    "pf.plot.MeshPlot(K.get_edges()).draw()\n",
    "\n",
    "# set up Nyström solver\n",
    "nyst = pf.NystromSolver(K, verbose=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Function with a gradient singularity\n",
    "\n",
    "Consider the function\n",
    "\\begin{align*}\n",
    "\tv(x) = r^\\alpha \\, \\sin(\\alpha \\theta)\n",
    "\\end{align*}\n",
    "where $x \\mapsto (r, \\theta)$ is given in polar coordinates,\n",
    "and $\\alpha = 1/2$ is a fixed parameter. \n",
    "Note that for $\\alpha < 1$, \n",
    "the gradient $\\nabla v$ has a singularity at the origin.\n",
    "However, $v$ is harmonic everywhere else. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get Cartesian coordinates of points on boundary\n",
    "x1, x2 = K.get_boundary_points()\n",
    "\n",
    "# convert to polar\n",
    "r = np.sqrt(x1**2 + x2**2)\n",
    "th = np.arctan2(x2, x1) % (2 * np.pi)\n",
    "\n",
    "# Dirichlet trace of v\n",
    "alpha = 1 / 2\n",
    "v_trace = r**alpha * np.sin(alpha * th)\n",
    "\n",
    "# Laplacian of v (harmonic function)\n",
    "v_lap = pf.Polynomial()\n",
    "\n",
    "# build local function\n",
    "v = pf.LocalFunction(nyst=nyst, lap_poly=v_lap, has_poly_trace=False)\n",
    "v.set_trace_values(v_trace)\n",
    "\n",
    "# compute all quantities needed for integration\n",
    "v.compute_all()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the normal derivative is unbounded near the origin.\n",
    "Let's take a look at the weighted normal derivative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pf.plot.TracePlot(\n",
    "    traces=v.harm_part_wnd,\n",
    "    title=\"Weighted normal derivative\",\n",
    "    fmt=\"k.\",\n",
    "    K=K,\n",
    "    quad_dict=quad_dict,\n",
    ").draw()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $H^1$ seminorm\n",
    "\n",
    "Let's try computing the square $H^1$ seminorm of $v$, \n",
    "\\begin{align*}\n",
    "\t\\int_K |\\nabla v|^2~dx\n",
    "\t&\\approx 1.20953682240855912\n",
    "\t\\pm 2.3929 \\times 10^{-18}\n",
    "\\end{align*}\n",
    "with an approximate value obtained with *Mathematica*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h1_norm_sq_computed = v.get_h1_semi_inner_prod(v)\n",
    "print(\"Computed square H^1 seminorm = \", h1_norm_sq_computed)\n",
    "\n",
    "h1_norm_sq_exact = 1.20953682240855912\n",
    "h1_norm_sq_error = abs(h1_norm_sq_computed - h1_norm_sq_exact)\n",
    "print(\"Error in square H^1 seminorm = %.4e\" % (h1_norm_sq_error))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $L^2$ norm\n",
    "\n",
    "Let's also try computing the $L^2$ norm\n",
    "\\begin{align*}\n",
    "\t\\int_K v^2 ~dx\n",
    "\t&\\approx 0.97793431492143971\n",
    "\t\\pm 3.6199\\times 10^{-19}\n",
    "\t~.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l2_norm_sq_computed = v.get_l2_inner_prod(v)\n",
    "print(\"Computed square L^2 seminorm = \", l2_norm_sq_computed)\n",
    "\n",
    "l2_norm_sq_exact = 0.977934314921439713\n",
    "l2_norm_sq_error = abs(l2_norm_sq_computed - l2_norm_sq_exact)\n",
    "print(\"Error in square L^2 seminorm = %.4e\" % l2_norm_sq_error)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convergence\n",
    "\n",
    "Using Kress parameter $p=7$\n",
    "\n",
    "|\tn\t|\tH1 error\t|\tL2 error\t|\n",
    "|-------|---------------|---------------|\n",
    "|\t4\t|\t7.2078e-02\t|\t2.1955e-02\t|\n",
    "|\t8\t|\t3.3022e-02\t|\t5.4798e-03\t|\n",
    "|\t16\t|\t1.2495e-03\t|\t1.0159e-04\t|\n",
    "|\t32\t|\t6.5683e-06\t|\t4.6050e-07\t|\n",
    "|\t64\t|\t4.6834e-08\t|\t2.1726e-09\t|"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus: Interior values\n",
    "\n",
    "Not included in the paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y1 = K.int_x1\n",
    "y2 = K.int_x2\n",
    "\n",
    "v.compute_interior_values()\n",
    "\n",
    "v_computed = v.int_vals\n",
    "v_x1_computed = v.int_grad1\n",
    "v_x2_computed = v.int_grad2\n",
    "\n",
    "plt.figure()\n",
    "plt.contourf(y1, y2, v_computed, levels=50)\n",
    "plt.colorbar()\n",
    "plt.title(\"Interior values of $v$\")\n",
    "\n",
    "plt.figure()\n",
    "plt.contourf(y1, y2, v_x1_computed, levels=50)\n",
    "plt.colorbar()\n",
    "plt.title(\"First component of grad $v$\")\n",
    "\n",
    "plt.figure()\n",
    "plt.contourf(y1, y2, v_x2_computed, levels=50)\n",
    "plt.colorbar()\n",
    "plt.title(\"Second component of grad $v$\")\n",
    "\n",
    "# convert to polar\n",
    "r = np.sqrt(y1**2 + y2**2)\n",
    "th = np.arctan2(y2, y1) % (2 * np.pi)\n",
    "\n",
    "# exact values\n",
    "cos_th = np.cos(th)\n",
    "sin_th = np.sin(th)\n",
    "cos_ath = np.cos(alpha * th)\n",
    "sin_ath = np.sin(alpha * th)\n",
    "\n",
    "v_exact = r**alpha * sin_ath\n",
    "v_x1_exact = alpha * r ** (alpha - 1) * (cos_th * sin_ath - sin_th * cos_ath)\n",
    "v_x2_exact = alpha * r ** (alpha - 1) * (sin_th * sin_ath + cos_th * cos_ath)\n",
    "\n",
    "# interior value errors\n",
    "v_error = np.log10(np.abs(v_computed - v_exact))\n",
    "plt.figure()\n",
    "plt.contourf(y1, y2, v_error, levels=50)\n",
    "plt.colorbar()\n",
    "plt.title(\"Interior errors ($\\log_{10}$)\")\n",
    "\n",
    "# first component of gradient errors\n",
    "v_x1_error = np.log10(np.abs(v_x1_computed - v_x1_exact))\n",
    "plt.figure()\n",
    "plt.contourf(y1, y2, v_x1_error, levels=50)\n",
    "plt.colorbar()\n",
    "plt.title(\"Gradient errors in $x_1$ ($\\log_{10}$)\")\n",
    "\n",
    "# second component of gradient errors\n",
    "v_x2_error = np.log10(np.abs(v_x2_computed - v_x2_exact))\n",
    "plt.figure()\n",
    "plt.contourf(y1, y2, v_x2_error, levels=50)\n",
    "plt.colorbar()\n",
    "plt.title(\"Gradient errors in $x_2$ ($\\log_{10}$)\")\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.11.6"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}