Stokes Equtions with pure Dirichlet boundary

Main Forum Old-Forum - do not post here
1

Hello everyone,

I am trying to solve Stokes equations on a unit square with pure Dirichlet boundary conditions. I know in this case I need to modify the bilinear form, otherwise the solution is unstable. I have taken a look to the tutorials but they all consider homogenous Neumann conditions in on of the boundaries, none of them is pure Dirichlet. I am trying something like this (adding a lagrange multiplier to the bilinear form to incorporate pure Dirichlet bcs)

from ngsolve import *
import netgen.gui

[code]from netgen.geom2d import SplineGeometry
geo = SplineGeometry()
geo.AddRectangle( (-1, -1), (1, 1), bcs = (“wall” , “outlet”, “wall”, “inlet”))
mesh = Mesh( geo.GenerateMesh(maxh=0.25))
mesh.Curve(3)
Draw (mesh)

V = VectorH1(mesh, order=2, dirichlet=“wall|inlet|outlet”)
Q = H1(mesh, order=1)
X = FESpace([V,Q])

u,p = X.TrialFunction()
v,q = X.TestFunction()

r = Q.TrialFunction()
s = Q.TestFunction()

a = BilinearForm(X)
a += SymbolicBFI(InnerProduct(grad(u),grad(v))+div(u)q+div(v)p)
a += SymbolicBFI(ps+qr)
a.Assemble()

gfu = GridFunction(X)

res = gfu.vec.CreateVector()
res.data = -a.mat * gfu.vec
inv = a.mat.Inverse(freedofs=X.FreeDofs())
gfu.vec.data += inv * res[/code]

But the results makes no sense… Anyone knows how to impose this kind of boundary condition?
Thanks in advance.

2

You implemented Dirchlet BC on the whole domain correctly. I am not sure what’s your intention behind adding the term a += SymbolicBFI(p*s+q*r).

However, of course you need to add a driving force, otherwise your solution is zero. Your code without your mysterious term and an example force:

[code]from netgen.geom2d import SplineGeometry
geo = SplineGeometry()
geo.AddRectangle( (-1, -1), (1, 1), bcs = (“wall” , “outlet”, “wall”, “inlet”))
mesh = Mesh( geo.GenerateMesh(maxh=0.25))
#mesh.Curve(3)
Draw (mesh)

V = VectorH1(mesh, order=2, dirichlet=“wall|inlet|outlet”)
Q = H1(mesh, order=1)
X = FESpace([V,Q])

u,p = X.TrialFunction()
v,q = X.TestFunction()

r = Q.TrialFunction()
s = Q.TestFunction()

a = BilinearForm(X)
a += SymbolicBFI(InnerProduct(grad(u),grad(v))+div(u)q+div(v)p)
#a += SymbolicBFI(ps+qr)
a.Assemble()

f = LinearForm(X)
f += SymbolicLFI(-x*v[1])
f.Assemble()
gfu = GridFunction(X)

res = gfu.vec.CreateVector()
res.data = f.vec -a.mat * gfu.vec
inv = a.mat.Inverse(freedofs=X.FreeDofs())
gfu.vec.data += inv * res

Draw(gfu.components[0], mesh, “vel”)

[/code]