# How to implementthe penalty FEM for solving the steady stokes equations

**URL:** <https://community.freefem.org/t/how-to-implementthe-penalty-fem-for-solving-the-steady-stokes-equations/3236>\
**Category:** General Discussion\
**Created:** [May 11, 2024, 3:12am UTC](https://community.freefem.org/t/how-to-implementthe-penalty-fem-for-solving-the-steady-stokes-equations/3236 "2024-05-11T03:12:51Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![hanweiwei](https://avatars.discourse-cdn.com/v4/letter/h/bbe5ce/32.png) [@hanweiwei](https://community.freefem.org/u/hanweiwei)\
**Post date:** [May 11, 2024, 3:12am UTC](https://community.freefem.org/t/how-to-implementthe-penalty-fem-for-solving-the-steady-stokes-equations/3236/1 "2024-05-11T03:12:51Z")

</div>

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](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/5/55c05076df5e1bf48eb83c28deb60970e30513ee.png)

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [May 11, 2024, 2:36pm UTC](https://community.freefem.org/t/how-to-implementthe-penalty-fem-for-solving-the-steady-stokes-equations/3236/2 "2024-05-11T14:36:14Z")

</div>

I think that the only proper way to implement the projection term is to use the original coupled formulation  
 ![stokesfe](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/4/4d005d0d13e37ea8a0b457f910c63505ce03291e.jpeg)  
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](https://perso.math.u-pem.fr/bouchut.francois/splitmesh6.cpp)

---

<div class="post-metadata">

**Author:** ![hanweiwei](https://avatars.discourse-cdn.com/v4/letter/h/bbe5ce/32.png) [@hanweiwei](https://community.freefem.org/u/hanweiwei)\
**Post date:** [May 11, 2024, 11:55pm UTC](https://community.freefem.org/t/how-to-implementthe-penalty-fem-for-solving-the-steady-stokes-equations/3236/3 "2024-05-11T23:55:16Z")

</div>

> [@fb77](#):
>
> lumped projection,

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.