% Solve the Poisson equation for a rectangular 
% domain with the following Dirichlet boundary values:
% Top and bottom boundaries have temperature 100C,
% Left and right boundaries have temperature 20C.
%
% Samuli Siltanen Dec 2015

% Choose the size of the problem
row = 3;
col = 4;

% Construct the system matrix
disp('Constructing system matrix')
A = FD_Laplace(row,col);
disp('System matrix ready')

% Construct the Dirichlet data
vec_t = 100*ones(col,1);
vec_b = 100*ones(col,1);
vec_l = 20*ones(row,1);
vec_r = 20*ones(row,1);

% Initialize the right-hand side vector
m = zeros(row*col,1);

% Construct the right-hand side vector
% by looping over elements
disp('Constructing rhs')
for iii = 1:row
    for jjj = 1:col
        ind = (jjj-1)*row+iii;
        if jjj==1
            m(ind) = m(ind)+vec_l(iii);
        end
        if jjj==col
            m(ind) = m(ind)+vec_r(iii);
        end
        if iii==1
            m(ind) = m(ind)+vec_t(jjj);
        end
        if iii==row
            m(ind) = m(ind)+vec_b(jjj);
        end
    end
end
disp('Rhs constructed')


% Solve the Poisson equation
disp('Solving')
x = A\m;
x = reshape(x,row,col);

% Take a look at the system matrix
%figure(1)
%clf
%spy(A)

% Take a look at the system matrix
figure(2)
clf
imagesc(x)
axis equal
axis off
colormap gray
