output_precision=3;
# define all the functions

clear dG

function res = dG(xy)
[n,m]=size(xy);
res=ones(n,1);
for k=1:n
  if(xy(n,1)>1) res(n)=0;endif
endfor
endfunction

[nodes,elem,edges]=ReadMesh("Fluss");
[A,b,n2d]=FEMEquation(nodes,elem,edges,1,0,0,'dG',0);
u=FEMSolveSym(nodes,A,b,n2d,'dG');

figure(1);
gset zrange [*:*]
gset nokey
ShowSolution(nodes,elem,u)


xv=linspace(0,5,41);
yv=linspace(0,2,21);

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

[values, grad] =FEMValue(xy,nodes,elem,u,0);

figure(2);
ShowVectorField(xy,-grad)

