clear *
makeplots=0;
R=1; % radius of circle
[nodes,elem,edges]=ReadMeshTriangle("./circle.1");
[la,vec]=FEMEig(nodes,elem,edges,1,0,1,2,1e-6);

gset nokey
figure(1);
ShowSolution(nodes,elem,vec(:,2))


nn=21;
v=linspace(-R,R,nn);

[xx,yy]=meshgrid(v,v);
xy=[xx(:),yy(:)];

[values, grad] =FEMValue(xy,nodes,elem,vec(:,2),0);
figure(2);
ShowVectorField(xy,grad)
gset xrange [-1.2:1.2]
gset yrange [-1.2:1.2]
replot

grad=FEMGradient(nodes,elem,vec(:,2));
g1=FEMValue(xy,nodes,elem,grad(:,1),0);
g2=FEMValue(xy,nodes,elem,grad(:,2),0);
figure(3);
ShowVectorField(xy,[g1';g2']')
gset xrange [-1.2:1.2]
gset yrange [-1.2:1.2]
replot

if makeplots print("gradient2.eps","-depsc")endif


