Hello everyone,

When using the penalty method for solving Stokes equation, we can decouple the computation of velocity field with the pressure. But how to implement projection operator?

I think that the only proper way to implement the projection term is to use the original coupled formulation

Otherwise a simplified way could be to use a lumped projection, but I do not know how to do that.

Nevertheless a particular case which is clean is if the spaces X_h, Q_h are such that \mathop{\rm div}X_h\subset Q_h, because then I_h is identity and you can directly solve

(\nabla u_{\varepsilon h},\nabla v_h)+\frac{1}{\varepsilon}(\nabla\cdot u_{\varepsilon h},\nabla\cdot v_h)=(f,v_h),\quad\forall v_h\in X_h.

Not that the space Q_h disappears then.

The couple of spaces X_h, Q_h need of course to satisfy an \inf\sup condition. This leads to so called Scott-Vogelius type finite elements.

Two simple cases exist:

- X_h=P2, Q_h=P1dc, with a Hsieh–Clough–Tocher mesh. Such a mesh Th can be obtained by a subdivision from an arbitrary mesh Th0 by the command

load “splitmesh3”

mesh Th=splitmesh3(Th0); - X_h=P1, Q_h=P_0, with a Powell-Sabin mesh. Such a mesh Th can be obtained by a subdivision from an arbitrary mesh Th0 by the command

load “splitmesh6”

mesh Th=splitmesh6PowellSabin(Th0);

This procedure is however only available in the developer branch of FreeFem++, but you can get it with the plugin

https://perso.math.u-pem.fr/bouchut.francois/splitmesh6.cpp

Dear professor,

Thanks for your answer. The case about div X_h \subset Q_h is interesting and I never consider this case. And i will try it.