Hi,
I want to solve pressure with the initial velocity in 3d, but it doesn’t work, the variational equation comes from the following N-S equation:
the program as follows:
//
load "msh3"
load "tetgen"
load "medit"
//verbosity=2;
real hs4=0.02;
real voltet4=(hs4^3)/6.;
real[int] domain = [0.,0.,0.5,53,voltet4];
real x0=0.,x1=1;
real y0=0.,y1=1;
real z0=0.,z1=1;
int [int,int] L=[[1,2],[3,4],[5,6]];
int[int] refx=[0,L(0,0)],refX=[0,L(0,1)];
int[int] refy=[0,L(1,0)],refY=[0,L(1,1)];
int[int] refz=[0,L(2,0)],refZ=[0,L(2,1)];
mesh Thx=square(50,50,[y0+(y1-y0)*x,z0+(z1-z0)*y]);
mesh Thy=square(50,50,[x0+(x1-x0)*x,z0+(z1-z0)*y]);
mesh Thz=square(50,50,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
mesh3 Thx0=movemesh23(Thx,transfo=[x0,x,y],label=refx);
mesh3 Thx1=movemesh23(Thx,transfo=[x1,x,y],label=refX);
mesh3 Thy0=movemesh23(Thy,transfo=[x,y0,y],label=refy);
mesh3 Thy1=movemesh23(Thy,transfo=[x,y1,y],label=refY);
mesh3 Thz0=movemesh23(Thz,transfo=[x,y,z0],label=refz);
mesh3 Thz1=movemesh23(Thz,transfo=[x,y,z1],label=refZ);
mesh3 ThH=Thx0+Thx1+Thy0+Thy1+Thz0+Thz1;
plot(ThH,wait=1);
mesh3 Th=tetg(ThH,switch="paAAYYQ",nbofregions=1,regionlist=domain);
savemesh(Th,"mesh.mesh");
medit("mesh",wait=1,Th);
// FE Space
fespace Xh(Th,P23d); //for velocity
fespace Mh(Th,P13d); //for pressure and stress
Xh u1old,u2old,u3old;
Xh p,q;
// Physical parameter
real re=1000,;
int k;
for(re=1000;re<=1000;re=re+1) // Loop on vicosity
{
real nu= 1./re;
// stop test for Newton
real eps=1e-6;
// intial guess with B.C.
u1old=(z==1);
u2old=0;
u3old=0;
verbosity=0;
solve Pold(p,q) =
int3d(Th) (dx(p)*dx(q)+dy(p)*dy(q)+dz(p)*dz(q))
-int3d(Th) (q*(dx(u1old)*dx(u1old)+dy(u2old)*dy(u2old)+dz(u3old)*dz(u3old)+2*dy(u1old)*dx(u2old)+2*dz(u1old)*dx(u3old)+2*dz(u2old)*dy(u3old)));
ofstream ppp(re+"+p.txt");
for (int m=0;m<Th.nv;m++)
{
ppp<<Th(m).x<<" "<<Th(m).y<< " "<<Th(m).z<< " "<<p(Th(m).x,Th(m).y,Th(m).z)<<endl;
}
}