-
Notifications
You must be signed in to change notification settings - Fork 4
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
Showing
12 changed files
with
207 additions
and
124 deletions.
There are no files selected for viewing
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
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -5,8 +5,28 @@ | |
| "id": "a08f0db0", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "# Embedded Trefftz-DG: Helmholtz\n", | ||
| "Standard polynomial Trefftz functions for the Helmholtz equation do not exist, to circumvent this problem, we weaken our condition in the Trefftz space. We introduce a projection $\\Pi$ that is yet to be defined and define the weak Trefftz space and the embedded weak Trefftz DG method:\n", | ||
| "# Embedded Trefftz-DG: Helmholtz" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "d6d31794", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "from ngsolve import *\n", | ||
| "from ngstrefftz import *\n", | ||
| "from netgen.occ import *\n", | ||
| "from ngsolve.webgui import Draw" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "4e554177", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We consider the Helmholtz equation with Robin boundary conditions\n", | ||
| "\n", | ||
| "$$\n", | ||
| "\\newcommand{\\Th}{{\\mathcal{T}_h}} \n", | ||
|
|
@@ -33,54 +53,114 @@ | |
| "\\newcommand{\\be}{\\mathbf{e}}\n", | ||
| "\\newcommand{\\bx}{{\\mathbf x}}\n", | ||
| "\\newcommand{\\inner}[1]{\\langle #1 \\rangle}\n", | ||
| "\\begin{align*}\n", | ||
| " \\begin{cases}\n", | ||
| " -\\Delta u - \\omega^2 u= 0 &\\text{ in } \\dom, \\\\\n", | ||
| " \\frac{\\partial u}{\\partial \\nx} + i u = g &\\text{ on } \\partial \\dom.\n", | ||
| " \\end{cases}\n", | ||
| "\\end{align*}\n", | ||
| "$$\n", | ||
| "\n", | ||
| "Standard *polynomial* Trefftz functions for the Helmholtz equation do not exist, to circumvent this problem, we weaken our condition in the Trefftz space. We introduce a projection $\\Pi$ that is yet to be defined and define the weak Trefftz space and the embedded weak Trefftz DG method:\n", | ||
| "\n", | ||
| "$$\n", | ||
| "\\begin{align}\n", | ||
| " \\text{Find }u_{\\IT}\\in \\IT^p(\\Th)&,~\\text{ s.t. }\n", | ||
| " a_h(u_{\\IT},v_{\\IT})=\\ell(v_{\\IT})\\qquad \\forall v_{\\IT}\\in \\IT^p(\\Th)\\quad \\text{ with } \\\\\n", | ||
| " \\IT^p(\\Th)&=\\{v\\in \\Vhp,\\ \\Pi \\calL v=0\\}. \\label{eq:weakTspace}\n", | ||
| " \\IT^p(\\Th)&=\\{v\\in V^k(\\mathcal T_h),\\ \\Pi \\calL v=0\\}. \\label{eq:weakTspace}\n", | ||
| " % \\\\ &\\IT^p(\\Th)=\\{v\\in L^2(\\dom) \\sst \\restr{v}{K}\\in\\IT(K),\\forall K\\in\\Th\\}\n", | ||
| "\\end{align}\n", | ||
| "$$\n", | ||
| "\n", | ||
| "This way, we can re-define the matrix $\\bW$ for the general case as \n", | ||
| "For the Helmholtz problem we choose $\\Pi:V^{k}(\\mathcal T_h)\\rightarrow V^{k-2}(\\mathcal T_h)$ the $L^2$ orthogonal projection. This way, we can define the matrix $\\bW$ as \n", | ||
| "\n", | ||
| "$$\n", | ||
| "\\begin{align} \\label{def:W3}\n", | ||
| " (\\bW)_{ij}&=\\inner{\\calL\\phi_j,\\tilde\\calL\\phi_i}_{0,h}. \n", | ||
| " (\\bW)_{ij}&=\\inner{\\calL\\phi_j,\\tilde\\phi_i}_{0,h}. \n", | ||
| "\\end{align}\n", | ||
| "$$\n", | ||
| "\n", | ||
| "with test functions $\\tilde\\phi_i\\in V^{k-2}(\\mathcal T_h)$ and $\\calL=-\\Delta u -\\omega^2 u$." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "4855289d", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "def TrefftzHelmholtzEmb(fes):\n", | ||
| " mesh = fes.mesh\n", | ||
| " k = fes.globalorder\n", | ||
| " u = fes.TrialFunction()\n", | ||
| " \n", | ||
| " Lap = lambda u : sum(Trace(u.Operator('hesse')))\n", | ||
| " fes2 = L2(mesh, order=order-2, dgjumps=True,complex=True)\n", | ||
| " v = fes2.TestFunction()\n", | ||
| " op = -u*v*dx - Lap(u)*v*dx\n", | ||
| " PP = TrefftzEmbedding(op,fes,test_fes=fes2)\n", | ||
| " return PP" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "263b8ff9", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "def EmbeddedBasisFunctionsAsMultiDim(Emb, fes):\n", | ||
| " gfshow = GridFunction(fes, multidim=0)\n", | ||
| " gf = GridFunction(fes)\n", | ||
| " coefvec = Emb.CreateRowVector()\n", | ||
| " for i in range(Emb.width):\n", | ||
| " coefvec[:] = 0\n", | ||
| " coefvec[i] = 1\n", | ||
| " gf.vec.data = Emb * coefvec\n", | ||
| " gfshow.AddMultiDimComponent(gf.vec)\n", | ||
| " return gfshow\n", | ||
| "\n", | ||
| "For the Helmholtz equation with Robin boundary conditions\n", | ||
| "\n", | ||
| "$$\n", | ||
| "\\begin{align*}\n", | ||
| " \\begin{cases}\n", | ||
| " -\\Delta u - \\omega^2 u= 0 &\\text{ in } \\dom, \\\\\n", | ||
| " \\frac{\\partial u}{\\partial \\nx} + i u = g &\\text{ on } \\partial \\dom.\n", | ||
| " \\end{cases}\n", | ||
| "\\end{align*}\n", | ||
| "$$\n", | ||
| "order = 3\n", | ||
| "mesh = Mesh(unit_square.GenerateMesh(maxh=1))\n", | ||
| "fes = L2(mesh, order=order, dgjumps=True,complex=True)\n", | ||
| "PP = TrefftzHelmholtzEmb(fes)\n", | ||
| "gfshow = EmbeddedBasisFunctionsAsMultiDim(PP,fes)\n", | ||
| "Draw (gfshow, mesh, interpolate_multidim=False, animate=False, autoscale=True, settings={\"subdivision\":20})" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "7ef6c524", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We compare this to the planewave space \n", | ||
| "\n", | ||
| "we choose the operators\n", | ||
| "\n", | ||
| "$$\n", | ||
| "\\begin{align*}\n", | ||
| "\\calL=-\\Delta u -\\omega^2 u && \\tilde\\calL=-\\Delta u\n", | ||
| " \\IT^p=\\{e^{-i\\omega(d_j\\cdot \\bx)},\\ j=-p,\\dots,p\\}.\n", | ||
| "\\end{align*}\n", | ||
| "$$" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "d6d31794", | ||
| "id": "9964431c", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "from ngsolve import *\n", | ||
| "from ngstrefftz import *\n", | ||
| "from netgen.geom2d import unit_square\n", | ||
| "from netgen.csg import unit_cube" | ||
| "order = 3\n", | ||
| "mesh = Mesh(unit_square.GenerateMesh(maxh=1))\n", | ||
| "fes = trefftzfespace(mesh, order=order, eq=\"helmholtz\",dgjumps=True,complex=True)\n", | ||
| "gfshow = GridFunction(fes, multidim=0)\n", | ||
| "gf = GridFunction(fes)\n", | ||
| "for i in range(fes.ndof):\n", | ||
| " gf.vec[:] = 0\n", | ||
| " gf.vec[i] = 1\n", | ||
| " gfshow.AddMultiDimComponent(gf.vec)\n", | ||
| "Draw (gfshow, mesh, interpolate_multidim=False, animate=False, autoscale=False, min=0.8,max=1, settings={\"subdivision\":20})" | ||
| ] | ||
| }, | ||
| { | ||
|
|
@@ -91,7 +171,8 @@ | |
| "outputs": [], | ||
| "source": [ | ||
| "mesh = Mesh(unit_square.GenerateMesh(maxh=.3))\n", | ||
| "fes = L2(mesh, order=4, complex=True, dgjumps=True)\n", | ||
| "order = 4\n", | ||
| "fes = L2(mesh, order=order, complex=True, dgjumps=True)\n", | ||
| "u,v = fes.TnT()" | ||
| ] | ||
| }, | ||
|
|
@@ -110,6 +191,24 @@ | |
| "eps = 10**-7" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "6111befd", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We consider the DG-scheme given by\n", | ||
| "\n", | ||
| "$$\n", | ||
| "\\begin{align}\n", | ||
| " a_h(u,v) &= \\sum_{K\\in\\Th}\\int_K \\nabla u\\nabla v-\\omega^2 uv\\ dV\n", | ||
| " -\\int_{\\Fh^\\text{int}}\\left(\\avg{\\nabla u}\\jump{v}+\\jump{u} \\avg{\\overline{\\nabla v}} \\right) dS \\nonumber \\\\ \n", | ||
| " &\\qquad+\\int_{\\Fh^\\text{int}} \\left( i\\alpha \\omega\\jump{u}\\jump{\\overline{v}} - \\frac{\\beta}{i\\omega}\\jump{\\nabla u}\\jump{\\overline{\\nabla v}} \\right) dS -\\int_{\\Fh^\\text{bnd}}\\delta\\left(\\nx\\cdot\\nabla u \\overline{v}+u \\overline{\\nx\\cdot\\nabla v}\\right) dS\\\\ \\nonumber\n", | ||
| " &\\qquad+\\int_{\\Fh^\\text{bnd}} \\left( i(1-\\delta)\\omega{u}{\\overline{v}} - \\frac{\\delta}{i\\omega}{\\nabla u}{\\overline{\\nabla v}} \\right) dS \\\\ \n", | ||
| " \\ell(v) &= \\int_{\\Fh^\\text{bnd}}\\left( (1-\\delta)g\\overline{v} - \\frac{\\delta}{i\\omega}g\\overline{\\nx\\cdot\\nabla v}\\right) dS\n", | ||
| "\\end{align}\n", | ||
| "$$" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
|
|
@@ -157,10 +256,7 @@ | |
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "Lap = lambda u : sum(Trace(u.Operator('hesse')))\n", | ||
| "op = (-u)*(Lap(v))*dx + (-Lap(u))*(Lap(v))*dx\n", | ||
| "eps = 10**-8\n", | ||
| "PP = TrefftzEmbedding(op,fes,eps)\n", | ||
| "PP = TrefftzHelmholtzEmb(fes)\n", | ||
| "PPT = PP.CreateTranspose()\n", | ||
| "with TaskManager():\n", | ||
| " TA = [email protected]@PP\n", | ||
|
|
||
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.