%codigo para resolucao do exercicio 3 da lista 6
%
% construcao do grafo de malha retangular uniforme
%  ("n" pontos horizontais "m" pontos verticais)
%
%   REDE
%
%   X-----X-----X-----X-----X-----X
%   |     |     |     |     |     |
%   X-----X-----X-----X-----X-----X
%   |     |     |     |     |     |
%   X-----X-----X-----X-----X-----X
%   |     |     |     |     |     |
%   n+1---X-----X-----X-----X-----X
%   |     |     |     |     |     |
%   1-----2-----3-----4-----5-----n
%      I    II    III
%
%   INDEXACAO da matriz conec:
%     As (n-1)*m primeiras linhas sao as conectividades horizontais.
%     Na linha (n-1)*m+1 iniciam-se as conectividades verticais,
%     dadas sequencialmente por:
%     (1,n+1), (2,n+2), ... , (n,n+n), ... , (n+1,n+n+1), ... , (n-1)*m + (m-1)*n.
n=50; m=50; 
nv = n*m;
nc = (n-1)*m + (m-1)*n;
% construimos a conectividade aresta-no do grafo
% conexoes horizontais
for j=1:m
  for i=1:n-1
  k = (j-1)*(n-1) + i;
  con(k,1) = (j-1)*n + i;
  con(k,2) = (j-1)*n + i+1;
  end
end
% conexoes verticais
for i=1:n
  for j=1:m-1
  k = (n-1)*m  +  (i-1)*(m-1) + j;
  con(k,1) = (i-1)*(m-1) + j;
  con(k,2) = (i-1)*(m-1) + j+n;
  end
end
% construcao da matriz do circuito
a = sparse(nv,nv); % mais rapido que zeros(nv,nv), matriz cheia
for ic=1:nc
    k1=con(ic,1);
    k2=con(ic,2);
    aux=1;
    a(k1,k1)=a(k1,k1)+aux;
    a(k2,k2)=a(k2,k2)+aux;
    a(k1,k2)=a(k1,k2)-aux;
    a(k2,k1)=a(k2,k1)-aux;
end
% imposicao de condicoes externas
atil = a;
b = zeros(nv,1);
% nos com probabilidade imposta
%tampa superior: metade ezquerda=0, metade direita=1
for k=1:n
n1=n*(m-1)+k;
atil(n1,:)=0; atil(n1,n1)=1; 
b(n1,1)=0; if (k>n/2) b(n1,1)=1;end
end
% resolucao do sistema
x = atil\b;
% reshape para grafico
x2 = reshape(x,n,m);
contourf(x2')
pause
% podemos ate mudar a parte da tampa que interessa
for kf=2:n-1
 for k=1:n
 n1=n*(m-1)+k;
 b(n1,1)=0; if (k>kf) b(n1,1)=1;end
 end
 x=atil\b;
 x2 = reshape(x,n,m);
 contourf(x2')
 pause
end