getGu

 calc G_u
 r=residual for numjac, rest as usual

0001 function Gu=getGu(p,u,lam,r) 
0002 % calc G_u
0003 % r=residual for numjac, rest as usual
0004 global pj lamj; % for call of numjac
0005 bc=p.bcf(p,u,lam); del=1e-6;
0006 if(p.jsw<2) % analytically
0007    [c,fu,flam,b]=p.jac(p,u,lam); zerov=zeros(p.neq,1); 
0008    [Gu,Glam]=assempde(bc,p.points,p.edges,p.tria,c,-fu,zerov); 
0009    if(any(b)) 
0010        Kadv=assemadv(p.points,p.tria,b); 
0011        Gu=Gu-Kadv; 
0012    end 
0013 else % numerically
0014   [S,F]=assempde(bc,p.points,p.edges,p.tria,0,ones(p.neq*p.neq,1),0,ones(p.np*p.neq,1));
0015   thresh=del*ones(size(u)); pj=p; lamj=lam; % global variable to pass to resinj
0016   [Gu,fac,G]=numjac('resinj',0,u,r,thresh,[],0,S,[]); % numerical G_u
0017 end