# Average condition for pressur in stokes equation

**URL:** <https://community.freefem.org/t/average-condition-for-pressur-in-stokes-equation/3870>\
**Category:** General Discussion\
**Created:** [April 18, 2025, 2:18pm UTC](https://community.freefem.org/t/average-condition-for-pressur-in-stokes-equation/3870 "2025-04-18T14:18:36Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![taha](https://avatars.discourse-cdn.com/v4/letter/t/ea5d25/32.png) [@taha](https://community.freefem.org/u/taha)\
**Post date:** [April 18, 2025, 2:18pm UTC](https://community.freefem.org/t/average-condition-for-pressur-in-stokes-equation/3870/1 "2025-04-18T14:18:36Z")

</div>

Hello everyone,i have a small question is it possible to add a condition that says the average of pressure that we want to calculate should equal to 0 ,here is my code that calculate the stokes equation in porous medium  
"load “msh3”  
load “iovtk”

real Mu = 1.;   
real pEps = 1.e-10;   
real r=0.1;

border C01(t=-r, r){x=-1; y=t; label=10;}  
border C02(t=-1, -r){x=t; y=r; label=1;}  
border C03(t=-1, -r){x=t; y=-r; label=2;}  
border C04(t=r, 1){x=-r; y=t; label=3;}  
border C05(t=-1, -r){x=-r; y=t; label=4;}  
border C06(t=-r, r){x=t; y=1; label=11;}  
border C07(t=-r, r){x=t; y=-1; label=22;}  
border C08(t=-1, -r){x=r; y=t; label=5;}  
border C09(t=r, 1){x=r; y=t; label=6;}  
border C10(t=r, 1){x=t; y=r; label=7;}  
border C11(t=r, 1){x=t; y=-r; label=8;}  
border C12(t=-r, r){x=1; y=t; label=36;}

int n = 100;

mesh Th = buildmesh(C01(-n)+C02(-n)+C04(-n)+ C06(-n)  
+C09(n)+C10(-n)+C12(n)+C11(n)+C08(n)+ C07(n)+C05(-n) + C03(n));

border C1(t=-1, 1){x=-1; y=t; label=44;}  
border C2(t=-1, 1){x=t; y=1; label=45;}  
border C3(t=1, -1){x=1; y=t; label=46;}  
border C4(t=1, -1){x=t; y=-1; label=47;}

mesh Th1 = buildmesh(C1(-n) + C2(-n) + C3(-n) + C4(-n));

plot(Th, wait=true);

fespace All(Th, [P2, P2, P1], periodic=[[10,y],[36,y],[11,x],[22,x]]);  
All [Ux, Uy, p];  
All [Uhx, Uhy, ph];

// Macros remain identical  
macro grad(u) [dx(u), dy(u)] //  
macro Grad(U) [grad(U#x), grad(U#y)] //  
macro div(ux, uy) (dx(ux) + dy(uy)) //  
macro Div(U) div(U#x, U#y) //

// Modified problem definition  
problem Stokes ([Ux, Uy, p], [Uhx, Uhy, ph])  
= int2d(Th)(  
Mu \* (Grad(U) : Grad(Uh))  
- p \* Div(Uh)  
- ph \* Div(U)  
- pEps \* p \* ph  
)  
+ int2d(Th)( 1 \* Uhx + 0\* Uhy)  
+ on(1,2,3,4,5,6,7,8, Ux=0, Uy=0)  
;

// Solve  
Stokes;

// Plot  
plot(Ux, fill=true, value=true, wait=true, cmm=“vitesse Ux”);  
plot(Uy, fill=true, value=true, wait=true, cmm=“Vitesse Uy”);  
plot(p, fill=true, value=true, wait=true, cmm=“Pression p”);

real area = int2d(Th1)(1.0);

real Kxx = -int2d(Th)(Ux) / area;  
real Kyy = -int2d(Th)(Uy) / area;

cout \<\< "Kxx = " \<\< Kxx \<\< endl;  
cout \<\< "Kyy = " \<\< Kyy \<\< endl;  
real ara = int2d(Th)(1.0);  
real pmoy = int2d(Th)(p) / ara;  
cout \<\< "Moyenne de p = " \<\< pmoy \<\< endl;

fespace Vh1(Th, P1);  
Vh1 Ux1 = Ux, Uy1 = Uy, p1 = p;

```
int[int] order = [1, 1, 1]; 
savevtk("StokesSolution.vtu", Th, Ux1, Uy1, p1, 
        dataname="Velocity_X Velocity_Y Pressure", order=order);

```

" and thank you so much for your help

---

<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:** [April 18, 2025, 5:33pm UTC](https://community.freefem.org/t/average-condition-for-pressur-in-stokes-equation/3870/2 "2025-04-18T17:33:10Z")

</div>

Hello,  
Your formulation indeed sets the average pressure to zero.  
To see that, take in the variational formulation Uhx=0, Uhy=0, ph=1.  
It gives \int\_\Omega \mathop{\rm div} u+\epsilon\int\_\Omega p=0.  
You have \int\_\Omega \mathop{\rm div} u=\int\_{\partial \Omega}u\cdot N=0 since u=0 on \partial\Omega, thus \int\_\Omega p=0.

---

<div class="post-metadata">

**Author:** ![taha](https://avatars.discourse-cdn.com/v4/letter/t/ea5d25/32.png) [@taha](https://community.freefem.org/u/taha)\
**Post date:** [April 19, 2025, 4:40pm UTC](https://community.freefem.org/t/average-condition-for-pressur-in-stokes-equation/3870/3 "2025-04-19T16:40:51Z")

</div>

> [@fb77](#):
>
> Hello,  
> Your formulation indeed sets the average pressure to zero.

oh ok thank you so much, but when i calculate the average in the end it always give something li 10^-5 or something like but never 0,is it normal to have value near to 0 but never 0 ??

---

<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:** [April 21, 2025, 10:22am UTC](https://community.freefem.org/t/average-condition-for-pressur-in-stokes-equation/3870/4 "2025-04-21T10:22:26Z")

</div>

The fact is that there are always rounding errors in computations. Here you see that  
one should have \int p=-\frac{1}{\epsilon}\int \mathop{\rm div}u. The integral of \mathop{\rm div} u is computed, say with an error of 1.e-16. Then when you divide by \epsilon you get something much bigger. It shows that you should not take \epsilon too small. A good value is \epsilon=1.e-8 to take the half of 1.e-16.  
Another way is, once p is computed, to subtract to it its average.

---

<div class="post-metadata">

**Author:** ![taha](https://avatars.discourse-cdn.com/v4/letter/t/ea5d25/32.png) [@taha](https://community.freefem.org/u/taha)\
**Post date:** [April 21, 2025, 8:03pm UTC](https://community.freefem.org/t/average-condition-for-pressur-in-stokes-equation/3870/5 "2025-04-21T20:03:19Z")

</div>

i got it ,thank you for your help