- Thank you received: 0
Implementing von Neumann Boundary Condition
4 years 2 weeks ago - 4 years 2 weeks ago #3290
by rgs
Implementing von Neumann Boundary Condition was created by rgs
Hello,
I am a mechanical engineer, and quite a beginner in FEM. I hope this question is not too basic. I am trying to implement the following problem in NGSolve:
[tex]{\nabla}^2 u = 0 [/tex]
with the von Neumann boundary condition
[tex] \nabla u . \textbf{n} = \textbf{e} . \textbf{n} [/tex]
on one surface (FSI in my script), with all other surfaces having periodic or symmetric conditions. 'e' is a unit vector directed along either the x, y or z axis and 'n' is the normal vector of the element face.
I tried to implement this by starting with the Poisson tutorial. This is the code that I have so far:
I would be grateful if someone can help me with is. Thank you in advance!
I am a mechanical engineer, and quite a beginner in FEM. I hope this question is not too basic. I am trying to implement the following problem in NGSolve:
[tex]{\nabla}^2 u = 0 [/tex]
with the von Neumann boundary condition
[tex] \nabla u . \textbf{n} = \textbf{e} . \textbf{n} [/tex]
on one surface (FSI in my script), with all other surfaces having periodic or symmetric conditions. 'e' is a unit vector directed along either the x, y or z axis and 'n' is the normal vector of the element face.
I tried to implement this by starting with the Poisson tutorial. This is the code that I have so far:
Code:
from netgen.csg import *
sphere=Sphere(Pnt(0,0,0),0.55).bc('FSI')
pX=Plane(Pnt(0.5,0,0),Vec(1,0,0))
nX=Plane(Pnt(-0.5,0,0),Vec(-1,0,0))
pY=Plane(Pnt(0,0.5,0),Vec(0,1,0))
nY=Plane(Pnt(0,-0.5,0),Vec(0,-1,0))
pZ=Plane(Pnt(0,0,0.5),Vec(0,0,1))
nZ=Plane(Pnt(0,0,-0.5),Vec(0,0,-1))
Zielinski=sphere*pX*nX*pY*nY*pZ*nZ
geo=CSGeometry()
geo.Add(Zielinski.mat("domain"), bcmod=[(pX,"xplus"),(nX,"xminus"),(pY,"yplus"),(nY,"yminus"),(pZ,"zplus"),(nZ,"zminus")])
geo.PeriodicSurfaces(pX,nX)
geo.PeriodicSurfaces(pY,nY)
geo.PeriodicSurfaces(pZ,nZ)
from ngsolve.comp import Mesh
ngmesh=geo.GenerateMesh(maxh=0.05)
# Converting netgen mesh to ngsolve mesh
mesh = Mesh(ngmesh)
#mesh.GetBoundaries()
from ngsolve import *
ngsglobals.msg_level = 1
# H1-conforming finite element space
fes = Periodic(H1(mesh, order=1 ))
# define trial- and test-functions
u = fes.TrialFunction()
v = fes.TestFunction()
# the right hand side
f = LinearForm(fes)
f += 0 * v * dx
# the bilinear-form
a = BilinearForm(fes, symmetric=True)
a += grad(u) *grad(v)*dx
a.Assemble()
f.Assemble()
e = CoefficientFunction( (0,0,1) )
# the solution field
gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs(), inverse="sparsecholesky") * f.vec
I would be grateful if someone can help me with is. Thank you in advance!
Last edit: 4 years 2 weeks ago by rgs.
Time to create page: 0.110 seconds