% Simple implementation of the Poisson equation on a square domain
clear;

domain_x = [0,1];
domain_y = [0,1];

n_x = 101;
n_y = 101;

coords_x = linspace( domain_x(1), domain_x(2), n_x);
coords_y = linspace( domain_y(1), domain_y(2), n_y);

n_vert = n_x*n_y;
vertices = zeros(n_vert,2);
k = 1;
for j = 1:n_y
    for i = 1:n_x
        coords = [coords_x(i), coords_y(j)];
        vertices(k, :) = coords;
        k = k+1;
    end
end

n_el = 2*(n_x-1)*(n_y-1);
elements = zeros(n_el, 3);
k = 1;
for j = 1:n_y-1
    for i = 1:n_x-1
        triangle_upper = [ j*n_x + i, (j-1)*n_x + i, j*n_x + i+1 ];
        triangle_lower = [ (j-1)*n_x + i+1, j*n_x + i+1, (j-1)*n_x + i ];
        elements(k, :) = triangle_upper;
        k = k+1;
        elements(k, :) = triangle_lower;
        k = k+1;
    end
end

is_boundary = zeros(n_vert,1);
bound_val = zeros(n_vert,1);
for j=1:n_x
    
    i = j;
    x = vertices(i,1);
    y = vertices(i,2);
    
    is_boundary(i) = 1;
    bound_val(i) = sin(2*pi*x);
    %bound_val(i) = 0;
end
for j=2:n_y-1
    
    i = j*n_x;
    x = vertices(i,1);
    y = vertices(i,2);
    
    is_boundary(i) = 1;
    bound_val(i) = 0;
    
    i = (j-1)*n_x+1;
    x = vertices(i,1);
    y = vertices(i,2);
    
    is_boundary(i) = 1;
    bound_val(i) = 0;
end
for j=1:n_x
    
    i = n_x*(n_y-1) + j;
    x = vertices(i,1);
    y = vertices(i,2);
    
    is_boundary(i) = 1;
    bound_val(i) = 0;
end

dof = zeros(n_vert);
k = 0;
for i = 1:n_vert
    if is_boundary(i) ~= 1
        k = k+1;
        dof(i) = k;
    end
end

n_dof = nnz(dof);
A = sparse(n_dof, n_dof);
f_Dir = zeros(n_dof,1);
f = zeros(n_dof, 1);

gradN1 = [-1; -1]; gradN2 = [1; 0]; gradN3 = [0; 1];
ref_gradients = [gradN1, gradN2, gradN3];

eps = 0.1;

for k = 1:n_el
    x1 = vertices( elements(k,1), 1); y1 = vertices( elements(k,1), 2);
    x2 = vertices( elements(k,2), 1); y2 = vertices( elements(k,2), 2);
    x3 = vertices( elements(k,3), 1); y3 = vertices( elements(k,3), 2);
    B = [x2-x1, x3-x1; y2-y1 , y3-y1];
    size = 0.5*det(B);
    center_x = 1/3*(x1+x2+x3); center_y = 1/3*(y1+y2+y3);
    
    for j = 1:3
        if is_boundary( elements(k, j) ) == 1, continue; end;
        
        dj = dof(elements(k,j));
        
        for i = 1:3
            
            gradi_ref = ref_gradients(:, i); gradj_ref = ref_gradients(:, j);
            
            gradi = (B')\gradi_ref; gradj = (B')\gradj_ref;
            integrand = eps*dot(gradi, gradj);
            integral = integrand*size;
            
            if is_boundary( elements(k, i) ) == 1
                c = bound_val(elements(k, i));
                f_Dir(dj) = f_Dir(dj) - c*integral;
            else
                di = dof(elements(k,i));
                A(di,dj) = A(di,dj) + integral;
            end
        end
        
        f(dj) = f(dj) + exp( -100*( (center_x-0.5)^2 +  (center_y-0.5)^2 ) ) * size;
        
        
    end
end

b = f+f_Dir;
u = A\b;

plot_z = zeros(n_vert,1);
for i = 1:n_vert
    if is_boundary(i) == 1
        plot_z(i) = bound_val(i);
    else
        plot_z(i) = u(dof(i));
    end
end
trisurf(elements, vertices(:,1),vertices(:,2), plot_z);
axis([0,1, 0,1, min(plot_z), max(plot_z)]);
xlabel x; ylabel y;



