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

CreateMeshTriangle("Fluss2",...
  [0 0 2;5 0 1;5 2 2;3 2 2;3 0.5 2;2 0.5 2;1 2 2; 0 2 1],0.002);
[nodes,elem,edges]=ReadMeshTriangle("Fluss2.1");

[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)
gset xrange [0:5.1]
gset yrange [0:2]
% print("FlussVoll.eps","-depsc")

xv=linspace(1.7,2.3,31);
yv=linspace(0.3,0.6,13);

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

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

figure(3);
ShowVectorField(xy,-grad)
gset xrange [1.7:2.35]
% print("FlussDetail.eps","-depsc")

npoints=31;
xv=linspace(0,5,npoints);
yv=0.25*ones(1,npoints);
[values,grad]=FEMValue([xv',yv'],nodes,elem,u,0);
figure(4);
gset yrange [0:0.5]
gset xrange [0:5]
plot(xv,-grad(:,1))
% print("FlussSpeedProfile.eps","-depsc")
gset autoscale

