# Stokes in Cubic Pipe

**URL:** <https://community.freefem.org/t/stokes-in-cubic-pipe/2744>\
**Category:** General Discussion\
**Created:** [October 10, 2023, 12:15am UTC](https://community.freefem.org/t/stokes-in-cubic-pipe/2744 "2023-10-10T00:15:32Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![FluidMan](https://avatars.discourse-cdn.com/v4/letter/f/e68b1a/32.png) [@FluidMan](https://community.freefem.org/u/FluidMan)\
**Post date:** [October 10, 2023, 12:15am UTC](https://community.freefem.org/t/stokes-in-cubic-pipe/2744/1 "2023-10-10T00:15:32Z")

</div>

Hello everyone,

I have been using the basic features of FreeFem for years, but now I am facing a problem that I can’t seem to make sense of. Specifically, I’m struggling to implement boundary conditions for a 3D Stokes solver.

For context, I have a cube with inflow and outflow conditions on the left and right faces, and I want to impose no-slip and no-penetration conditions on the other boundaries.

I’ve found that the documentation on this topic is possibly erroneous, and I’m having difficulty specifying the conditions I need, particularly on the left and right faces (Note that it NEEDS to be the left and right faces and I need to use `buildlayers` because I plan to make the problem more complex later on). I have no problem specify condition on the top and bottom face using labels 1 and 2. The issues seems to be with specifying labelmid (where the documentation says the other surfaces are labeled).

Does have any insights or suggestions on how to specify these boundary conditions? I would greatly appreciate any help.

-FluidMan

```auto
load "msh3" load "medit"
int nn=5;
int Wallt = 11;	//Pipe wall label
int Wallb = 22;	//Pipe wall label
int Inlet = 33;	//Pipe inlet label
int Outlet = 44;	//Pipe outlet label
int Res = 5; //Resolution of Solver

// Define mesh boundaries
border C01(t=0, 1){x=0; y=t; label=Inlet;}
border C02(t=0, 1){x=1-t; y=0; label=Wallb;}
border C03(t=0, 1){x=1; y=-t+1; label=Outlet;}
border C04(t=0, 1){x=t; y=1; label=Wallt;}

mesh Th2= buildmesh(C01(-Res) + C02(-Res) + C03(-Res) + C04(-Res));

fespace Vh2(Th2,P2); Vh2 ux,uz,p2;
int[int] rup=[0,2], rdown=[0,1], rmid=[11,3,22,4,33,5,44,6];
real zmin=0,zmax=1;
mesh3 Th=buildlayers(Th2,nn,
  zbound=[zmin,zmax], labelmid=rmid, 
  reffaceup = rup, reffacelow = rdown);
  
medit("c10x10x10",Th,wait=1);

// FFCS: testing 3d plots
plot(Th);

fespace VVh(Th,[P2,P2,P2,P1]);
fespace Vh(Th,P23d);
macro Grad(u) [dx(u),dy(u),dz(u)]// EOM
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM

VVh [u1,u2,u3,p];
VVh [v1,v2,v3,q];
func fup = (1-x)*(x)*y*(1-y)*50;
solve vStokes([u1,u2,u3,p],[v1,v2,v3,q]) = 
  int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) + Grad(u2)'*Grad(v2) + Grad(u3)'*Grad(v3) //)';
		       - div(u1,u2,u3)*q - div(v1,v2,v3)*p + 1e-10*q*p ) 
  + on(2,u1=0,u2=0,u3=0) + on(1,u1=0,u2=0,u3=0) + on(LEFT FACE,u1=fup,u2=0,u3=0)+ on(RIGHT FACE,u1=fup,u2=0,u3=0)+ on(BACK FACE,u1=0,u2=0,u3=0)+ on(FRONT FACE,u1=0,u2=0,u3=0);

for(int i=1;i<10;i++)
  {
    real yy=i/10.;
    ux= u1(x,yy,y);
    uz= u3(x,yy,y);
    p2= p(x,yy,y);
    plot([ux,uz],p2,cmm=" cut y = "+yy,wait= 1);
  }

```

---

<div class="post-metadata">

**Author:** ![julienG](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/julieng/32/222_2.png) [@julienG](https://community.freefem.org/u/julienG)\
**Post date:** [October 10, 2023, 2:32am UTC](https://community.freefem.org/t/stokes-in-cubic-pipe/2744/2 "2023-10-10T02:32:00Z")

</div>

hmm, I’m not sure it does what you want, but maybe this will help.

You can `buildlayers` and then rotate with `movemesh`

so

```auto
load "msh3" load "medit"

real Lx = 5.0, Ly=Lx, Lz=2*Lx;
int NN=5;
int[int] Labels=[0,1,2,3];
mesh Th2 = square(NN,NN,[Lx*x,Ly*y],label=Labels);

int Inlet = 33;	//Pipe inlet label
int Outlet = 44;	//Pipe outlet label
int Faces = 66;	//Pipe outlet label

int[int] rup=[0,Outlet], rdown=[0,Inlet], rmid=[0,Faces,1,Faces,2,Faces,3,Faces];
real zmin=-10,zmax=10;
mesh3 Th=buildlayers(Th2,int(NN*Lz/Lx),
 zbound=[-0.5*Lz,0.5*Lz], labelmid=rmid, 
 reffaceup = rup, reffacelow = rdown);

plot(Th,wait=1);

real alpha=0.5*pi;
Th = movemesh(Th,[x*cos(alpha)-z*sin(alpha) ,y, x*sin(alpha)+z*cos(alpha)]);
plot(Th, wait=1);

fespace Vh2(Th2,P2); Vh2 ux,uz,p2;

fespace VVh(Th,[P2,P2,P2,P1]);
fespace Vh(Th,P23d);
macro Grad(u) [dx(u),dy(u),dz(u)]// EOM
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM

VVh [u1,u2,u3,p];
VVh [v1,v2,v3,q];
func fup = (Lx-x)*(x)*y*(Ly-y)*50;
solve vStokes([u1,u2,u3,p],[v1,v2,v3,q]) = 
 int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) + Grad(u2)'*Grad(v2) + Grad(u3)'*Grad(v3) //)';
   	       - div(u1,u2,u3)*q - div(v1,v2,v3)*p + 1e-10*q*p ) 
 + on(Faces,u1=0,u2=0,u3=0) + on(Inlet,u1=fup,u2=0,u3=0)+ on(Outlet ,u1=fup,u2=0,u3=0);

for(int i=1;i<10;i++)
 {
   real yy=i/10.;
   ux= u1(x,yy,y);
   uz= u3(x,yy,y);
   p2= p(x,yy,y);
   plot([ux,uz],p2,cmm=" cut y = "+yy,wait= 1);
 }

```

does this do as you want?

---

<div class="post-metadata">

**Author:** ![FluidMan](https://avatars.discourse-cdn.com/v4/letter/f/e68b1a/32.png) [@FluidMan](https://community.freefem.org/u/FluidMan)\
**Post date:** [October 10, 2023, 11:22am UTC](https://community.freefem.org/t/stokes-in-cubic-pipe/2744/3 "2023-10-10T11:22:37Z")

</div>

Thank you for the help. I am trying to work toward prescribing general Neumann conditions on all surfaces so I need to be able to specify velocity on all 6 surfaces of the cube. This code is closer but it still seems like conditions are initiable only on the bottom and top faces (then you rotated them to be the left and right faces).

---

<div class="post-metadata">

**Author:** ![julienG](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/julieng/32/222_2.png) [@julienG](https://community.freefem.org/u/julienG)\
**Post date:** [October 10, 2023, 4:22pm UTC](https://community.freefem.org/t/stokes-in-cubic-pipe/2744/4 "2023-10-10T16:22:44Z")

</div>

You can separate them also

```auto
load "msh3" load "medit"

real Lx = 5.0, Ly=Lx, Lz=2*Lx;
int NN=5;
int[int] Labels=[0,1,2,3];
mesh Th2 = square(NN,NN,[Lx*x,Ly*y],label=Labels);

int Inlet = 33;	//Pipe inlet label
int Outlet = 44;	//Pipe outlet label
int Face1 = 65;	//Pipe outlet label
int Face2 = 66;	//Pipe outlet label
int Face3 = 67;	//Pipe outlet label
int Face4 = 68;	//Pipe outlet label

int[int] rup=[0,Outlet], rdown=[0,Inlet], rmid=[0,Face1,1,Face2,2,Face3,3,Face3];
real zmin=-10,zmax=10;
mesh3 Th=buildlayers(Th2,int(NN*Lz/Lx),
 zbound=[-0.5*Lz,0.5*Lz], labelmid=rmid, 
 reffaceup = rup, reffacelow = rdown);

plot(Th,wait=1);

real alpha=0.5*pi;
Th = movemesh(Th,[x*cos(alpha)-z*sin(alpha) ,y, x*sin(alpha)+z*cos(alpha)]);
plot(Th, wait=1);

fespace Vh2(Th2,P2); Vh2 ux,uz,p2;

fespace VVh(Th,[P2,P2,P2,P1]);
fespace Vh(Th,P23d);
macro Grad(u) [dx(u),dy(u),dz(u)]// EOM
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM

VVh [u1,u2,u3,p];
VVh [v1,v2,v3,q];
func fup = (Lx-x)*(x)*y*(Ly-y)*50;
solve vStokes([u1,u2,u3,p],[v1,v2,v3,q]) = 
 int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) + Grad(u2)'*Grad(v2) + Grad(u3)'*Grad(v3) //)';
   	       - div(u1,u2,u3)*q - div(v1,v2,v3)*p + 1e-10*q*p ) 
 + on(Face1,u1=0,u2=0,u3=0) 
 + on(Face2,u1=0,u2=0,u3=0) 
 + on(Face3,u1=0,u2=0,u3=0) 
 + on(Face4,u1=0,u2=0,u3=0) 
+ on(Inlet,u1=fup,u2=0,u3=0)+ on(Outlet ,u1=fup,u2=0,u3=0);

for(int i=1;i<10;i++)
 {
   real yy=i/10.;
   ux= u1(x,yy,y);
   uz= u3(x,yy,y);
   p2= p(x,yy,y);
   plot([ux,uz],p2,cmm=" cut y = "+yy,wait= 1);
 }

```