5_Benchmark_ HOSQandGPLS_torus-summer-school.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "453a92e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Imports\n",
    "import numpy as np\n",
    "import scipy.io\n",
    "import matplotlib.pyplot as plt\n",
    "from math import pi\n",
    "from time import time\n",
    "import minterpy_levelsets as ls\n",
    "\n",
    "# Local imports\n",
    "import surfgeopy as sp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1acd9786",
   "metadata": {},
   "source": [
    "# Gauss Bonnet theorem on torus\n",
    "This is a benchmark of computing surface integrals using high-order surface quadrature (HOSQ) method and GPLS method for algebraic varieties.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68c55eab",
   "metadata": {},
   "source": [
    "# Step 1: Load and extract the vertices of a torus mesh composed of          triangles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d4ae8bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh_path=\"../meshes/torus_N=2512.mat\"\n",
    "mesh_mat = scipy.io.loadmat(mesh_path)\n",
    "pointcloud= mesh_mat[\"xs\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "506a928b",
   "metadata": {},
   "outputs": [],
   "source": [
    "R=0.5\n",
    "r=0.3\n",
    "def phi(x: np.ndarray):\n",
    "    ph = np.sqrt(x[0]*x[0] + x[1]*x[1])\n",
    "    return (ph - R)*(ph - R) + x[2]*x[2] - r*r\n",
    "def dphi(x: np.ndarray):\n",
    "    ph = np.sqrt(x[0]*x[0] + x[1]*x[1])\n",
    "    return np.array([-2*R*x[0]/ph + 2*x[0],-2*R*x[1]/ph + 2*x[1],2*x[2]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e280a58",
   "metadata": {},
   "source": [
    "# Step 2: Perform surface fitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a68b029f",
   "metadata": {},
   "outputs": [],
   "source": [
    "result = np.zeros(pointcloud.shape[0])  # initialize an empty array to store the results\n",
    "for i in range(pointcloud.shape[0]):\n",
    "    result[i] = phi(pointcloud[i])\n",
    "print(f\"The accuracy of the given mesh is:{(result).max()}\")\n",
    "if (result).max()>1e-8:\n",
    "    for i in range(pointcloud.shape[0]):\n",
    "        pointcloud[i,:]=sp.SimpleImplicitSurfaceProjection(phi, dphi,pointcloud[i,:])\n",
    "newt_poly = ls.LevelsetPoly(pointcloud, method='BK',verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ac0e8a4",
   "metadata": {},
   "source": [
    "# Step 3: We execute  HOSQ for the torus $\\mathbb{T}^2_{r,R}$.\n",
    "\n",
    "\n",
    "This step involves utilizing the zero level set provided by GPLS and computing the curvature with the help of GPLS:\n",
    "\n",
    "If $M$ is orienatble, then  the gradient $\\nabla Q_M = (\\partial_x Q_M,\\, \\partial_y Q_M,\\, \\partial_z Q_M) \\in \\mathbb{R}^3$ never vanishes on $M$ and provides, together with the Hessian $H_M =\\nabla(\\nabla Q_M) \\in \\mathbb{R}^{3 \\times 3}$ of $Q_M$\n",
    "\\begin{equation*}\n",
    " H_M  = \\left(\\begin{array}{ccc}\n",
    " \\frac{\\partial^2 Q_M}{\\partial_{x}^2} & \\frac{\\partial^2 Q_M}{\\partial_{x}\\partial_y} & \\frac{\\partial^2 Q_M}{\\partial_{x}\\partial_z} \\\\\n",
    " \\frac{\\partial^2 Q_M}{\\partial_{y}\\partial_x} & \\frac{\\partial^2 Q_M}{\\partial_y^2} & \\frac{\\partial^2 Q_M}{\\partial_{y}\\partial_z} \\\\\n",
    " \\frac{\\partial^2 Q_M}{\\partial_z\\partial_{x}} & \\frac{\\partial^2 Q_M}{\\partial_{z}\\partial_y} & \\frac{\\partial^2 Q_M}{\\partial_z^2} \\\\\n",
    " \\end{array}\\right) ,\n",
    "\\end{equation*}\n",
    "the main ingredients for the following computations. Both Gauss and mean curvature can be computed from these quantities \\cite{goldman2005} as:\n",
    "\\begin{align}\n",
    "  K_{\\mathrm{Gauss}} &= \\frac{\\det \\left(\\begin{array}{cc}H_M  & \\nabla Q_M^T \\\\  \\nabla Q_M & 0 \\end{array}\\right)}{\\|\\nabla Q_M\\|^4} \\label{eq:GC}\\\\\n",
    "  K_{\\mathrm{mean}} &= \\frac{\\nabla Q_M H_M \\nabla Q_M^T - \\|\\nabla Q_M\\|^2\\mathrm{trace}(H_M)}{2\\|\\nabla Q_M\\|^3}\\,. \\label{eq:MC}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b11488e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def err_t(intp_degree,lp_degree,func,grad_fun,mesh_path,Refinement):\n",
    "    mesh_mat = scipy.io.loadmat(mesh_path)\n",
    "    xs =mesh_mat[\"xs\"]\n",
    "    elems = mesh_mat[\"surfs\"] - 1\n",
    "    t0 = time()\n",
    "    pnts, ws, offset = sp.compute_surf_quadrature(func,grad_fun, xs, elems,intp_degree,lp_degree,Refinement)\n",
    "    # number of quadrature points\n",
    "    nqp = pnts.shape[0]\n",
    "    #compute the mean and Gauss curvature using GPLS\n",
    "    mean_curvature, gauss_curvature = newt_poly.compute_curvatures_at(pnts)\n",
    "    # function value at each quadrature points\n",
    "    fs_qp = np.array(\n",
    "            [\n",
    "                    gauss_curvature[qpid]\n",
    "\n",
    "                for qpid in range(nqp)\n",
    "            ]\n",
    "        )\n",
    "    # numerical integrations\n",
    "    nints = np.matmul(ws, fs_qp)\n",
    "    exact_area =0\n",
    "    t1 = time()\n",
    "    print(\"Relative error: \", abs( nints - exact_area))\n",
    "    print (\"The main function takes:\",{(t1-t0)})\n",
    "    error=abs( nints - exact_area)\n",
    "    return error\n",
    "\n",
    "# here is the error obtained using Dune Curved Grid\n",
    "error_dune2_14=np.array([9.262908e-05, 4.163018e-05, 6.594602e-08, 1.587536e-07,\n",
    "                         5.133794e-10, 9.538793e-10, 1.248652e-11, 9.704872e-11, \n",
    "                         1.370585e-09, 1.559943e-08, 1.181970e-07,3.655481e-07, 4.165369e-06])\n",
    "\n",
    "eror_HOSQ_GPLS=np.array([6.0068600790345785e-05,1.4799681049813568e-05,4.8092809326518576e-08,\n",
    "                        8.489890112109433e-09,1.9726209229443958e-10,7.845483377538764e-12,\n",
    "                        2.545225315231381e-13,1.8339431535896455e-13,6.941821249772939e-14,\n",
    "                       7.100700236128965e-14,7.192987525050931e-14,5.520822514426538e-14,5.989306273157524e-14])\n",
    "# Degree of Polynomial\n",
    "Nrange = list(range(2,15))\n",
    "lp_degree=float(\"inf\")\n",
    "refinement=0\n",
    "error1=[] \n",
    "for n in Nrange:\n",
    "    if n%1==0:print(n)\n",
    "    erro1 = err_t(int(n),lp_degree,newt_poly,newt_poly.compute_gradients_at,mesh_path,refinement)\n",
    "    error1.append(erro1)\n",
    "\n",
    "plt.semilogy(Nrange, eror_HOSQ_GPLS, '-or')\n",
    "plt.semilogy(Nrange, error_dune2_14, '-ob')\n",
    "plt.xlabel(\"Degree of Polynomial\",fontsize=13)\n",
    "plt.ylabel(\"Relative Error\",fontsize=13)\n",
    "plt.legend(['HOSQ+GPLS','DCG'],prop={'size': 13},loc='upper center')\n",
    "plt.xticks(np.arange(min(Nrange), max(Nrange)+1, 1.0))\n",
    "plt.ylim([2.758195177427762e-16,3.9514540203871754e-04])\n",
    "plt.grid()\n",
    "plt.savefig(\"../images/dune_vs_HOSQ+GPLS_Torus.pdf\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}