Ctrl+K

Computing the gradient of a Gaussian peak in dual cell spaces

3.3. Computing the gradient of a Gaussian peak in dual cell spaces#

by M. Wess, 2024#

This Notebook is part of the dualcellspaces documentation for the addon package implementing the Dual Cell method in NGSolve.

We verify our implementation by projecting a Gaussian peak

\[ f(\mathbf x) = \frac{1}{2}\exp(-100 \|\mathbf x - (\tfrac{1}{2},\tfrac{1}{2})^\top\|^2), \]

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})\).

We import the packages and define the necessary spaces

from ngsolve import *
import dualcellspaces as dcs
from ngsolve.webgui import Draw
from time import time

mesh = Mesh(unit_square.GenerateMesh(maxh = 0.03))

order = 4
h1 = dcs.H1DualCells(mesh, order = order)
hdiv = dcs.HDivPrimalCells(mesh, order = order)

print("DoFs H1Primal: {}".format(h1.ndof))
print("DoFs HDivDual: {}".format(hdiv.ndof))

DoFs H1Primal: 154870
DoFs HDivDual: 344250

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})\)

\[ (p,q)_h = (f,q)_h, \]

where \((\cdot,\cdot)_h\) is the lumped approximation of the \(L^2\) inner product.

mass_h1_inv = h1.Mass(1).Inverse()

dx_h1 = dx(intrules = h1.GetIntegrationRules())
peak = CF( 0.5 * exp(-100*( (x-0.5)**2 + (y-0.5)**2 ))  )

p,q = h1.TnT()
rhs = LinearForm(peak*q*dx_h1).Assemble().vec

gfp = GridFunction(h1)
gfp.vec.data = mass_h1_inv * rhs

Draw(gfp, order = 2, points = dcs.GetWebGuiPoints(2), deformation = True, euler_angles = [-40,-4,-150]);

Next we need to assemble the differential operator (including the boundary terms) $\( b(p,v) = \sum_{T\in\mathcal T} -\int_T p \,\mathrm{div} v dx +\int_{\partial T} p v\cdot n ds. \)$

n = specialcf.normal(2)

dSw = dx(element_boundary = True, intrules = dcs.GetIntegrationRules(2*order - 1))
dxw = dx(intrules = dcs.GetIntegrationRules(2*order -1))

v = hdiv.TestFunction()
grad = BilinearForm(-p*div(v)*dxw + p*(v*n)*dSw, geom_free = True).Assemble().mat

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).

Lastly we solve the weak problem to find \(u\in X^{\mathrm{div}}_P(\mathcal T)\) such that $\( (u,v)_h = b(p,v) \)\( for all \)v$.

gfu = GridFunction(hdiv)

mass_hdiv_inv = hdiv.Mass().Inverse()

gfu.vec.data = mass_hdiv_inv @ grad * gfp.vec

Draw(gfu, order = 2, points = dcs.GetWebGuiPoints(2), vectors = True, euler_angles = [-40,-4,-150]);

Exercises#

  • Add/Remove the flag geom_free = True and study the application and setup times of the discrete differential operator.

Contents