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


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


f = LinearForm(fes)

f += -sourcevdx


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.


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,