function err=model_laplace_test

n1=10; n2=5; h1=0.5;h2=2;
a=-1; b=5;c=3;

jmf=model_laplace(n1,n2,h1,h2);

K=assemble(jmf.Kloc,jmf.ien,jmf.id);

disp('test of model_laplace with linear solution, result should be zero')

x=jmf.xyz(:,1);
y=jmf.xyz(:,2);

u=a*x+b*y+c;
f=K*u;
load=zeros(n1+1,n2+1);
load(2:n1,1)=-a*h2;
load(2:n1,1+n2)=a*h2;
load(1,2:n2)=-b*h1;
load(1+n1,2:n2)=b*h1;
load(1,1)=0.5*(-a*h2-b*h1);
load(1+n1,1)=0.5*(-a*h2+b*h1);
load(1+n1,1+n2)=0.5*(a*h2+b*h1);
load(1,1+n2)=0.5*(a*h2-b*h1);
err=norm(f-load(:),1);