[Aliciemath]

Hi all,
I have written uzawa scheme for Bingham stokes flow on freefem. It is supposed to change when the bn number changes. However when I compile with different bn numbers it does not change. Could anyone help me?
Here is the code
real Bn = 10.;real Mu=1.;

real r = 1.; real r2=5e-6;
real L = 1; //Pipe length
real D = 1.; //Pipe height
real L1=.1;
real L2=.9;
real rho=1.; real alpha=1.;
macro Gradient(u) [dx(u), dy(u)] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
macro UgradV(ux,uy,vx,vy) [ [ux,uy]'*[dx(vx),dy(vx)] , [ux,uy]'*[dx(vy),dy(vy)] ]//
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))] //

func f1=1.;
func f2=1.;
border b1(t=0., 1.){x=t; y=0.;};
border b2(t=0., 1.){x=L; y=D*t; };
border b3(t=L, L2){x=t; y=D;};
border b4(t=L2, L1){x=t; y=D;};
border b5(t=L1, 0.){x=t; y=D;};
border b6(t=0., 1.){x=0.; y=D-D*t;};
int nnL = 50;
int nnD = 50;
mesh Th=buildmesh(b1(nnL)+b2(nnD)+b3(5)+b4(5)+b5(5)+b6( nnD));
//Th = adaptmesh(Th, 1./128., IsMetric=1, nbvx=10000);
fespace Wh(Th,P2);
fespace Qh(Th,P1);
fespace Xh(Th,P2);

Wh vx,vy,ux1,uy1,dux,duy;
Wh ux,uy;
Qh p, q, dp,p1;

fespace Lh(Th,P1);
Lh TmGxx, TmGyy, TmGxy;
Lh Txx, Tyy, Txy;
Lh Gdotxx, Gdotyy, Gdotxy;
Lh dTxx=0., dTyy=0., dTxy=0.;
Lh TmGxx1, TmGyy1, TmGxy1;
Lh TmGxx0, TmGyy0, TmGxy0;
func u5= 1.;
func u4=1.;
func u3=1.;


TmGxx=0.;
TmGyy=0.;

TmGxy=0.;

for(int j=0; j<100; ++j)
{
solve bing([ux,uy],[vx,vy])
=int2d(Th)((Mu)*(Gradient(ux)'*Gradient(vx)
+ Gradient(uy)'*Gradient(vy)
)
)

+int2d(Th)( Bn*[TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)
)
-int2d(Th)([f1 , f2]' *[vx ,vy]
)

+ on(b1, ux=0., uy=0.)
+ on(b2, ux=0., uy=0.)
+ on(b3, ux=u3, uy=0.)
+ on(b4, ux=u4, uy=0.)
+ on(b5, ux=u5, uy=0.);

Gdotxx = 2.*dx(ux);
Gdotyy = 2.*dy(uy);
Gdotxy = dx(uy)+dy(ux);

for(int k=0; k<Lh.ndof; ++k){
real TaGxx=TmGxx[][k]+r2*Gdotxx[][k],
TaGyy=TmGyy[][k]+r2*Gdotyy[][k],
TaGxy=TmGxy[][k]+r2*Gdotxy[][k];
real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+2.*TaGxy^2) );
real res = TaGnorm;
if( 1.<res ){

real q = 1./res;

TmGxx[][k] = TaGxx*q;
TmGyy[][k] = TaGyy*q;
TmGxy[][k] = TaGxy*q;
}
else {

TmGxx[][k] = TaGxx;
TmGyy[][k] = TaGyy;
TmGxy[][k] = TaGxy;
}
}ux0[]=ux[];
uy0[]=uy[];



}
fespace Rh(Th, P2);
Rh dyg,dyg2;
solve streamlines (dyg, dyg2,solver=UMFPACK)
= int2d(Th)(
dx(dyg)*dx(dyg2)
+ dy(dyg)*dy(dyg2)
)
+ int2d(Th)(
- dyg2*(dy(ux) - dx(uy))
)
+ on(1,2,3,4,5,6, dyg=0);
plot(dyg,ps="ldsre1500.eps",nbiso=10);
Here is the uzawa bingham scheme
laplace u^n+1+divergence lamda^n=f
lamdaa^n+1=P(lamda^n+D(n^n+1))
P(q)=q if |q|<1 p(q)=q/|q| if |q|>1