The problem of solving pressure(p) in 3D

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:
image
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;
	}  
}

I think I have make it