clc;clear all;close all;

bcdialog = [' Lower x boundary condition:'
    ' Upper x boundary condition:'
    ' Lower y boundary condition:'
    ' Upper y boundary condition:'];

bc = zeros(4,3); % zero the boundary conditions
disp('       Solution of the Laplace or Poisson Equation')
disp('(Two-dimensional elliptic partial differential equation)')
disp(' '); disp(' ')
disp('                        Upper y boundary                    '             )
disp('                                                            '             )
disp('         __________________________________________________ '             )
disp('        |                                                  |'             )
disp('        |                                                  |'             )
disp('        |                                                  |'             )
disp('        |                                                  |'             )
disp('Lower x |                                                  | Upper x'     )
disp('boundary|                                                  | boundary'    )
disp('        |                                                  |'             )
disp('        |                                                  |'             )
disp('        y                                                  |'             )
disp('        |__x_________________length (L)____________________|'             )
disp('       0                                                    '             )
disp('                        Lower y boundary                    '             )
disp(' ');
L = 1;%input(' Length of the object (x-direction) (cm) = ');
W = 1;%input(' Width of the object  (y-direction) (cm) = ');
nx = 10;%input(' Number of divisions in x-direction = ');
ny = 10;%input(' Number of divisions in y-direction = ');
disp(' ');
f = -0.1;%rho/T--input('Right-hand side of the Poisson Equation (f) = ');
disp(' '); guess = 0.02;%input('Initial starting guess = ');
disp(' '); disp(' Boundary conditions:')

for k = 1:4
    disp(' ')
    disp(bcdialog(k,:))
    disp(' 1 - Dirichlet')
    disp(' 2 - Neumann')
    disp(' 3 - Robbins')
    bc(k,1) = 1;%input(' Enter your choice : ');
    
    if bc(k,1) == 1
        bc(k,2) = 0;%input(' Value of the Dirichlet boundary = ');
    end
    
    if bc(k,1) == 2
        bc(k,2) = input(' Value of Neumann boundary = ');
    end
    
    if bc(k,1) == 3
        disp(' Value of Robbins boundary:');
        disp(' u'' = (beta) + (gamma)*u')
        bc(k,2) = input('Constant (beta) = ');
        bc(k,3) = input('Coefficient (gamma) = ');
    end
end

dx = L/nx; dy = W/ny;
x = [0:nx]*dx; y = [0:ny]*dy;

u = guess*ones(nx+1,ny+1);
U = u;
count = 0; iter = 0;
total_points = (nx+1)*(ny+1);
while count < (total_points)
    % Set the boundary conditions
    % Lower x
    i=1;
    switch bc(1,1)
        case 1
            u(i,:) = bc(1,2);
        case {2,3}
            u(i,:) = (4*u(i+1,:)-u(i+2,:)-2*dx*bc(1,2))/(3+2*dx*bc(1,3));
    end
    % Upper x
    i=nx+1;
    switch bc(2,1)
        case 1
            u(i,:) = bc(2,2);
        case {2,3}
            u(i,:) = (4*u(i-1,:)-u(i-2,:)+2*dx*bc(2,2))/(3-2*dx*bc(2,3));
    end
    % Lower y
    j=1;
    switch bc(3,1)
        case 1
            u(:,j) = bc(3,2);
        case {2,3}
            u(:,j) = (4*u(:,j+1)-u(:,j+2)-2*dy*bc(3,2))/(3+2*dy*bc(3,3));
    end
    % Upper y
    j=ny+1;
    switch bc(4,1)
        case 1
            u(:,j) = bc(4,2);
        case {2,3}
            u(:,j) = (4*u(:,j-1)-u(:,j-2)+2*dy*bc(4,2))/(3-2*dy*bc(4,3));
    end
    for i = 2:nx
        for j = 2:ny
            u(i,j) = ((u(i+1,j)+u(i-1,j))/dx^2+(u(i,j+1)...
                +u(i,j-1))/dy^2-f)/(2*(1/dx^2+1/dy^2));
        end
    end
    
    % Examine convergence
    for i = 1:nx+1
        for j = 1:ny+1
            if u(i,j)~=0
                if abs((U(i,j)-u(i,j))/u(i,j)) <= 1e-12
                    count = count+1;
                end
            else
                if abs((U(i,j)-u(i,j))) <= 1e-12
                    count = count+1;
                end
            end
        end
    end
    U = u;
end

fprintf('\n Iterations = %g ',iter)
fprintf('\n Points converged = %g/%g \n\n', count, total_points)
disp(' The matrix of results is:');
u'
%Plot the final results
figure;surf(x,y,u');zlim([-max(u(:)),max(u(:))])
xlabel('Length (cm)');ylabel('Width (cm)');
zlabel('Displacement (cm)')
shading interp
title('Displacement of the Stretched Membrane');pause(1);
pause;
%% Hyperpolic PDE to evaluate what happens when force is removed from the membrane in a vaccum (it will oscillate for infinity!)
un_1 = reshape(u,[numel(u),1]);
un = un_1;
dt = 0.05;
nt = 200;
t = transpose(0:nt-1)*dt;
U_hyper = zeros(numel(un),nt);
U_hyper(:,1) = un;
U_hyper(:,2) = un;
asq = 2;%tension of membrane divided by density of membrane (higher values will vibrate faster)
A = zeros(numel(un));
for n = 1:numel(un)
    j_n = floor(n/size(u,1));
    i_n = n-size(u,1)*(j_n-1);
    for m = 1:numel(un)
        j_m = floor(m/size(u,1));
        i_m = m-size(u,1)*(j_m-1);
%         if i_n == 1 || i_n == numel(un) || j_n == 1 || j_n == numel(un)
%             A(n,m) = 1;
        if i_n == i_m && j_n == j_m
            A(n,m) = 2*asq/dx^2+1/dt^2;
        elseif (i_n == i_m && j_n == j_m+1) || (i_n == i_m && j_n == j_m-1) || (i_n == i_m+1 && j_n == j_m) || (i_n == i_m-1 && j_n == j_m)
            A(n,m) = -asq/(2*dx^2);
        end
    end
end

B = zeros(numel(un));
for n = 1:numel(un)
    j_n = floor(n/size(u,1));
    i_n = n-size(u,1)*(j_n-1);
    for m = 1:numel(un)
        j_m = floor(m/size(u,1));
        i_m = m-size(u,1)*(j_m-1);
%         if i_n == 1 || i_n == numel(un) || j_n == 1 || j_n == numel(un)
%             B(n,m) = 1;
        if i_n == i_m && j_n == j_m
            B(n,m) = -2*asq/dx^2-1/dt^2;
        elseif (i_n == i_m && j_n == j_m+1) || (i_n == i_m && j_n == j_m-1) || (i_n == i_m+1 && j_n == j_m) || (i_n == i_m-1 && j_n == j_m)
            B(n,m) = asq/(2*dx^2);
        end
    end
end
figure;
for n = 3:nt
    rhs = B*U_hyper(:,n-2)+2/dt^2*U_hyper(:,n-1);
    temp = reshape(A\rhs,size(u));
    for i = 1:size(u,1)
        for j = 1:size(u,2)
            if i == 1 || i == size(u,1) || j == 1 || j == size(u,2)
                temp(i,j) = 0;
            end
        end
    end
    U_hyper(:,n) = reshape(temp,[numel(temp),1]);
    frame = transpose(reshape(U_hyper(:,n),size(u)));
    surf(x,y,frame);zlim([-1.5*max(u(:)),1.5*max(u(:))])
    xlabel('Length (cm)');ylabel('Width (cm)');
    zlabel('Displacement (cm)')
    shading interp
    title('Displacement of the Stretched Membrane when Released from pressure (no friction)');pause(.1);
end
% figure;
% for n = 1:size(U_hyper,2)
%     frame = transpose(reshape(U_hyper(:,n),size(u)));
%     surf(x,y,frame);zlim([-max(u(:)),max(u(:))])
%     xlabel('Length (cm)');ylabel('Width (cm)');
%     zlabel('Displacement (cm)')
%     shading interp
%     title('Displacement of the Stretched Membrane');pause(.5);
% end