Hi, when I try to solve the generallized Stokes problem, I find some Strange phenomena： The pressure error in the cornor is so large. And I don’t know where I was wrong. Who can help me point out the mistakes? Thanks a lot.

The function is

`````` -2*mu*div(0.5*(grad(u)+grad(u)^T))+grad(ksi)=f
-div(u)+1/lambda*ksi=fp
``````

and the code is

``````int n = 64;

mesh Th =square(n,n);
//plot(Th,wait=1);

real E = 1000.0, nu = 0.499;
real lambda = E*nu/((1+nu)*(1-2*nu)), mu= E/(2*(1+nu)); //

real sqrt2=sqrt(2.);
macro epsilon(u1,u2) [dx(u1),dy(u2),(dy(u1)+dx(u2))/sqrt2] // EOM
macro div(u,v) ( dx(u)+dy(v) ) // EOM

// Test case 6
func uxtrue = sin(x);
func uytrue = sin(y);
func ksitrue=cos(x+y);

func fux=2*mu*sin(x)-sin(x+y);
func fuy=2*mu*sin(y)-sin(x+y);
func fp=-(cos(x)+cos(y)+1/lambda*cos(x+y));

fespace Vh(Th,[P2,P2,P1]);

// Define FE variables
Vh [ux,uy,ksi],[vx,vy,phi],[uxerr,uyerr,ksierr];

varf GStokes([ux,uy,ksi],[vx,vy,phi])
=int2d(Th)(2*mu*(epsilon(ux,uy)'*epsilon(vx,vy))-ksi*div(vx,vy))
+int2d(Th)(fux*vx+fuy*vy)  //rhs
+int2d(Th)(-div(ux,uy)*phi-1/lambda*ksi*phi)
+int2d(Th)(fp*phi) //rhs
+on(1,2,3,4,ux=uxtrue,uy=uytrue);

matrix MatB1=GStokes(Vh,Vh,solver=UMFPACK);
real[int] rhs1=GStokes(0,Vh);
ux[]=MatB1^-1*rhs1;

[uxerr,uyerr,ksierr]=[ux,uy,ksi]-[uxtrue,uytrue,ksitrue];
plot(ksi,value=1);
plot(ksierr,value=1);
``````

the Strange places is