How to implementthe penalty FEM for solving the steady stokes equations

General Discussion
1

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?

penaltymethod
penaltymethod1196×941 62.5 KB

2

I think that the only proper way to implement the projection term is to use the original coupled formulation
stokesfe
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:

  1. 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);
  2. 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
3

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.