-
Notifications
You must be signed in to change notification settings - Fork 3
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
9dfdb12
commit bb911ff
Showing
30 changed files
with
3,504 additions
and
185 deletions.
There are no files selected for viewing
Binary file added
BIN
+11.7 KB
_images/62830f45f3668151e15c51437c689304b6404196d28410d15f53bb76ad5c9b67.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added
BIN
+11.8 KB
_images/9afaa03fb2f71b53130a264ce3e652f6d636eb575b7beaa9ca6e748272bd6a6a.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Large diffs are not rendered by default.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,234 @@ | ||
| { | ||
| "cells": [ | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "8fd56fc3", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "# Computing the gradient of a Gaussian peak in dual cell spaces\n", | ||
| "\n", | ||
| "### by M. Wess, 2024\n", | ||
| "*This Notebook is part of the `dualcellspaces` [documentation](https://ngsolve.github.io/dcm) for the addon package implementing the Dual Cell method in [NGSolve](https://ngsolve.org).*" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "0b129f3c", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We verify our implementation by projecting a Gaussian peak\n", | ||
| "\n", | ||
| "$$\n", | ||
| "f(\\mathbf x) = \\frac{1}{2}\\exp(-100 \\|\\mathbf x - (\\tfrac{1}{2},\\tfrac{1}{2})^\\top\\|^2),\n", | ||
| "$$\n", | ||
| "\n", | ||
| " in $\\mathbb R^2$ into the discrete space $\\tilde X^{\\mathrm{grad}}_P(\\tilde{\\mathcal T})$ and computing the discrete gradient in $X^{\\mathrm{div}}_P({\\mathcal T})$.\n" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "ee49f25a", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We import the packages and define the necessary spaces" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 1, | ||
| "id": "60f2ed24", | ||
| "metadata": {}, | ||
| "outputs": [ | ||
| { | ||
| "name": "stdout", | ||
| "output_type": "stream", | ||
| "text": [ | ||
| "DoFs H1Primal: 154870\n", | ||
| "DoFs HDivDual: 344250\n" | ||
| ] | ||
| } | ||
| ], | ||
| "source": [ | ||
| "from ngsolve import *\n", | ||
| "import dualcellspaces as dcs\n", | ||
| "from ngsolve.webgui import Draw\n", | ||
| "from time import time\n", | ||
| "\n", | ||
| "mesh = Mesh(unit_square.GenerateMesh(maxh = 0.03))\n", | ||
| "\n", | ||
| "order = 4\n", | ||
| "h1 = dcs.H1DualCells(mesh, order = order)\n", | ||
| "hdiv = dcs.HDivPrimalCells(mesh, order = order)\n", | ||
| "\n", | ||
| "print(\"DoFs H1Primal: {}\".format(h1.ndof))\n", | ||
| "print(\"DoFs HDivDual: {}\".format(hdiv.ndof))" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "97e8613e", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We obtain the (lumped) mass matrix by using `Mass`. We assemble the right hand side for the projection using the lumped integration rule and solve the projection problem to find $p\\in \\tilde X^{\\mathrm{grad}}_P(\\tilde {\\mathcal T})$\n", | ||
| "\n", | ||
| "$$\n", | ||
| "(p,q)_h = (f,q)_h,\n", | ||
| "$$\n", | ||
| "\n", | ||
| "where $(\\cdot,\\cdot)_h$ is the lumped approximation of the $L^2$ inner product." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 3, | ||
| "id": "5f43a2c0", | ||
| "metadata": {}, | ||
| "outputs": [ | ||
| { | ||
| "data": { | ||
| "application/vnd.jupyter.widget-view+json": { | ||
| "model_id": "7e38079058ea4fae9ce05e0b895c631d", | ||
| "version_major": 2, | ||
| "version_minor": 0 | ||
| }, | ||
| "text/plain": [ | ||
| "WebGuiWidget(layout=Layout(height='50vh', width='100%'), value={'gui_settings': {'camera': {'euler_angles': [-…" | ||
| ] | ||
| }, | ||
| "metadata": {}, | ||
| "output_type": "display_data" | ||
| } | ||
| ], | ||
| "source": [ | ||
| "mass_h1_inv = h1.Mass(1).Inverse()\n", | ||
| "\n", | ||
| "dx_h1 = dx(intrules = h1.GetIntegrationRules())\n", | ||
| "peak = CF( 0.5 * exp(-100*( (x-0.5)**2 + (y-0.5)**2 )) )\n", | ||
| "\n", | ||
| "p,q = h1.TnT()\n", | ||
| "rhs = LinearForm(peak*q*dx_h1).Assemble().vec\n", | ||
| "\n", | ||
| "gfp = GridFunction(h1)\n", | ||
| "gfp.vec.data = mass_h1_inv * rhs\n", | ||
| "\n", | ||
| "Draw(gfp, order = 2, points = dcs.GetWebGuiPoints(2), deformation = True, euler_angles = [-40,-4,-150]);" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "e495bd22", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "Next we need to assemble the differential operator (including the boundary terms)\n", | ||
| "$$\n", | ||
| "b(p,v) = \\sum_{T\\in\\mathcal T} -\\int_T p \\,\\mathrm{div} v dx +\\int_{\\partial T} p v\\cdot n ds.\n", | ||
| "$$" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 11, | ||
| "id": "7cb4ebf7", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "n = specialcf.normal(2)\n", | ||
| "\n", | ||
| "dSw = dx(element_boundary = True, intrules = dcs.GetIntegrationRules(2*order - 1))\n", | ||
| "dxw = dx(intrules = dcs.GetIntegrationRules(2*order -1))\n", | ||
| "\n", | ||
| "v = hdiv.TestFunction()\n", | ||
| "grad = BilinearForm(-p*div(v)*dxw + p*(v*n)*dSw, geom_free = True).Assemble().mat" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "974ab5bb", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "The flag `geom_free = True` tells the `BilinearForm` that the element contributions of the matrix are independent of specific element and the geometric quantities (i.e., all Jacobian matrices of the transformation cancel out)." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "e5a38e20", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "Lastly we solve the weak problem to find $u\\in X^{\\mathrm{div}}_P(\\mathcal T)$ such that\n", | ||
| "$$\n", | ||
| "(u,v)_h = b(p,v)\n", | ||
| "$$\n", | ||
| "for all $v$.\n" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 13, | ||
| "id": "e781f39d", | ||
| "metadata": {}, | ||
| "outputs": [ | ||
| { | ||
| "data": { | ||
| "application/vnd.jupyter.widget-view+json": { | ||
| "model_id": "cfe8e4d589ae4d228c0ed86cab6c41e8", | ||
| "version_major": 2, | ||
| "version_minor": 0 | ||
| }, | ||
| "text/plain": [ | ||
| "WebGuiWidget(layout=Layout(height='50vh', width='100%'), value={'gui_settings': {'camera': {'euler_angles': [-…" | ||
| ] | ||
| }, | ||
| "metadata": {}, | ||
| "output_type": "display_data" | ||
| } | ||
| ], | ||
| "source": [ | ||
| "gfu = GridFunction(hdiv)\n", | ||
| "\n", | ||
| "mass_hdiv_inv = hdiv.Mass().Inverse()\n", | ||
| "\n", | ||
| "gfu.vec.data = mass_hdiv_inv @ grad * gfp.vec\n", | ||
| "\n", | ||
| "Draw(gfu, order = 2, points = dcs.GetWebGuiPoints(2), vectors = True, euler_angles = [-40,-4,-150]);" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "4c548f6b", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "### Exercises\n", | ||
| "* Add/Remove the flag `geom_free = True` and study the application and setup times of the discrete differential operator." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "24c64bf9", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [] | ||
| } | ||
| ], | ||
| "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.10.12" | ||
| } | ||
| }, | ||
| "nbformat": 4, | ||
| "nbformat_minor": 5 | ||
| } |
Large diffs are not rendered by default.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.