Skip to content

Add plots, remove extra changes #20

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 8 commits into
base: master
Choose a base branch
from
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
100 changes: 99 additions & 1 deletion chapter1/poisson.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -257,7 +257,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "d829e410",
"id": "827c3e69-77ac-4312-a4dd-a1c26a40b27c",
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that this id can be changed by hand to its original value using a text editor (or from the GitHub web interface). If it causes a problem, our CI build will spot it.

"metadata": {},
"outputs": [],
"source": [
Expand All @@ -273,6 +273,104 @@
"# print the result\n",
"print(h1_error)"
]
},
{
"cell_type": "markdown",
"id": "14d92a71-8c0f-4a91-8f56-c405f3fc76d1",
"metadata": {},
"source": [
"### Visualization\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2a33eab4-f4b8-428d-ae35-2d086d5645c3",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
"from sympy import lambdify\n",
"\n",
"def refine_array_1d(breaks, N):\n",
" \"\"\"Refine a 1D array by inserting N points between each pair of original values.\"\"\"\n",
" refined = np.concatenate([np.linspace(breaks[i], breaks[i+1], N+1)[:-1] for i in range(len(breaks)-1)])\n",
" return np.append(refined, breaks[-1]) # Add the last point\n",
"\n",
"def plot_solutions(Vh, uh, phi_exact, N=10):\n",
" \"\"\"\n",
" Refine the grid, evaluate solutions, compute errors, and plot results.\n",
"\n",
" Parameters:\n",
" Vh: The finite element space (must have `spaces` attribute with `breaks`).\n",
" uh: The numerical solution function.\n",
" phi_exact: The exact solution function (callable).\n",
" N: Number of points to insert between breaks (default: 10).\n",
" \"\"\"\n",
" # Generate refined grid for visualization\n",
" eta1 = refine_array_1d(Vh.spaces[0].breaks, N)\n",
" eta2 = refine_array_1d(Vh.spaces[1].breaks, N)\n",
"\n",
" # Evaluate numerical solution on the refined grid\n",
" num = np.array([[uh(e1, e2) for e2 in eta2] for e1 in eta1])\n",
"\n",
" # Compute exact solution\n",
" ex = np.array([[phi_exact(e1, e2) for e2 in eta2] for e1 in eta1])\n",
" err = num - ex\n",
"\n",
" # Generate physical grid coordinates\n",
" xx, yy = np.meshgrid(eta1, eta2, indexing='ij')\n",
"\n",
" # Create figure with 3 subplots\n",
" fig, axes = plt.subplots(1, 3, figsize=(12.8, 4.8))\n",
"\n",
" def add_colorbar(im, ax):\n",
" \"\"\"Add a colorbar to the given axis.\"\"\"\n",
" divider = make_axes_locatable(ax)\n",
" cax = divider.append_axes(\"right\", size=0.2, pad=0.2)\n",
" cbar = ax.get_figure().colorbar(im, cax=cax)\n",
" return cbar\n",
"\n",
" # Plot exact solution\n",
" ax = axes[0]\n",
" im = ax.contourf(xx, yy, ex, 40, cmap='jet')\n",
" add_colorbar(im, ax)\n",
" ax.set_xlabel(r'$x$', rotation='horizontal')\n",
" ax.set_ylabel(r'$y$', rotation='horizontal')\n",
" ax.set_title(r'$\\phi_{exact}(x,y)$')\n",
" ax.set_aspect('equal')\n",
"\n",
" # Plot numerical solution\n",
" ax = axes[1]\n",
" im = ax.contourf(xx, yy, num, 40, cmap='jet')\n",
" add_colorbar(im, ax)\n",
" ax.set_xlabel(r'$x$', rotation='horizontal')\n",
" ax.set_ylabel(r'$y$', rotation='horizontal')\n",
" ax.set_title(r'$\\phi(x,y)$')\n",
" ax.set_aspect('equal')\n",
"\n",
" # Plot numerical error\n",
" ax = axes[2]\n",
" im = ax.contourf(xx, yy, err, 40, cmap='jet')\n",
" add_colorbar(im, ax)\n",
" ax.set_xlabel(r'$x$', rotation='horizontal')\n",
" ax.set_ylabel(r'$y$', rotation='horizontal')\n",
" ax.set_title(r'$\\phi(x,y) - \\phi_{exact}(x,y)$')\n",
" ax.set_aspect('equal')\n",
"\n",
" # Show figure\n",
" plt.tight_layout()\n",
" fig.show()\n",
"\n",
"# Define the exact solution as a symbolic expression (example: phi_exact = x^2 + y^2)\n",
"ue = sin(np.pi*x)*sin(np.pi*y)\n",
"\n",
Comment on lines +367 to +369
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ue was already defined in the section titled "Computing the L^2 norm", as it was needed to define the error norms. Please remove its definition from this place, and move the comment line (which is useful) to the previous definition.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

By the way, I have just realized that the same section redefines u, which was previously defined in the section titled "Formal Model". Could you please remove the second definition?

"# Convert the symbolic expression to a callable function using lambdify\n",
"phi_exact = lambdify((x, y), ue, 'numpy')\n",
"plot_solutions(Vh, uh, phi_exact, N=10)"
]
}
],
"metadata": {},
Expand Down