clear;
format long;
fc=1;    % f constant 
bc=0.;     % The Dirichlet boundary condition, constant
W=1;   % the width of the rectangle
L=2.0;    % the length of the rectangle
N=20;      % the number of elements in Y direction
M=40;       % the number of elements in X direction
dx = L/M;    % dx 
dy = W/N;    % dy 
D=(N+1)*(M+1);    % The dimension of the global stiffniss matrix or number of nodes
s=zeros(D,D);      % global stiffness matrix
f=zeros(D,1);       % load factor
v=zeros(N+1,M+1);    % the collocated results for plotting
xindex=0;
yindex=0;
x1=0; x2=0; x3=0; x4=0;   % (x1,x2,x3) and (x4,x2,x3) are summits of the two trangles in one rectangle element
j1=0; j2=0;
% generating the global stiffness matrix and load vectors 
for yindex = 1:N
    for xindex = 1:M
        x1 = xindex + (yindex-1)*(M+1);
        x2 = 1+ x1;
        x3 = xindex + yindex*(M+1);
        x4 = 1+ x3;
        
        j1 = dx*dy;  % the Jacobian of the first trangle
        s(x1,x1) = s(x1,x1)+(dx^2+dy^2)/(2*j1);
        s(x1,x2) = s(x1,x2)-dy^2/(2*j1);
        s(x1,x3) = s(x1,x3)-dx^2/(2*j1);
        s(x2,x1) = s(x2,x1)-dy^2/(2*j1);
        s(x2,x2) = s(x2,x2)+dy^2/(2*j1);
        s(x2,x3) = s(x2,x3);
        s(x3,x1) = s(x3,x1)-dx^2/(2*j1);
        s(x3,x2) = s(x3,x2);
        s(x3,x3) = s(x3,x3)+dx^2/(2*j1);
        f(x1,1) = f(x1,1)+ fc*j1/6;
        f(x2,1) = f(x2,1)+ fc*j1/6;
        f(x3,1) = f(x3,1)+ fc*j1/6;   %  Process the first trangle in one block
        
        j2 = dx*dy;   % the Jacobian of the second trangle
        s(x4,x4) = s(x4,x4)+(dx^2+dy^2)/(2*j2);
        s(x4,x2) = s(x4,x2)-dx^2/(2*j2);
        s(x4,x3) = s(x4,x3)-dy^2/(2*j2);
        s(x2,x4) = s(x2,x4)-dx^2/(2*j2);
        s(x2,x2) = s(x2,x2)+dx^2/(2*j2);
        s(x2,x3) = s(x2,x3);
        s(x3,x4) = s(x3,x4)-dy^2/(2*j2);
        s(x3,x2) = s(x3,x2);
        s(x3,x3) = s(x3,x3)+dy^2/(2*j2);
        f(x4,1) = f(x4,1)+ fc*j2/6;
        f(x2,1) = f(x2,1)+ fc*j2/6;
        f(x3,1) = f(x3,1)+ fc*j2/6;  %  Process the second trangle 
    end
end

% applying BC 
for xindex=1:M+1     
    s(xindex,:) = zeros(1,D); s(xindex,xindex)=1; f(xindex,1)=bc;
    s(xindex+N*(M+1),:) = zeros(1,D);
    s(xindex+N*(M+1),xindex+N*(M+1))=1;
    f(xindex+N*(M+1),1)=bc;      
end         % the bottome and upper BC 
for yindex=1:N-1
    s(1+yindex*(M+1),:)=zeros(1,D);
    s(1+yindex*(M+1),1+yindex*(M+1))=1;
    f(1+yindex*(M+1),1)=bc;
    s((1+yindex)*(M+1),:)=zeros(1,D);
    s((1+yindex)*(M+1),(1+yindex)*(M+1))=1;
    f((1+yindex)*(M+1),1)=bc;
end         % the left and right side BC

q=s\f;      % solve q

% generating visulized results
for yindex=1:N+1
    for xindex=1:M+1
        v(yindex,xindex) = q((M+1)*(yindex-1)+xindex);
    end
end
[X,Y] = meshgrid(0:dx:L,0:dy:W);
[C,h] = contourf(X,Y,v);
clabel(C,h)
colormap cool
%axis equal tight