I use velocity-gradient-velocity-pressure formulation to solve the Stokes equations based on the pairs (P^k-RT^k-P^k), where it represents the polynomial order for the velocity gradient (discontinuous, continuity is enforced via Lagrange multiplier), velocity ( H(div)-conforming space) and tangential facet space.
I use the following code to do the convergence test and homogeneous Dirichlet boundary condition is used. But the solution is not so correct, it seems that some minor issues occur. I guess it is due to the arrangement of the velocity gradient. Could you please help for me? I am a beginner of NGsolve and am still on my way to study. Thanks a lot.
from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.125))
order=1
S = L2(mesh, order=order, dgjumps=True)
P = L2(mesh, order=order, dgjumps=True)
V = HDiv(mesh, order=order,dirichlet=“bottom|left|right|top”,RT=True)
Fu = TangentialFacetFESpace(mesh, order = order, dirichlet = “bottom|left|right|top”)
N = NumberSpace(mesh)
fes = FESpace([S,S,S,S,V,P,Fu,N])
sigmaxx,sigmaxy,sigmayx,sigmayy,u,p,uhat, p0= fes.TrialFunction()
tauxx,tauxy,tauyx,tauyy,v,q,vhat,q0= fes.TestFunction()
sigma = CoefficientFunction( (sigmaxx, sigmaxy,
sigmayx, sigmayy), dims=(2,2))
tau = CoefficientFunction( (tauxx, tauxy,
tauyx, tauyy), dims=(2,2))
def tang(vec):
return vec - (vec*n)*n
mu=1
source = CoefficientFunction((cos(x)cos(y) + 2mupisin(2piy)(6x
- 6x**2 + 2x2.*pi2 - 4x**3pi2 + 2*x4pi**2 - 1),
-sin(x)sin(y) - mu(4x - 2)(- 2pi2cos(2y*pi)*x2 +
2pi**2cos(2ypi)x + 3cos(2ypi) - 3)))
n = specialcf.normal(mesh.dim)
a = BilinearForm(fes, symmetric=True)
a += (InnerProduct((sigma)/mu,tau) + InnerProduct(grad(u),tau) +
InnerProduct(grad(v),sigma) -pdiv(v)-qdiv(u))*dx
dS = dx(element_boundary=True)
a += (InnerProduct(sigman,tang(vhat-v))+ InnerProduct(tang(uhat-u),taun))*dS
a += SymbolicBFI (p0*q + p * q0)
a.Assemble()
f = LinearForm(fes)
f += -sourcevdx
f.Assemble()