function err=model_laplace_test(n1,n2,h1,h2)

n1=10; n2=5; h=0.5;h2=2;
a=-1; b=5;
jmf=model_laplace(n1,n2,h1,h2);

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

disp('test with known exact solution u=a*x+b*y, result should be zero')

x=jmf
u=a*x+b*y;
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);