function poi2D( )
%
% -div(grad u) = f, in [0,1]x[0,1], 
%            u = 0, on boundary             
%
clear all, clc
% triangulation
g = [2 0 1 0 0 1 0; 
     2 1 1 0 1 1 0; 
     2 1 0 1 1 1 0; 
     2 0 0 1 0 1 0]';
[p,e,t] = initmesh(g,'hmax',0.05);
figure(1); clf
pdemesh(p,e,t)
% assemble
[A,b] = assemble(p,e,t,'f');
% solve 
U = A\b;
% visualize
figure(2); clf
pdesurf(p,t,U)
xlabel('x'), ylabel('y'), zlabel('U(x,y)')
% subroutines -----------------------------------------------------------------

function z = f(x,y)
z = 2*pi^2*sin(pi*x).*sin(pi*y);

function [A,b] = assemble(p,e,t,f)
Nt = size(t,2);
Np = size(p,2);
Ne = size(e,2);
A = sparse(Np,Np);
b = zeros(Np,1);
for i = 1:Nt
  n = t(1:3,i);
  x = p(1,n); 
  y = p(2,n); 
  dx = [y(2)-y(3); y(3)-y(1); y(1)-y(2)]; 
  dy = [x(3)-x(2); x(1)-x(3); x(2)-x(1)];
  area = 0.5*abs(x(2)*y(3)-y(2)*x(3)-x(1)*y(3)+y(1)*x(3)+x(1)*y(2)-y(1)*x(2));
  A(n,n) = A(n,n) + (dx*dx'+dy*dy')/4/area; 
  b(n) = b(n) + area/12*[2 1 1; 1 2 1; 1 1 2]*feval('f',x,y)';
end
% BC
for i = 1:Ne
  n = e(1,i);
  A(n,n) = 1e6;
  b(n) = 0;
end