clear all
close all
precon = input('Preconditioner? (0 = No, 1 = Yes) ');
reord = input('Reordering of the points? (0 = No, 1 = Yes) ');
m1 = 40;
[x,y,bcidx] = meshsquare(m1);
m = length(x);
ell = delaunay(x,y);
l = size(ell,1);
triplot(ell,x,y)
hold on
j = 1;
plot(x(ell(j,:)),y(ell(j,:)),'r*')
gfun = @(x,y) 2*(1-x).*x+2*(1-y).*y; % NON e` lineare, no esattezza
gell = zeros(l,1);
delta = zeros(l,1);
dummy = spalloc(m,m,9*m);
for j = 1:l
  for k = 1:3
    dummy(ell(j,k),ell(j,k)) = 1;
    for h = k+1:3
      dummy(ell(j,k),ell(j,h)) = 1;
      dummy(ell(j,h),ell(j,k)) = 1;
    end
  end
end
figure
spy(dummy)
title('sparsity pattern before reordening')
A = spalloc(m,m,nnz(dummy));
P = spalloc(m,m,nnz(dummy));
if (reord)
% riordinamento
  p = symrcm(dummy);
  q = 1:m;
  q(p) = q;
  for j = 1:l
    for k = 1:3
      ell(j,k) = q(ell(j,k));
    end
  end
  x = x(p);
  y = y(p);
  for j = 1:length(bcidx)
    bcidx(j) = q(bcidx(j));
  end
end
a = zeros(l,3);
b = zeros(l,3);
c = zeros(l,3);
g = zeros(m,1);
for j = 1:l
  gell(j) = (gfun(x(ell(j,1)),y(ell(j,1)))+...
	     gfun(x(ell(j,2)),y(ell(j,2)))+...
	     gfun(x(ell(j,3)),y(ell(j,3))))/3;
  delta(j) = det([1,1,1;...
		  x(ell(j,1)),x(ell(j,2)),x(ell(j,3));...
		  y(ell(j,1)),y(ell(j,2)),y(ell(j,3))]) / 2;
  a(j,1) = x(ell(j,2))*y(ell(j,3))-x(ell(j,3))*y(ell(j,2));
  b(j,1) = -y(ell(j,3))+y(ell(j,2));
  c(j,1) = x(ell(j,3))-x(ell(j,2));		 
  a(j,2) = -x(ell(j,1))*y(ell(j,3))+x(ell(j,3))*y(ell(j,1));
  b(j,2) = y(ell(j,3))-y(ell(j,1));
  c(j,2) = x(ell(j,1))-x(ell(j,3));		 
  a(j,3) = x(ell(j,1))*y(ell(j,2))-x(ell(j,2))*y(ell(j,1));
  b(j,3) = -y(ell(j,2))+y(ell(j,1));
  c(j,3) = x(ell(j,2))-x(ell(j,1));		 
end
for j = 1:l
  for k = 1:3
    A(ell(j,k),ell(j,k)) = A(ell(j,k),ell(j,k))+...
	b(j,k)*b(j,k)/(4 * abs(delta(j)))+...
	c(j,k)*c(j,k)/(4 * abs(delta(j)));
    P(ell(j,k),ell(j,k)) = P(ell(j,k),ell(j,k))+...
	abs(delta(j))/6;
    for h = k+1:3
      A(ell(j,k),ell(j,h)) = A(ell(j,k),ell(j,h))+...
	  b(j,k)*b(j,h)/(4 * abs(delta(j)))+...
	  c(j,k)*c(j,h)/(4 * abs(delta(j)));
      P(ell(j,k),ell(j,h)) = P(ell(j,k),ell(j,h))+abs(delta(j))/12;
      A(ell(j,h),ell(j,k)) = A(ell(j,k),ell(j,h));
      P(ell(j,h),ell(j,k)) = P(ell(j,k),ell(j,h));
    end
    g(ell(j,k)) = g(ell(j,k))+gell(j)*abs(delta(j))/3;
  end
end
figure
spy(A)
title('sparsity pattern after reordening')
%Dirichlet b.c.
for i = bcidx
  A(i,i) = 1e16;
  g(i) = 0;
end
if (precon)
  R = cholinc(A,'0');
  disp('norm of R''R-A')
  normest(R'*R-A)
  u = pcg(A,g,0.1/m1^2,100,R',R);
else
  u = pcg(A,g,0.1/m1^2,100);
end
figure
uexact = x.*(1-x).*y.*(1-y);
plot3(x,y,u,'*')
norm(u-uexact,inf)