Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Test function

More
2 years 2 months ago #4196 by GAURAV
Test function was created by GAURAV
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.


 
More
2 years 2 months ago #4204 by joachim
Replied by joachim on topic Test function
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
 
Time to create page: 0.151 seconds