% 
clear; clc;
format short e

[Gauss_point_reference_1D, Gauss_coefficient_reference_1D] = generate_Gauss_reference_1D(3);

%% Domain [0,1]
left = 0; right = 1;

%% generate mesh
N_array = 2.^(2:7);
error = zeros(size(N_array));
order = zeros(size(error));
tic;
for kk = 1:length(N_array)
    N = N_array(kk);
    h = (right - left) / N;
    P = left + h * (0:N);
    T = [1:N
        2:N+1];

    %% linear FE mesh
    Pb = P;
    Tb = T;
    boundarynodes = [-1 -1
        1 N+1];

    Nb = N+1; % number of the basis function
    Nlb = 2; % number of the local basis function for linear FE

    %% Assemble the stiffness matrix A
    A = sparse(Nb,Nb);
    for n = 1:N
        vertices = P(T(:,n));
        [Gauss_point_local_1D, Gauss_coefficient_local_1D] = generate_Gauss_local_1D(vertices(1),vertices(2),Gauss_point_reference_1D, Gauss_coefficient_reference_1D);
        local_S = zeros(Nlb,Nlb);
        for alpha = 1:Nlb
            for beta = 1:Nlb
                %             r = integral(@(x) coe(x) .* FE_basis_function(x, vertices, alpha, 1) .* FE_basis_function(x, vertices, beta, 1),vertices(1),vertices(2));
                %             A(Tb(beta, n), Tb(alpha, n)) = A(Tb(beta, n), Tb(alpha, n)) + r;
%                 local_S(beta, alpha) = integral(@(x) coe(x) .* FE_basis_function(x, vertices, alpha, 1) .* FE_basis_function(x, vertices, beta, 1),vertices(1),vertices(2));
                local_S(beta, alpha) = Gauss_int_vol_trial_test('coe', Gauss_point_local_1D, Gauss_coefficient_local_1D, vertices, alpha, 1,beta, 1);
            end
        end
        A(Tb(:, n), Tb(:, n)) = A(Tb(:, n), Tb(:, n)) + local_S;
    end

    %% Assemble the load vector b
    b = zeros(Nb,1);
    for n = 1:N
        vertices = P(T(:,n));
        [Gauss_point_local_1D, Gauss_coefficient_local_1D] = generate_Gauss_local_1D(vertices(1),vertices(2),Gauss_point_reference_1D, Gauss_coefficient_reference_1D);
        local_d = zeros(Nlb,1);
        for beta = 1:Nlb
            %         r = integral(@(x) f(x) .* FE_basis_function(x, vertices, beta, 0),vertices(1),vertices(2));
            %         b(Tb(beta, n), 1) = b(Tb(beta, n), 1) + r;
%             local_d(beta,1) = integral(@(x) f(x) .* FE_basis_function(x, vertices, beta, 0),vertices(1),vertices(2));
            local_d(beta,1) = Gauss_int_vol_test('f', Gauss_point_local_1D, Gauss_coefficient_local_1D, vertices, beta, 0);
        end
        b(Tb(:, n), 1) = b(Tb(:, n), 1) + local_d;
    end

    %% Deal with the Dirichlet boundary condition
    nbn = size(boundarynodes,2);
    for k = 1:nbn
        if -1 == boundarynodes(1, k)
            i = boundarynodes(2, k);
            A(i,:) = 0;
            A(i,i) = 1;
            b(i) = g(Pb(i));
        end
    end

    %% Solve linear system
    u = A \ b;

    error(kk) = max(abs(u - exact(Pb.')));
end
toc;

%% compute order of convergence and print the numerical results:
order(2:end) = log(error(1:end-1) ./ error(2:end)) ./ log(N_array(2:end)./N_array(1:end-1));
fprintf('N_array  error      order\n');
for kk = 1:length(N_array)
    fprintf('%4d    %.4e  %.4f\n',N_array(kk), error(kk), order(kk));
end

%% results:
% Elapsed time is 0.444839 seconds.
% N_array  error      order
%    4    2.3340e-03  0.0000
%    8    5.8317e-04  2.0008
%   16    1.4645e-04  1.9936
%   32    3.6675e-05  1.9975
%   64    9.1700e-06  1.9998
%  128    2.2929e-06  1.9997

function r = g(x)
r = exact(x);
end

function r = exact(x)
r = x .* cos(x);
end

function r = coe(x)
r = exp(x);
end

function r = f(x)
r = -exp(x).*(cos(x) - 2*sin(x) - x .* cos(x) - x .*sin(x));
end

function r = Gauss_int_vol_test(f, Gauss_point_local_1D, Gauss_coefficient_local_1D, vertices, test_index, test_der_order)
r = 0;
Gpn = length(Gauss_point_local_1D);
for k = 1:Gpn
    r = r + Gauss_coefficient_local_1D(k) * feval(f,Gauss_point_local_1D(k)) .* FE_basis_function(Gauss_point_local_1D(k), vertices, test_index, test_der_order);
end
end

function r = Gauss_int_vol_trial_test(coe, Gauss_point_local_1D, Gauss_coefficient_local_1D, vertices, trial_index, trial_der_order, test_index, test_der_order)
r = 0;
Gpn = length(Gauss_point_local_1D);
for k = 1:Gpn
    r = r + Gauss_coefficient_local_1D(k) * feval(coe,Gauss_point_local_1D(k)) .* FE_basis_function(Gauss_point_local_1D(k), vertices, trial_index, trial_der_order) .* FE_basis_function(Gauss_point_local_1D(k), vertices, test_index, test_der_order);
end
end

function r = FE_basis_function(x, vertices, basis_index, der_order)
h = vertices(2) - vertices(1);
if 1 == basis_index
    if 0 == der_order
        r = (vertices(2) - x) / h;
    elseif 1 == der_order
        r = -1 / h;
    elseif 2 == der_order
        r = 0;
    else
        warning = 'derivate order is wrong'
    end
elseif 2 == basis_index
    if 0 == der_order
        r = (x - vertices(1)) / h;
    elseif 1 == der_order
        r = 1 / h;
    elseif 2 == der_order
        r = 0;
    else
        warning = 'derivate order is wrong'
    end
else
    warning = 'basis_index is wrong'
end
end


function [Gauss_point_local_1D, Gauss_coefficient_local_1D] = generate_Gauss_local_1D(a,b,Gauss_point_reference_1D, Gauss_coefficient_reference_1D)
Gauss_point_local_1D = (b-a)/2 * Gauss_point_reference_1D + (b+a)/2;
Gauss_coefficient_local_1D = (b-a)/2 * Gauss_coefficient_reference_1D;
end

function [Gauss_point_reference_1D, Gauss_coefficient_reference_1D] = generate_Gauss_reference_1D(n)
syms x;
Poly = diff((x^2-1)^(n+1),x,n+1) / factorial(n+1) / 2^(n+1);
Gauss_point_reference_1D = solve(Poly);
dPoly = diff(Poly, x, 1);
Gauss_coefficient_reference_1D = zeros(size(Gauss_point_reference_1D));
for k = 1:length(Gauss_coefficient_reference_1D)
    Gauss_coefficient_reference_1D(k) = 2 / ((1-Gauss_point_reference_1D(k).^2) .* subs(dPoly,x,Gauss_point_reference_1D(k)).^2);
end
end