4_Benchmark_ HOSQandGPLS_genus_2-summer_school.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "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 genus two surface\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 genus 2 mesh composed of          triangles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d4ae8bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh_path=\"../meshes/genus_N=15632.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": [
    "\n",
    "def phi(x: np.ndarray):\n",
    "    \n",
    "    return 2*x[1]*(x[1]*x[1] - 3*x[0]*x[0])*(1 - x[2]*x[2])+(x[0]*x[0] + x[1]*x[1])**2-(9*x[2]*x[2] - 1)*(1 - x[2]*x[2]);\n",
    "\n",
    "\n",
    "def dphi(x: np.ndarray):\n",
    "    \n",
    "    return np.array([4*x[0]*(x[0]*x[0] + x[1]*x[1] + 3*x[1]*(x[2]*x[2] - 1)) ,\n",
    "        4*x[1]*(x[0]*x[0] + x[1]*x[1]) + 4*x[1]*x[1]*(1 - x[2]*x[2]) + 2*(3*x[0]*x[0] - x[1]*x[1])*(x[2]*x[2] - 1) ,\n",
    "        4*x[2]*(x[1]*(3*x[0]*x[0] - x[1]*x[1]) + 9*x[2]*x[2] - 5)])"
   ]
  },
  {
   "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 a genus two surface.\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 =-4*np.pi\n",
    "    t1 = time()\n",
    "    print(\"Relative error: \", abs( nints - exact_area)/abs(exact_area))\n",
    "    print (\"The main function takes:\",{(t1-t0)})\n",
    "    error=abs( nints - exact_area)/abs(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"
   ]
  }
 ],
 "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
}