function convergence_2dFEM()

a = 1; hs = 0.5*2.^[0:-1:-6];

load GaussQuadratureForTriangles.mat
order = 7;
w = weights{order};
ksi = points{order};
NGP = length(w);

uanalyt = @(x1,x2) (x1.*(1-x1).*x2.*(1-x2));
duanalytdx1 = @(x1,x2) ((1-2*x1).*x2.*(1-x2));
duanalytdx2 = @(x1,x2) (x1.*(1-x1).*(1-2*x2));

for i=1:length(hs)
    h = hs(i);
    [P,T,ID,IN] = discretizeRectangle2(a,a,h);
    N = size(P,2)
    m = size(T,2);
    nD = length(ID); uhat = zeros(nD,1);
    mN = size(IN,2); tauhat = zeros(mN,1);

    massP = (1/3) * (P(:,T(1,:))+P(:,T(2,:))+P(:,T(3,:)));
    k = ones(1,m);
    f = 2 * ( massP(1,:).*(1-massP(1,:)) + massP(2,:).*(1-massP(2,:)) );
    [A,b,utilde] = FEM(P,T,{'aloc_2dLaplace'},{k},{'bloc_2dsource'},{f},ID,uhat,IN,tauhat);

    uH = A\b;
    u = utilde;
    nonDirNodes = setdiff(1:N,ID);
    u(nonDirNodes) = u(nonDirNodes) + uH;

%     close all
%     trisurf(T',P(1,:),P(2,:),u')
%     hold on; trisurf(T',P(1,:),P(2,:),uanalyt(P(1,:),P(2,:)))

    l2err(i) = 0;
    h1err(i) = 0;
    for e=1:m
        Te = T(:,e);
        Re = [P(:,Te(2))-P(:,Te(1)),P(:,Te(3))-P(:,Te(1))]; P1e = P(:,Te(1));
        xe = P1e*ones(1,NGP) + Re*ksi;
        ue = u(Te);
        integrande = ( uanalyt(xe(1,:),xe(2,:)) - (ue(1)*(1-ksi(1,:)-ksi(2,:))+ue(2)*ksi(1,:)+ue(3)*ksi(2,:)) ).^2;
        l2err(i) = l2err(i) + abs(det(Re)) * w * integrande';
        gradue = inv(Re)'*[-1 1 0;-1 0 1]*ue;
        integrande = ( duanalytdx1(xe(1,:),xe(2,:)) - gradue(1) ).^2 + ( duanalytdx2(xe(1,:),xe(2,:)) - gradue(2) ).^2;
        h1err(i) = h1err(i) + abs(det(Re)) * w * integrande';
    end
    l2err(i) = sqrt(l2err(i))
    h1err(i) = sqrt(h1err(i))
end

plot(log10(hs),log10(l2err),'r',log10(hs),log10(hs.^2),'r--',log10(hs),log10(h1err),'b',log10(hs),log10(hs.^1),'b--')

