clc; clear
resultFolder = 'result_M_convergence\';
D = load(['result\','FEMdata.mat'], 'data');
[X_fem, Y_fem] = meshgrid(linspace(0,1,1201));
U1_fem = D.data.U1;
U2_fem = D.data.U2;
W3_fem = D.data.W3;
deltas = [0.05, 0.1, 0.15, 0.2];
mvalues = [2, 3, 4, 5, 6,7];

ErrU = zeros(numel(mvalues), numel(deltas));
ErrU2 = zeros(numel(mvalues), numel(deltas));
ErrW3 = zeros(numel(mvalues), numel(deltas));
dx = zeros(numel(mvalues), numel(deltas));
gOrder = 2;

for j = 1:numel(deltas)
    delta = deltas(j);
    % ndivs = ceil(1/(delta./mvalues));
    fileName = [resultFolder,sprintf('Model_II_M_Conv_Order%d_delta%.2f.mat', gOrder, delta)];
    
    if(~exist(fileName, 'file'))
        continue;
    end
    A = load(fileName);
    modelData = A.Model_II_M_Conv;
    
    result = modelData(numel(modelData)).result;
    X_fem = result.X;
    Y_fem = result.Y;
    VOL_fem = result.VOL;
    U1_fem = result.U1;
    U2_fem = result.U2;
    W3_fem =  result.W3;
    
    for i = numel(modelData)-1:-1:1
        if(isempty(modelData(i).result))
            continue;
        end
        result = modelData(i).result;
        X = result.X;
        Y = result.Y;
        VOL = result.VOL;
        U1 = result.U1;
        U2 = result.U2;
        W3 =  result.W3;
        mvalue = modelData(i).mvalue;
        ndiv =  modelData(i).ndiv;
        gfunOrder = modelData(i).gfunOrder;
       
        U1_f = interp2(X_fem,Y_fem,U1_fem,X,Y);
        U2_f = interp2(X_fem,Y_fem,U2_fem,X,Y);
        ErrU(i, j) = sqrt(sum(sum((U1-U1_f).^2+(U2-U2_f).^2))/sum(sum(U1_f.^2+U2_f.^2)));
         %ErrU(i, j) = sqrt(sum(sum((U1.^2+U2.^2))))/numel(U1);
%         figure
%         pcolor(X,Y,U2_f)
%         colormap(jet(20)); colorbar
%         figure
%         pcolor(X,Y,U2)
%         colormap(jet(20)); colorbar
        
%         W3_f = W3_fem(1:1200/ndiv:end,1:1200/ndiv:end);
%         ErrW3(i,j) = sqrt(sum(sum((W3-W3_f).^2))/sum(sum(W3_f.^2)));
%         dx(j,i) = 1 / ndiv;
    end
end

markers = {'o', 'h', 'd', '>', 'p'};
lines = {'-', '--', '-.', '-', ':'};
colors = {'r', 'g', 'b', 'k', 'm'};

for j = 1:numel(deltas)
    figure(1);
    plot(mvalues, (ErrU(:,j)), [colors{j}, markers{j}, lines{j}], 'markerfacecolor', colors{j}, 'markersize', 6)
    hold on
    xlabel('$\it{m} \rm{value}$', 'interpreter', 'latex', 'fontname', 'times new roman',...
        'fontsize', 16);
    ylabel('$L_2(u)$', 'interpreter', 'latex', 'fontname', 'times new roman',...
        'fontsize', 16)
    set(gca, 'fontname', 'times new roman', 'fontsize', 16, 'xminortick', 'on', 'yminortick', 'on')
    
    % p1 = polyfit(log10(dx(:,j)), log10(ErrU(:,j)), 1)

    figure(2);
    plot(log10(dx(:,j)), log10(ErrW3(:,j)), [colors{j}, markers{j}, lines{j}], 'markerfacecolor', colors{j}, 'markersize', 6)
    hold on
    xlabel('$\rm{log_{10}}{\rm{\Delta}\it{x}}$', 'interpreter', 'latex', 'fontname', 'times new roman',...
        'fontsize', 16);
    ylabel('$\rm{log_{10}}{||\it{\varepsilon(\it{\theta_{\rm{3}}})}||_{\it{L}_{\rm{2}}}}$', 'interpreter', 'latex', 'fontname', 'times new roman',...
        'fontsize', 16)
    set(gca, 'fontname', 'times new roman', 'fontsize', 16, 'xminortick', 'on', 'yminortick', 'on')
p1 = polyfit(log10(dx(:,j)), log10(ErrW3(:,j)), 1);

end
if(gOrder==1)
    lgds = {
        '$m=1$', '$m=2$', '$m=3$', '$m=4$', '$m=5$'
        };
else
    lgds = {
        '$m=2$', '$m=3$', '$m=4$', '$m=5$'
        };
end
if(gOrder==2)
    figure(1);
    axis equal
    xmin = -2.4;   xmax = -0.8; xinc = 0.4; xarr = xmin:xinc:xmax;
    ymin = -2.7;   ymax = -0.7; yinc = 0.5; yarr = ymin:yinc:ymax;
    axis([xmin, xmax, ymin, ymax])
    legend(lgds, 'interpreter', 'latex', 'Location', 'northwest');
    set(gca, 'xtick', xarr, 'ytick', yarr)
    figure(2);
    axis equal
    xmin = -2.5;   xmax = -0.5; xinc = 0.5; xarr = xmin:xinc:xmax;
    ymin = -2.3;   ymax = 0.2;  yinc = 0.5; yarr = ymin:yinc:ymax;
    axis([xmin, xmax, ymin, ymax])
    legend(lgds, 'interpreter', 'latex', 'Location', 'northwest');
    set(gca, 'xtick', xarr, 'ytick', yarr)
end