How to define a partial differentiation on boundaries

Main Forum Old-Forum - do not post here
1

Hello Everyone!
I have written a code for solving a scattering problem. I have used PML boundaries (may not be properly used since I am still working on it) on three sides and on the top boundary I have some mixed boundary conditions. The code is like this:

import netgen.gui
from ngsolve import *
import netgen.geom2d as geom2d

geo = geom2d.SplineGeometry()

p1,p2,p3,p4,p5,p6,p7,p8 = [geo.AppendPoint(x,y) for x,y in [(0,0),(4,0),(4,2),(3.5,2),(3.5,0.5),(0.5,0.5),(0.5,2),(0,2)] ]

pp1 = geo.AppendPoint(1.5,2)

geo.Append ([“line”, p4, pp1],bc=“bc1”,leftdomain=1,rightdomain=0)
geo.Append ([“line”, pp1, p7],bc=“bc2”,leftdomain=1,rightdomain=0)

geo.Append ([“line”, p1, p2],leftdomain=2,rightdomain=0)
geo.Append ([“line”, p2, p3],leftdomain=2,rightdomain=0)
geo.Append ([“line”, p3, p4],leftdomain=2,rightdomain=0)
geo.Append ([“line”, p4, p5],leftdomain=2,rightdomain=1)
geo.Append ([“line”, p5, p6],leftdomain=2,rightdomain=1)
geo.Append ([“line”, p6, p7],leftdomain=2,rightdomain=1)
geo.Append ([“line”, p7, p8],leftdomain=2,rightdomain=0)
geo.Append ([“line”, p8, p1],leftdomain=2,rightdomain=0)

geo.SetMaterial(1, “inner”)
geo.SetMaterial(2, “pmlregion”)

#geo.AddRectangle ( (0,0), (4,2),leftdomain=1,rightdomain=0)
mesh=Mesh(geo.GenerateMesh (maxh=0.1))
Draw(mesh)

mesh.SetPML(pml.Cartesian((0,0),(0,2),1j),“pmlregion”)

fes = H1(mesh, order=5, complex=True)
u, v = fes.TnT()

mu = 0.1
muA = 0.01
muB = 0.02
rho = 1
rhoA = 2
rhoB = 1.1
omega = 10
kA = 10
gammaA = 1

f1 = -exp(gammaAy-1jkAx)(rhoomega**2+mugammaA2-mu*kA2)
f2 = -exp(-1jkAx)(gammaAmu+kA2muA-rhoAomega2)
f3 = -exp(-1jkAx)(gammaAmu+kA2muB-rhoBomega2)

a = BilinearForm(fes)
a += grad(u)grad(v)dx + rhoomega**2uvdx
a += -muAgrad(u)[0]grad(v)[0]ds(“bc1”)
a += rhoAomega**2uvds(“bc1”)
a += -muBgrad(u)[0]grad(v)[0]ds(“bc2”)
a += rhoBomega**2uvds(“bc2”)
a.Assemble()

f = LinearForm(fes)

f += f1vdx
f += f2vds(“bc1”)
f += f3vds(“bc2”)
f.Assemble()

gfu = GridFunction(fes, name=“u”)
gfu.vec.data = a.mat.Inverse() * f.vec
Draw(gfu)

I am getting this error: Trialfunction does not support BND-forms, maybe a Trace() operator is missing

How to deal with error?

Thank you.

2

Hi Gaurav,

you have to write

a += -muB*grad(u).Trace()*grad(v).Trace()*ds("bc2")

to get the surface-gradient, which is the projection to the tangential component of the full gradient.

It can be evaluated from the boundary values, only.

you can find some docu for it here:
https://docu.ngsolve.org/latest/i-tutorials/unit-6.1.2-surfacepde/surface_pdes.html#

Joachim

3

Thank you Joachim!