/usr/share/doc/freefem++/examples/examples++-3d/NSI3d-carac.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 | load "msh3"
real nu=0.01,dt=0.3;
real alpha=1./dt,alpha2=sqrt(alpha);
int nn=5;
mesh Th2=square(nn,nn);
fespace Vh2(Th2,P2);
Vh2 ux,uz,p2;
int[int] rup=[0,2], rdown=[0,1], rmid=[1,1,2,1,3,1,4,1];
real zmin=0,zmax=1;
mesh3 Th=buildlayers(Th2,nn,
zbound=[zmin,zmax],
// region=r1,
labelmid=rmid,
reffaceup = rup,
reffacelow = rdown);
fespace VVh(Th,[P23d,P23d,P23d,P13d]);
fespace Vh(Th,P23d);
fespace Ph(Th,P13d);
macro Grad(u) [dx(u),dy(u),dz(u)]// EOM
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM
varf 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=1.,u2=0,u3= 0)
+ on(1,u1=0,u2=0,u3=0)
;
cout << "b mat " << endl;
matrix A=vStokes(VVh,VVh);
cout << "e mat " << endl;
set(A,solver=UMFPACK);
cout << "e fac mat " << endl;
real[int] b= vStokes(0,VVh);
VVh [u1,u2,u3,p];
VVh [X1,X2,X3,Xp];
VVh [x1,x2,x3,xp]=[x,y,z,0];
u1[]= A^-1 * b;
ux= u1(x,0.5,y);
uz= u3(x,0.5,y);
p2= p(x,0.5,y);
plot([ux,uz],p2,cmm=" cut y = 0.5",wait=1);
macro XX1() (x-u1*dt)//
macro XX2() (y-u2*dt)//
macro XX3() (z-u3*dt)//
varf vNS([uu1,uu2,uu3,p],[v1,v2,v3,q]) =
int3d(Th)( alpha*(uu1*v1+uu2*v2+uu3*v3) + nu*(Grad(uu1)'*Grad(v1) + Grad(uu2)'*Grad(v2) + Grad(uu3)'*Grad(v3))
- div(uu1,uu2,uu3)*q - div(v1,v2,v3)*p + 1e-10*q*p )
+ on(2,uu1=1,uu2=0,uu3=0)
+ on(1,uu1=0,uu2=0,uu3=0)
+ int3d(Th,optimize=1,qforder=4)( alpha*( convect([u1,u2,u3],-dt,u1)*v1 + convect([u1,u2,u3],-dt,u2)*v2 + convect([u1,u2,u3],-dt,u3)*v3 ) ) ;
// + int3d(Th,optimize=1)( alpha*( u1(X1,X2,X3)*v1 + u2(X1,X2,X3)*v2 + u3(X1,X2,X3)*v3 ) ) ;
// + int3d(Th,optimize=1)( alpha*( u1(XX1,XX2,XX3)*v1 + u2(XX1,XX2,XX3)*v2 + u3(XX1,XX2,XX3)*v3 ) ) ;
//+ int3d(Th,optimize=1)( alpha*( u1(x,y,z)*v1 + u2(x,y,z)*v2 + u3(x,y,z)*v3 ) ) ;
//+ int3d(Th,optimize=1)( alpha*( u1*v1 + u2*v2 + u3*v3 ) ) ;
cout << " build A" << endl;
A = vNS(VVh,VVh);
cout << " fac A" << endl;
set(A,solver=UMFPACK);
real t=0;
for(int i=0;i<10;++i)
{
t += dt;
cout << " iteration " << i << " t = " << t << endl;
X1[]=x1[]+u1[]*(-dt);
// verbosity=1000;
b=vNS(0,VVh);
verbosity=2;
u1[]= A^-1 * b;
ux= u1(x,0.5,y);
uz= u3(x,0.5,y);
p2= p(x,0.5,y);
plot([ux,uz],p2,cmm=" cut y = 0.5, time ="+t,wait=0);
if(i%5==6)
{
exec("mkdir dd");
string prefu="dd/u-"+(100+i);
string prefp="dd/p-"+(100+i);
savemesh(Th,prefu+".mesh");
savemesh(Th,prefp+".mesh");
ofstream file(prefu+".bb");
ofstream filep(prefp+".bb");
Ph up1=u1,up2=u2,up3=u3,pp=p;
file << "3 1 3 "<< up1[].n << " 2 \n";
filep << "3 1 1 "<< pp[].n << " 2 \n";
for (int j=0;j<up1[].n ; j++)
{
file << up1[][j] <<" " <<up2[][j] <<" "<< up3[][j] <<"\n";
filep << pp[][j] << endl;
}
}
}
plot([ux,uz],p2,cmm=" cut y = 0.5, time ="+t,wait=1);
|