{ "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 }