clear fF gNF one zero
h=logspace(-0.5,-2.0,11);
error=zeros(size(h));
xy=[0,0,1; 1,0,1; 1,1,2;0,1,1];

function res=exactF(xy)
  res=sin(xy(:,1)).*sin(xy(:,2));
endfunction

function res = fF(xy)
  res=-2*sin(xy(:,1)).*sin(xy(:,2));
endfunction


function res = gNF(xy)
  res=sin(xy(:,1)).*cos(xy(:,2));
endfunction

function res = one(xy)
  res=ones(size(xy(:,1)));
endfunction

function res = zero(xy)
  res=zeros(size(xy(:,1)));
endfunction


for k=1:length(h);
area=h(k)**2;

CreateMeshTriangle("test",xy,area)
[nodes,elem,edges]=ReadMeshTriangle("./test.1");


[A,b,n2d]=FEMEquation(nodes,elem,edges,'one','zero','fF','exactF','gNF');
u=FEMSolveSym(nodes,A,b,n2d,'exactF');
# ShowSolution(nodes,elem,u);

error(k)=max(abs(u-exactF(nodes)));
endfor


gset title
gset nokey
gset xlabel "log h"
gset ylabel "log error"
plot(log(h),log(error))



