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.

How to define a partial differentiation on boundaries

More
2 years 3 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 3 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 3 months ago #4113 by GAURAV
Thank you Joachim!
Time to create page: 0.126 seconds