%% Approximating a function using piecewise elements
%
% So we want to approximate the function
% u\left(x\right)=2x\;\sin \left(2\pi x\right)+3\;\mathrm{for}\;x\in \left\lbrack 0,1\right\rbrack
% using piecewise elements of order 1 or 2.
% This means finding all nodal values u_i to the approximation (the blue circles in the figure above). This problem can be formulated using Galerkin as such: Find u_i such that
% \int_0^1 \left(\sum_i \varphi_i \left(x\right){\;u}_i -u\left(x\right)\right)\varphi_j \;d\;x\;=0
% with U\left(x\right)=\sum_i \varphi_i \left(x\right){\;u}_i
% expanding the expression in the integral we get
% \int_0^1 \left(\sum_i \varphi_i \left(x\right)\varphi_j \left(x\right){\;u}_i -u\left(x\right)\varphi_j \left(x\right)\right)\;d\;x\;=0
% \int_0^1 \sum_i \varphi_i \left(x\right)\varphi_j \left(x\right){\;u}_i \;d\;x-\int_0^1 u\left(x\right)\varphi_j \left(x\right)\;d\;x=0
% move u_i \;out
% \sum_i \int_0^1 \varphi_i \left(x\right)\varphi_j \left(x\right)\;d\;x\;\;{\;u}_i =\int_0^1 u\left(x\right)\varphi_j \left(x\right)\;d\;x
% which is equivalent to the matrix system
% M_{\textrm{ij}} u_i =f_i or \mathit{\mathbf{M}}\;\mathit{\mathbf{u}}=\mathit{\mathbf{f}}
% Solving this system gets us all nodal values for u_i.
% Utilizing a piecewise approximation leads to the element mass matrix and element load vector:
% M_e =\int_{x_1 }^{x_2 } {\underline{\varphi} }^T \underline{\varphi} \;d\;x
% and
% f_e =\int_{x_1 }^{x_2 } {\underline{\varphi} }^T u\left(x\right)\;d\;x
% which need to be assembled into the global mass matrix and load vectors.
%
% We will do this in code now using linear element and symbolic integration.

clear

syms x
ue(x) = 2*x*sin(2*pi*x)+3;

N = 4;
nnod = N+1;
xnod = linspace(0,1,N+1);
nodes = [(1:nnod-1)', (2:nnod)'];

Phi = @(x,x1,x2) [(x2-x)/(x2-x1), (x-x1)/(x2-x1)];

M = zeros(nnod,nnod); F = zeros(nnod,1);
for iel = 1:N
    inod = nodes(iel,:);
    xc = xnod(inod);
    x1 = xc(1); x2 = xc(2);
    phi = @(x)Phi(x,x1, x2);

    Me = int(phi(x).'*phi(x), x,x1,x2);
    fe = int(phi(x).'*ue(x), x,x1,x2);
    
    M(inod,inod) = M(inod,inod) + Me;
    F(inod) = F(inod) + fe;
end

u = M\F;

figure
fplot(ue,[0,1],'r','LineWidth',2); hold on
plot(xnod,u,'-bo')

% No we will compute the average error between the exact function and our Galerkin approximation, we shall use the L_2 norm as a measure of error. This basicaly means squaring the difference between the functions over the entire domain and taking the square root.
% L_2 \textrm{error}:=\sqrt{\int_0^1 {\left(U-u\right)}^2 \;d\;x}
% So we need to compute the difference between the function for each element and adding it all up.
% L_2 \textrm{error}:=\sqrt{\int_0^1 {\left(U-u\right)}^2 \;d\;x}=\sqrt{\sum_e^N \int_e^{\;} {\left(U_e -u\right)}^2 \;d\;x}
% where
% U_e \left(x\right)=\underline{\varphi} \left(x\right){\mathit{\mathbf{u}}}_e \overset{\textrm{lin}\;\textrm{ele}}{=} \left\lbrack \begin{array}{cc}
% \varphi_1  & \varphi_1
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% u_1 \\
% u_2
% \end{array}\right\rbrack

L2err = 0;
for iel = 1:N
    inod = nodes(iel,:);
    x1 = xnod(inod(1)); x2 = xnod(inod(2));
    phi = @(x)Phi(x,x1, x2);
    U = phi(x)*u(inod);
    L2err = L2err + double(int((U-ue)^2, x,x1,x2));
end
L2err = sqrt(L2err)

%
%% Error convergence study
% Now we can start changing N to see how the error goes down. If we plot N versus the error we get a convergence plot.

clear
syms x
ue(x) = 2*x*sin(2*pi*x)+3;

NN = [2,4,8,16,32];
L2Err = [];
Time = [];

figure
fplot(ue,[0,1],'r','LineWidth',2); hold on
hp = plot(NaN,NaN,'-b.');

for N = 4
    tic
    nnod = N+1;
    xnod = linspace(0,1,N+1);
    nodes = [(1:nnod-1)', (2:nnod)'];
    
    Phi = @(x,x1,x2) [(x2-x)/(x2-x1), (x-x1)/(x2-x1)];
    
    M = zeros(nnod,nnod); F = zeros(nnod,1);
    for iel = 1:N
        inod = nodes(iel,:);
        xc = xnod(inod);
        x1 = xc(1); x2 = xc(2);
        phi = @(x)Phi(x,x1, x2);
    
        Me = vpaintegral(phi(x).'*phi(x), x,x1,x2);
        fe = vpaintegral(phi(x).'*ue(x), x,x1,x2);
        
        M(inod,inod) = M(inod,inod) + Me;
        F(inod) = F(inod) + fe;
    end
    
    u = M\F;
    t1 = toc;
    Time(end+1) = t1;

    L2err = 0;
    for iel = 1:N
        inod = nodes(iel,:);
        x1 = xnod(inod(1)); x2 = xnod(inod(2));
        phi = @(x)Phi(x,x1, x2);
        U = phi(x)*u(inod);
        L2err = L2err + double(vpaintegral((U-ue)^2, x,x1,x2));
    end
    L2err = sqrt(L2err);

    L2Err(end+1) = L2err;

    set(hp,'XData',xnod, 'YData',u);
    title(sprintf("Approximation using %d linear elements",N))
    pause(0.1)
    drawnow
end

%
% Quadratic approximation

syms x1 x2 x3 x
A = [1 x1 x1^2
     1 x2 x2^2
     1 x3 x3^2];
C = A\eye(3);
Phi = simplify([1 x x^2]*C);
Phi = matlabFunction(Phi);

Time2 = []; L2Err2 = [];

figure;
for N = 4
    tic
    nnod = N*2+1;
    xnod = linspace(0,1,nnod);
    nodes = [(1:2:nnod-2)', (2:2:nnod-1)', (3:2:nnod)'];
    neq = nnod*1;
    
    M = zeros(nnod,nnod); F = zeros(nnod,1);
    for iel = 1:N
        inod = nodes(iel,:);
        xc = xnod(inod);
        x1 = xc(1); x2 = xc(2); x3 = xc(3);
        phi = @(x)Phi(x,x1, x2, x3);
    
        Me = vpaintegral(phi(x).'*phi(x), x,x1,x2);
        fe = vpaintegral(phi(x).'*ue(x), x,x1,x2);
    
        M(inod,inod) = M(inod,inod) + Me;
        F(inod) = F(inod) + fe;
    end
    
    u = M\F;
    t1 = toc;
    Time2(end+1) = t1;
    
    L2err = 0;
    for iel = 1:N
        inod = nodes(iel,:);
        xc = xnod(inod);
        x1 = xc(1); x2 = xc(2); x3 = xc(3);
        phi = @(x)Phi(x,x1, x2, x3);
        U = phi(x)*u(inod);
        L2err = L2err + double(vpaintegral((U-ue)^2, x,x1,x2));
    end
    L2err = sqrt(L2err);
    L2Err2(end+1) = L2err;
    
    
    clf
    fplot(ue,[0,1],'r','LineWidth',1); 
    hold on
    
    for iel = 1:N
        inod = nodes(iel,:);
        xc = xnod(inod);
        x1 = xc(1); x2 = xc(2); x3 = xc(3);
        phi = @(x)Phi(x,x1, x2, x3);
        U = phi(x)*u(inod);
    
        fplot(U, [x1,x3])
        plot(xc,u(inod),'.b');
    end
    title(sprintf("Approximation using %d quadratic elements",N))
%     if N == 2
%         exportgraphics(gca,"PiecewiseQuadratic.gif","Append",false)
%     end
%     for i = 1:10
%         exportgraphics(gca,"PiecewiseQuadratic.gif","Append",true)
%     end
end



%


figure; hold on
plot(NN,L2Err,'-ko','DisplayName','Linear')
plot(NN,L2Err2,'-ks','DisplayName','Quadratic')
set(gca,'xscale','log','yscale','log')
xlabel("N"); ylabel("L_2 error")
xticks(NN)
xticklabels({'2','4','8','16','32'})
grid on
legend show

figure; hold on
plot(NN,Time,'-ko','DisplayName','Linear')
plot(NN,Time2,'-ks','DisplayName','Quadratic')
set(gca,'xscale','log','yscale','log')
xlabel("N"); ylabel("Time (s)")
xticks(NN)
xticklabels({'2','4','8','16','32'})
yticks([0.5,1,2,4,8,16])
grid on
legend("show",location="northwest")