How to define a partial differentiation on boundaries

More
2 years 11 months ago #4109 by GAURAV
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(gammaA*y-1j*kA*x)*(rho*omega**2+mu*gammaA**2-mu*kA**2)
f2 = -exp(-1j*kA*x)*(gammaA*mu+kA**2*muA-rhoA*omega**2)
f3 = -exp(-1j*kA*x)*(gammaA*mu+kA**2*muB-rhoB*omega**2)

a = BilinearForm(fes)
a += grad(u)*grad(v)*dx + rho*omega**2*u*v*dx
a += -muA*grad(u)[0]*grad(v)[0]*ds("bc1")
a += rhoA*omega**2*u*v*ds("bc1")
a += -muB*grad(u)[0]*grad(v)[0]*ds("bc2")
a += rhoB*omega**2*u*v*ds("bc2")
a.Assemble()

f = LinearForm(fes)

f += f1*v*dx
f += f2*v*ds("bc1")
f += f3*v*ds("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.
More
2 years 11 months ago #4111 by joachim
Hi Gaurav,

you have to write
Code:
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:
docu.ngsolve.org/latest/i-tutorials/unit...e/surface_pdes.html#

Joachim
More
2 years 11 months ago #4113 by GAURAV
Thank you Joachim!
Time to create page: 0.099 seconds