# Velocity gradient-velocity-pressure for the Stokes

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(4
x - 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) +

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()

The `grad` is used inconsistently, row-wise or column-wise.

Use `Grad` instead, then it is guaranteed to be the row-wise gradient.

The method looks like the hybridized version of MCS. Then you can use the space `HCurlDiv` for the (deviatoric part of the) stress tensor.

Joachim

Thank you so much, Joachim. It is working now. Yes, I am using the hybridized version of MCS. Thank you for your suggestions.

Best regards,
Lina