Decoupling-Simulation-Accuracy-from-Mesh-Quality.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Decoupling Simulation Accuracy from Mesh Quality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simple python notebook to regenrate the teaser image.\n",
    "\n",
    "The data and scripts to regenerate the other figures can be found \n",
    "[here](https://github.com/polyfem/Decoupling-Simulation-Accuracy-from-Mesh-Quality).\n",
    "\n",
    "Note, we used [polyfem](https://polyfem.github.io) [C++ version](https://github.com/polyfem/polyfem/)  for all figures, with direct solver Pardiso, this notebook uses the default iterative algebraic multigrid Hypre."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this can be run direcly with binder!\n",
    "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/polyfem/Decoupling-Simulation-Accuracy-from-Mesh-Quality.git/master?filepath=Decoupling-Simulation-Accuracy-from-Mesh-Quality.ipynb)\n",
    "\n",
    "Importing libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import polyfempy as pf\n",
    "import json\n",
    "\n",
    "#plotting:\n",
    "import meshplot as mp\n",
    "import plotly.offline as plotly\n",
    "import plotly.graph_objs as go\n",
    "#Necessary for the notebook\n",
    "plotly.init_notebook_mode(connected=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this notebook is designed to be run in the root folder.\n",
    "If you want to reproduce the results of the teaser you can find the images [here](https://github.com/polyfem/Decoupling-Simulation-Accuracy-from-Mesh-Quality/tree/master/fig1-teaser/meshes)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setup problem, the exact definition of the function can be found [here](https://github.com/polyfem/polyfem/blob/d598e1764512da28d62ecbbbd1a832407ffb42de/src/problem/MiscProblem.cpp#L23)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "settings = pf.Settings()\n",
    "\n",
    "#Teaser uses laplacian\n",
    "settings.set_pde(pf.PDEs.Laplacian)\n",
    "problem = pf.Problem()\n",
    "#Use quartic function f = (2*y-0.9)^4+0.1 for bounary conditions\n",
    "#f''=\\Deta f = 48*(2*y-.9)^2 for rhs\n",
    "problem.add_dirichlet_value(1, \"(2*y-0.9)^4+0.1\")\n",
    "problem.exact = \"(2*y-0.9)^4+0.1\"\n",
    "problem.exact_grad = [\"0\", \"8*(2*y-.9)^3\"]\n",
    "problem.rhs = \"48*(2*y-.9)^2\"\n",
    "\n",
    "\n",
    "settings.problem = problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we solve using standard $P_1$ elements and our method which we enable with `settings.set_advanced_option(\"use_p_ref\", True)`.\n",
    "\n",
    "**Note** this will take some time, the last mesh quality is really low!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_meshes = 8\n",
    "\n",
    "total_solutions = {}\n",
    "total_errors = {}\n",
    "\n",
    "for use_pref in [False, True]:\n",
    "    settings.set_advanced_option(\"use_p_ref\", use_pref)\n",
    "    \n",
    "    solver = pf.Solver()\n",
    "    solver.set_log_level(6)\n",
    "    solver.settings(settings)\n",
    "\n",
    "    solutions = []\n",
    "    errors = []\n",
    "\n",
    "    #iterating over input files, named from 0 to 8\n",
    "    for i in range(n_meshes+1):\n",
    "        solver.load_mesh_from_path(\"fig1-teaser/meshes/large_angle_strip_{}.obj\".format(i), normalize_mesh=True)\n",
    "        \n",
    "        #we dont distinguis bc contitions, all sidesets to 1\n",
    "        solver.set_boundary_side_set_from_bary(lambda v: 1)\n",
    "        \n",
    "        solver.solve()\n",
    "        solver.compute_errors()\n",
    "\n",
    "        #getting and storing solution\n",
    "        v, f, sol = solver.get_sampled_solution()\n",
    "        pts = np.append(v, sol, 1)\n",
    "        solutions.append({\"pts\": pts, \"f\": f})\n",
    "\n",
    "        #getting error:\n",
    "        log = solver.get_log()\n",
    "        log = json.loads(log)\n",
    "        errors.append(log['err_l2'])\n",
    "        \n",
    "    #store these list in the total one\n",
    "    total_solutions[\"ours\" if use_pref else \"p1\"] = solutions\n",
    "    total_errors[\"ours\" if use_pref else \"p1\"] = errors\n",
    "\n",
    "print(\"Done!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plz wait until done!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot 3D results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For nice plot we append the 2 meshses, ours on the left"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_v_f_c(sols, i):\n",
    "    #append points and offset x by 1.1\n",
    "    pts_ours = np.array(total_solutions[\"ours\"][i][\"pts\"])\n",
    "    pts_p1 = np.array(total_solutions[\"p1\"][i][\"pts\"])\n",
    "    pts_p1[:, 0] += 1.1\n",
    "    v = np.append(pts_ours, pts_p1, 0)\n",
    "    \n",
    "    #colors\n",
    "    c = np.array(v[:, 2])\n",
    "    \n",
    "    #append faces\n",
    "    f_ours = np.array(total_solutions[\"ours\"][i]['f'])\n",
    "    f_p1 = np.array(total_solutions[\"p1\"][i]['f'])\n",
    "    f_p1 += np.max(f_ours)+1\n",
    "    f = np.append(f_ours, f_p1, 0)\n",
    "    \n",
    "    return v, f, c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our results are on the left!\n",
    "\n",
    "You can rotate the view with the mouse!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "v, f, c = get_v_f_c(total_solutions, 0)\n",
    "p = mp.plot(v, f, c, return_plot=True)\n",
    "\n",
    "@mp.interact(mesh=[('mesh_{}'.format(i), i) for i in range(n_meshes+1)])\n",
    "def ff(mesh):\n",
    "    v, f, c = get_v_f_c(total_solutions, mesh)\n",
    "    mp.plot(v, f, c, plot=p)\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error plor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = {'ours': 'rgb(85, 239, 196)', 'p1': 'rgb(255, 201, 26)'}\n",
    "marker_shape = 'circle'\n",
    "marker_size = 6\n",
    "\n",
    "def get_plot_trace(errors, method):\n",
    "    x = []\n",
    "    y = []\n",
    "    \n",
    "    for i in range(len(errors[method])):\n",
    "        x.append(i)\n",
    "        y.append(errors[method][i])\n",
    "        \n",
    "    x, y = zip(*sorted(zip(x, y))[:])\n",
    "    trace = go.Scatter(x=x, y=y,\n",
    "        mode='lines+markers',\n",
    "        name=method,\n",
    "        line=dict(color=colors[method]),\n",
    "        marker=dict(symbol=marker_shape, size=marker_size)\n",
    "    )\n",
    "    return trace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_data = [get_plot_trace(total_errors, \"ours\"), get_plot_trace(total_errors, \"p1\")]\n",
    "\n",
    "layout = go.Layout(\n",
    "    legend=dict(x=0.1, y=0.98),\n",
    "    xaxis=dict(\n",
    "        title=\"Mesh\",\n",
    "        nticks=5,\n",
    "        ),\n",
    "    yaxis=dict(\n",
    "        title=\"Error L2\",\n",
    "        range=[0.0, 0.3]\n",
    "    ),\n",
    "    hovermode='closest')\n",
    "\n",
    "fig = go.Figure(data=plot_data, layout=layout)\n",
    "plotly.iplot(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}