Periodicity in 3D

General Discussion
1

Hello everyone, i hope you are doing welle, i juste have a probleme in periodicity for this stokes probleme,"load “msh3”
load “iovtk”
load “medit”

real pEps = 1.e-8;
real r=0.05;

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 = 10;

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));

int nz = 10;
real h = 1; // hauteur du cube

mesh3 Th3 = buildlayers(Th, nz, zbound=[-h/2, h/2]);

// Vérification visuelle
medit(“maillage3D”, Th3);

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

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

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

// Solve
Stokes;
p -= int2d(Th)(p) / int2d(Th)(1.0); // p ← p - moyenne(p)
p -= int2d(Th)(p) / int2d(Th)(1.0); // p ← p - moyenne(p)
p -= int2d(Th)(p) / int2d(Th)(1.0); // p ← p - moyenne(p)
p -= int2d(Th)(p) / int2d(Th)(1.0); // p ← p - moyenne(p)

// 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”);" it was a 2d probleme at first and now i need it 3D,but periodicity don’t work any more,thank you so much for your help

2

Dear Taha,
There is a problem of matching boundaries. This can be seen on the medit image. You see that the 2d boundaries 11 and 22 from the 3d mesh will not match when putting them face to face. The same problem may occur for boundaries 10,36.
I have changed the 2d numbering of vertices at the problematic locations to that it works, but this approach can only be done at hand, seeing the problem on the 3d mesh.
I hope this will be enough for you!
Apart from that, I have put label 51 for the “hole” boundaries,
label 101 for the bottom and 102 for the top.
3dperiodic.edp (3.1 KB)