clear;

// Finite Elements Method in 2D-dimension space for elliptic PDE's - Delta u(x) + lambda u(x) = f(x). Delta: Laplacian operator
// Domain: 2D unit-disk


// Parameters

h = 0.1 // Step size for numerical integration
bound = "Dirichlet"  // Boundary conditions (Dirichlet or Neumann)
lambda = 0

disp("===================================")
disp("2D Finite Element Method")
disp(["Boundary condition:", bound])
disp("===================================")

// Functions and vectors

function z=Test(x)
    z = 1-x(1)-x(2) + (1e-12)
endfunction

function z=f(x)
    // Second member function
    // Inputs:
    // - x: Space variable, size: (2,1)
    z = 4 + (1e-12)*x(1)
endfunction





function z=grad(g,x)
    // Gradient of a function with evaluation at point x
    // Inputs:
    // - g: Function with 2D input variable
    // - x: Space variable (evaluation point), size: (2,1)
    z(1) = (g([x(1)+h,x(2)]) - g([x(1)-h,x(2)]))/(2*h)
    z(2) = (g([x(1),x(2)+h]) - g([x(1),x(2)-h]))/(2*h)
endfunction

function I = Int_Delta(g, X1, X2, X3)
    // Integration of function g on triangle with verticies X1, X2 and X3. Numerical approximation with existing numerical integration
    // Inputs:
    // - g: Function with 2D input variable
    // - X1, X2, X3: Vertices of the triangle. Size: (1,2)

    function zz = gg(x1,x2)
        // Adapted function for numerical integration
        x = [x1, x2]
        zz = g(x)
    endfunction
    
    I = int2d([X1(1);X2(1);X3(1)], [X1(2);X2(2);X3(2)], gg)
endfunction

function aaa = aa(u, v, X1, X2, X3)
    // Bilinear form in the weak formulation of the PDE on the triangle with vertices X1, X2 and X3
    // Inputs:
    // - u, v: Functions with 2D inputs and scalar outputs
    // - X1, X2, X3: Vertices of the triangle. Shape: (2,1)
    function z = g(x)
        // Function which is integrated on the triangle
        //z = (u([x(1)+h,x(2)]) - u([x(1)-h,x(2)]))/(2*h)
        z = sum(grad(u,x).*grad(v,x)) + lambda*u(x)*v(x)
    endfunction
    aaa = Int_Delta(g, X1, X2, X3)
endfunction

function LLL = LL(v, X1, X2, X3)
    // Linear form in the weak formulation of the PDE on the triangle with vertices X1, X2 and X3
    // Inputs:
    // - v: Function with 2D input and scalar output
    // - X1, X2, X3: Vertices of the triangle. Shape: (2,1)
    function z = g(x)
        // Function which is integrated on the triangle
        z = f(x)*v(x)
    endfunction
    LLL = Int_Delta(g, X1, X2, X3)
endfunction








function [C,V] = FE(J,K)
    // Gives mesh (cells and vertices) for finite elements
    // Inputs:
    // - J: Number of verticies at the edge of the domain
    // - K: Number of layers of points inside the unit disk
    V = zeros(J*K+1, 3)
    for k=0:K-1
        J0 = (1-(-1)^k)/4 + int(k/2)
        V(1+k*J:1+(k+1)*J-1,1) =  (1-k/K)*cos(2*%pi*(J0:1:J+J0-1)/J)'
        V(1+k*J:1+(k+1)*J-1,2) =  (1-k/K)*sin(2*%pi*(J0:1:J+J0-1)/J)'
    end

    V((J)*K+1,1) = 0
    V((J)*K+1,2) = 0
    
    V(1:J,3) = ones(J,1) // 1 for belong to bound
    
    
    C1 = []
    C2 = []
    C3 = []
    for k = 0:K-2
        C1 = [C1 ; [(J*k+1:J*k+J)']]
        C2 = [C2 ; [(J*k+2:J*k+J)';J*k+1]]
        C3 = [C3 ; [(J*k+J+1:J*k+2*J)']]

        C1 = [C1 ; [(J*k+J+1:J*k+2*J)']]
        C2 = [C2 ; [(J*k+J+2:J*k+2*J)';J*k+J+1]]
        C3 = [C3 ; [(J*k+2:J*k+J)';J*k+1]]
    end
    
    C1 = [C1 ; [(J*K-J+1:J*K)']]
    C2 = [C2 ; [(J*K-J+2:J*K)';J*K-J+1]]
    C3 = [C3 ; [(J*K+1)*ones(J,1)]]
    
    C = [C1,C2,C3]

endfunction

function Plot_Domain(J,K,save_fig)
    // Plots mesh domain (cells and verticies)
    // Inputs:
    // - J: Number of verticies at the edge of the domain
    // - K: Number of layers of points inside the unit disk
    // - save_fig: Saves the figure or not. Boolean.
    
    [C, V] = FE(J,K)
    scf()
    scatter(V(:,1),V(:,2),50,"red",".")
    for j = 1:size(V)(1)
        xstring(V(j,1),V(j,2),string(j))
    end
    isoview("on")
    
    //disp(C)
    
    for k = 1:size(C)(1)
        j1 = C(k,1)
        j2 = C(k,2)
        j3 = C(k,3)
        plot([V(j1,1),V(j2,1),V(j3,1),V(j1,1)],[V(j1,2),V(j2,2),V(j3,2),V(j1,2)],"color",13,"linewidth",2)
    end
    
    if save_fig == %T then
        filename="C:\Users\Utilisateur\Documents\Scripts_Scilab\Finite_Elements_2D\Figures_Mesh_Domain\Finite_Element_Mesh_Domain_J="+string(J)+"_K="+string(K)+"_bound="+bound;
        xs2pdf(0,filename);
        xs2pdf(gcf(),filename);
    end
    
endfunction


function Plot_Local_Basis(save_fig)
    // Plots local basis associated to a triangle cell of the mesh
    // Inputs:
    // - save_fig: Saves the figure or not. Boolean
    
    AA = 45
    TT = 15
    
    X = [[-1,0];[2,0];[0,2]]
    M = [ones(3,1),X]
    a = (M\eye(3,3))'
    
    Z = a(:,1) + a(:,2).*X(:,1) + a(:,3).*X(:,2)
    
    scf()
    
    subplot(1,3,1)
    vertex = [X,[Z(1);0;0]]
    face = [1,2,3]
    xvf = matrix(vertex(face,1),size(face,1),length(vertex(face,1))/size(face,1))'
    yvf = matrix(vertex(face,2),size(face,1),length(vertex(face,1))/size(face,1))'
    zvf = matrix(vertex(face,3),size(face,1),length(vertex(face,1))/size(face,1))'
    tcolor = [2]
    plot3d(xvf,yvf,list(zvf,tcolor),alpha=AA,theta=TT)
    gce().hiddencolor = -1
    isoview("on")
    title("Function 1")
    gca().grid=[1 1 1]
    
    subplot(1,3,2)
    vertex = [X,[0;Z(2);0]]
    face = [1,2,3]
    xvf = matrix(vertex(face,1),size(face,1),length(vertex(face,1))/size(face,1))'
    yvf = matrix(vertex(face,2),size(face,1),length(vertex(face,1))/size(face,1))'
    zvf = matrix(vertex(face,3),size(face,1),length(vertex(face,1))/size(face,1))'
    tcolor = [3]
    plot3d(xvf,yvf,list(zvf,tcolor),alpha=AA,theta=TT)
    gce().hiddencolor = -1
    isoview("on")
     title("Function 2")
    gca().grid=[1 1 1]
    
    subplot(1,3,3)
    vertex = [X,[0;0;Z(3)]]
    face = [1,2,3]
    xvf = matrix(vertex(face,1),size(face,1),length(vertex(face,1))/size(face,1))'
    yvf = matrix(vertex(face,2),size(face,1),length(vertex(face,1))/size(face,1))'
    zvf = matrix(vertex(face,3),size(face,1),length(vertex(face,1))/size(face,1))'
    tcolor = [5]
    plot3d(xvf,yvf,list(zvf,tcolor),alpha=AA,theta=TT)
    gce().hiddencolor = -1
    isoview("on")
    title("Function 3")
    gca().grid=[1 1 1]
    
    h = get('current_figure')
    h.figure_size=[900,500]
    if save_fig == %T then
        filename='C:\Users\Utilisateur\Documents\Scripts_Scilab\Finite_Elements_2D\Figures\Local_Basis';
        xs2pdf(0,filename);
        xs2pdf(gcf(),filename);
    end
endfunction
       
function FEM(J, K, save_fig)
    // Finite Elements Method to solve elliptic PDE on the unit disk in two-dimensional space, with homogeneous Dirichlet or Neumann boundary conditions
    // Inputs:
    // - J: Number of verticies at the edge of the domain
    // - K: Number of layers of points inside the unit disk
    // - save_fig: Saves the figure or not. Boolean
    
    if bound == "Dirichlet" then
        A = zeros(J*K+1-J,J*K+1-J)
        B = zeros(J*K+1-J,1)
    elseif bound == "Neumann" then
        A = zeros(J*K+1,J*K+1)
        B = zeros(J*K+1,1)
    end
    
    
    disp('Mesh creation.')
    [C, V] = FE(J,K)
    
    disp('Matrices constructions.')
    for k = 1:size(C)(1)
        disp("Loading [%]")
        disp(int(1000*k/size(C)(1))/10)
        X1 = V(C(k,1),1:2)
        X2 = V(C(k,2),1:2)
        X3 = V(C(k,3),1:2)
        M = [[1,X1];[1,X2];[1,X3]]
        
        for betaa = 1:3
            i = C(k, betaa)
            e_beta = eye(3,3)(:,betaa)
            a_beta = M\e_beta
            
            function z = phi_beta(x)
                // Function of local basis with index alpha
                // Inputs
                // - x: Space variable. Shape: (2,1)
                z = a_beta(1) + a_beta(2)*x(1) + a_beta(3)*x(2)
            endfunction
            
            if bound == "Dirichlet" then
                if (V(i,3) == 0) then
                    B(i-J) = B(i-J) + LL(phi_beta, X1, X2, X3)
                end
            elseif bound == "Neumann" then
                B(i) = B(i) + LL(phi_beta, X1, X2, X3)
            end   
                
            for alphaa = 1:3 
                j = C(k, alphaa)
                e_alpha = eye(3,3)(:,alphaa)
                a_alpha = M\e_alpha      
                 
                function z = phi_alpha(x)
                    // Function of local basis with index alpha
                    // Inputs
                    // - x: Space variable. Shape: (2,1)
                    z = a_alpha(1) + a_alpha(2)*x(1) + a_alpha(3)*x(2)
                endfunction
                               
                if bound == "Dirichlet" then
                    if (V(i,3) == 0) && (V(j,3) == 0) then
                        A(i-J,j-J) = A(i-J,j-J) + aa(phi_alpha, phi_beta, X1, X2, X3)
                    end
                elseif bound == "Neumann" then
                    A(i,j) = A(i,j) + aa(phi_alpha, phi_beta, X1, X2, X3)
                end             
            end
        end
    end
    
    //disp(A,B)
    
    disp('Resolution of linear system.')
    U = A\B
    
    //disp(U)
    
    if bound == "Dirichlet" then
        U = [zeros(J,1);U]
    end
    
    Error = abs(U - (1-V(:,1).^2-V(:,2).^2))
    
    disp('Plotting solution.')
    vertex = [V(:,1:2),U]
    vertex_error = [V(:,1:2),Error]
    face = C
    xvf = matrix(vertex(face,1),size(face,1),length(vertex(face,1))/size(face,1))'
    yvf = matrix(vertex(face,2),size(face,1),length(vertex(face,1))/size(face,1))'
    zvf = matrix(vertex(face,3),size(face,1),length(vertex(face,1))/size(face,1))'
    zvf_error = matrix(vertex_error(face,3),size(face,1),length(vertex(face,1))/size(face,1))'
        
    tcolor = []
    Color = []
    
    for k = 1:size(C)(1)
        Delta_U = max(U) - min(U)
        col = int(255*((U(C(k,1)) + U(C(k,2)) + U(C(k,3)))/(3) - min(U))/Delta_U)
        Color = [Color; col, 0, 255 - col]
        col = color(col, 0, 255 - col)
        tcolor = [tcolor; col]
    end

    scf()
    Color = Color/255
    h = get('current_figure')
    h.figure_size=[900,500]
    
    
    
    subplot(1,2,1);
    //plot3d(xvf,yvf,list(zvf,tcolor),alpha=55,theta=45)
    plot3d1(xvf,yvf,zvf,alpha=55,theta=45)
    gce().hiddencolor = -1
    isoview("on")
    title("Finite Element Method")

    subplot(1,2,2);
    //plot3d(xvf,yvf,list(zvf,tcolor),alpha=55,theta=45)
    plot3d1(xvf,yvf,zvf_error,alpha=0,theta=-90)
    gce().hiddencolor = -1
    isoview("on")    
    h.color_map = unique(Color,'r')
    colorbar(min(Error),max(Error),[0,length(unique(tcolor))])
    gca().grid=[1 1 1]
    title("Error")
    
    
    if save_fig == %T then
        filename="C:\Users\Utilisateur\Documents\Scripts_Scilab\Finite_Elements_2D\Figures\Finite_Element_Method_2D_J="+string(J)+"_K="+string(K)+"_bound="+bound;
        xs2pdf(0,filename);
        xs2pdf(gcf(),filename);
    end
endfunction



