Test function

Hello Everyone!
I am trying to solve a problem involving interface between two materials. My code is attached below:

geo = geom2d.SplineGeometry()

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

pp1 = geo.AppendPoint(0,0)

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

# f = 100*exp(-20**2*((x-0.0)*(x-0.0)+y*y))
# Draw(f,mesh, 'f')
Draw(mesh)
mesh.SetPML(pml.Cartesian((-2,-2),(-2,0),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 = 0
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).Trace()*grad(v).Trace()*ds("bc1")
a += rhoA*omega**2*u*v*ds("bc1")
a += -muB*grad(u).Trace()*grad(v).Trace()*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 have attached the weak form also as an image. I am having trouble in understanding how can we deal with terms containing v(0,0) (shown in the weak form). Any help would do. Thank you.


 

Hi,

you are talking about point-sources ?

You can actually do it via

f = LinearForm(fes)
f += (17*v)(0.5, 0.5)
f.Assemble()

Joachim