- Thank you received: 0
How to define a partial differentiation on boundaries
2 years 11 months ago #4109
by GAURAV
How to define a partial differentiation on boundaries was created 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.
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.
2 years 11 months ago #4111
by joachim
Replied by joachim on topic How to define a partial differentiation on boundaries
Hi Gaurav,
you have to write
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
you have to write
Code:
a += -muB*grad(u).Trace()*grad(v).Trace()*ds("bc2")
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
Time to create page: 0.099 seconds