% 2D P1 triangular elements
% Fernando Guevara Vasquez 2010

f = @(x,y) 2*(y.*(1-y) + x.*(1-x)) ;
utrue = @(x,y) x.*(1-x).*y.*(1-y);

% build grid of 0,1
k=1; errs=[]; hs=[];
for Nx=[10,20,30,40,50,60],
%for Nx=10,
Ny=Nx;
x=linspace(0,1,Nx+1); y=linspace(0,1,Ny+1);
[X,Y]=meshgrid(x,y);

px=X(:); py=Y(:);
tri = delaunay(px,py);
% find the nodes that are on the boundary
ib=find( abs(px-0)<1e-10 | abs(px-1)<1e-10 ...
      | abs(py-0)<1e-10 | abs(py-1)<1e-10 );

% assemble stiffness matrix and right hand side
A = sparse(length(px),length(px));
F = zeros(length(px),1);
G = [-1 1 0; -1 0 1];
Mloc = (1/24)*[2 1 1; 1 2 1; 1 1 2];
for e=1:size(tri,1),
 xloc = px(tri(e,:)); yloc = py(tri(e,:));
 B = [xloc(2:3)'-xloc(1); yloc(2:3)'-yloc(1)];
 A(tri(e,:),tri(e,:)) = A(tri(e,:),tri(e,:)) + 0.5* det(B)*G'*((B'*B)\G);
 F(tri(e,:)) = F(tri(e,:)) + det(B)*Mloc*f(xloc,yloc);
end;


for i=1:length(ib), A(ib(i),:) = 0; A(ib(i),ib(i))=1; end;
F(ib) = 0;

% solve problem
U=A\F;

Utrue=utrue(px,py);

% plot solution
trimesh(tri,px,py,U)

L2err=0;
for e=1:size(tri,1),
 errloc = U(tri(e,:)) - Utrue(tri(e,:));
 L2err = L2err + det(B)*errloc'*Mloc*errloc;
end;

H1err=0;
for e=1:size(tri,1),
 errloc = U(tri(e,:)) - Utrue(tri(e,:));
 H1err = H1err + det(B)*errloc'*Mloc*errloc + 0.5* det(B)*errloc'*G'*((B'*B)\G)*errloc;
end;

L2errs(k)=sqrt(L2err);
H1errs(k)=sqrt(H1err);
hs(k)=1/Nx;
k=k+1;
pause(0.1)
end;
