% use finite elements to solve -Lap u= fn on arb domain with 0 BC
%    discrete eqn  B u = F; loop by triangles
function fe3
h=.1;
fd=@(p) p(:,1).^2/4+p(:,2).^2-1; % define ellipse
tic; [p,t]=distmesh2d(fd,@huniform,h,[-2,-1;2,1],[]); toc
tic, n=size(p,1), nt=size(t,1), geps=1e-3;
ind=findind(p,geps,fd); % reorder nodes so bdry pts at end
p=p(ind,:); t(:,1)=ind(t(:,1)); t(:,2)=ind(t(:,2)); t(:,3)=ind(t(:,3));
for i=n:-1:1  if abs(fd(p(i,:)))>geps nint=i, break; end, end % nint=# of interior nodes
B=sparse(nint,nint); F=zeros(nint,1); interior=1:nint;
fn=inline('-exp(p(:,1)+p(:,2)).*(p(:,1)+4*p(:,2)+.5+p(:,1).^2/2+2*p(:,2).^2)','p'); % RHS
ue=(p(interior,1).^2/4+p(interior,2).^2-1).*exp(p(interior,1)+p(interior,2)); % exact soln
area=abs(((p(t(:,3),1)-p(t(:,1),1)).*(p(t(:,2),2)-p(t(:,1),2))...
    -(p(t(:,2),1)-p(t(:,1),1)).*(p(t(:,3),2)-p(t(:,1),2))))/2;
centroid=(p(t(:,1),:)+p(t(:,2),:)+p(t(:,3),:))/3;
for j=1:nt
    nodes=sort(t(j,:)); % node # of triangle j
    M=[p(nodes,1), p(nodes,2) ones(3,1)];    
    abc=M\eye(3); % coef of basis fn for node i on triangle j,   phi_i=abc(:,i)'*[x;y;1]    
    abc=abc(1:2,:); maxnode=3; 
    if nodes(3)>nint maxnode=2; end, if nodes(2)>nint maxnode=1; end, if nodes(1)>nint maxnode=0; end
    for i=1:maxnode % calculate contribution of triangle j to diag entries
        ni=nodes(i); B(ni,ni)=B(ni,ni)+sum(abc(:,i).^2)*area(j);
    end    
    % calculate contribution to off-diag entries
    for i=2:maxnode B(nodes(1),nodes(i))=B(nodes(1),nodes(i))+abc(:,1)'*abc(:,i)*area(j); end
    if maxnode==3 B(nodes(2),nodes(3))=B(nodes(2),nodes(3))+abc(:,2)'*abc(:,3)*area(j); end
    c=area(j)/3*fn(centroid(j,:)); % calculate contribution to Load vtr
    F(nodes(1:maxnode))=F(nodes(1:maxnode))+c;
end
B=B+triu(B,1)';
u=B\F; toc
figure(1), trisurf(t,p(:,1),p(:,2),[u; zeros(n-nint,1)]); 
[sqrt((u-ue)'*B*(u-ue)), norm(u-ue)*h],

function ind=findind(p,geps,fd) 
% permute index of nodes so that bdry nodes occur at end
n=size(p,1); i=1; in=n; ind=[1:n]';
while i<=in
    if abs(fd(p(i,:)))<geps 
        ind(i)=i; found=0;
        for j=in:-1:i+1
            if abs(fd(p(j,:)))>geps ind(i)=j; ind(j)=i; in=j-1; found=1; break; end
        end
        if found==0 return; end
    end
    i=i+1;
end