0001 function [c,a,f,b]=vgpf(p,u,lam)
0002
0003 mu1=p.mu1; mu2=p.mu2; c=[1;0;0;1;1;0;0;1;1;0;0;1;1;0;0;1];
0004 x=p.points(1,:)'; y=p.points(2,:)'; r=x.^2+y.^2;
0005 pot=p.pa*pdeintrp(p.points,p.tria,r);
0006 ui=pdeintrp(p.points,p.tria,u);
0007 u1=ui(1,:); v1=ui(2,:);u2=ui(3,:); v2=ui(4,:);
0008 a1=u1.^2+v1.^2; a2=u2.^2+v2.^2;
0009 f1=(a1+0.5*a2).*u1;f2=(a1+0.5*a2).*v1;
0010 f3=(a2+0.5*a1).*u2;f4=(a2+0.5*a1).*v2;
0011 f=[f1;f2;f3;f4];
0012 xi=pdeintrp(p.points,p.tria,x);yi=pdeintrp(p.points,p.tria,y);
0013 a=[-mu1+pot;-mu1+pot;-mu2+pot;-mu2+pot];
0014 b=zeros(p.neq*p.neq*2,p.nt);
0015 b(3,:)=yi;b(4,:)=-xi;b(9,:)=-yi;b(10,:)=xi;
0016 b(23,:)=yi;b(24,:)=-xi;b(29,:)=-yi;b(30,:)=xi;
0017 b=lam*b;