%% CONVERGENCE: 3 CASES ON THE SAME PLOT
% Each .mat file contains a table with columns:
% Mesh, Deflection, Stress
%
% Output:
% 1) Deflection vs Mesh with all 3 cases
% 2) Stress vs Mesh with all 3 cases

clear; clc; close all;

%% INPUT FILES
matFiles = { ...
    'TAB_tweel_convergence.mat', ...
    'TAB_tweel_conv_caso2.mat', ...
    'TAB_tweel_conv_caso3.mat'};

caseNames = { ...
    'Case 1 - Best overall', ...
    'Case 2 - Most critical', ...
    'Case 3 - Lightest feasible'};

nCases = numel(matFiles);

allMesh   = cell(nCases,1);
allDefl   = cell(nCases,1);
allStress = cell(nCases,1);

for ic = 1:nCases
    %% 1) Load the .mat file
    S = load(matFiles{ic});
    vars = fieldnames(S);

    T = [];
    for k = 1:numel(vars)
        if istable(S.(vars{k}))
            T = S.(vars{k});
            break;
        end
    end

    if isempty(T)
        error('Nel file %s non ho trovato nessuna tabella.', matFiles{ic});
    end

    fprintf('\nFile letto: %s\n', matFiles{ic});
    disp(T.Properties.VariableNames')

    %% 2) Find the columns
    meshCol   = pickVarName(T, {'Mesh'});
    deflCol   = pickVarName(T, {'Deflection'});
    stressCol = pickVarName(T, {'Maximum stress'});

    %% 3) Extract the data
    meshVals   = T.(meshCol);
    deflVals   = asNumeric(T.(deflCol));
    stressVals = asNumeric(T.(stressCol));

    %% 4) Sort by mesh
    if isnumeric(meshVals)
        [meshSorted, idx] = sort(double(meshVals), 'ascend');
        xVals = meshSorted;
        xTickLabels = string(meshSorted);
    else
        meshStr = string(meshVals);
        idx = sortMeshLabels(meshStr);
        xVals = 1:numel(idx);
        xTickLabels = meshStr(idx);
    end

    allMesh{ic}   = xVals;
    allDefl{ic}   = deflVals(idx);
    allStress{ic} = stressVals(idx);

    % also save the mesh labels of the first case
    if ic == 1
        commonX = xVals;
        commonLabels = xTickLabels;
    end
end

%% ===== PLOT 1: DEFLECTION =====
figure('Color','w');
hold on; grid on; box on;

for ic = 1:nCases
    plot(commonX, allDefl{ic}, '-o', ...
        'LineWidth', 3, ...
        'MarkerSize', 8, ...
        'DisplayName', caseNames{ic});
end

xlabel('Mesh',Interpreter='latex',FontSize=20);
ylabel('Vertical deflection [m]',Interpreter='latex',FontSize=20);
% title('Mesh convergence of the three selected cases - Deflection');
ax = gca;
ax.FontSize = 25;      % larger axis numbers
ax.LineWidth = 1;    % slightly thicker axes
ax.TickLabelInterpreter = "latex";

legend('Location','northwest',Interpreter='latex',FontSize=20);

%% ===== PLOT 2: STRESS =====
figure('Color','w');
hold on; grid on; box on;

for ic = 1:nCases
    plot(commonX, allStress{ic}, '-o', ...
        'LineWidth', 3, ...
        'MarkerSize', 8, ...
        'DisplayName', caseNames{ic});
end

xlabel('Mesh', Interpreter='latex',FontSize=20);
ylabel('Maximum stress [Pa]',Interpreter='latex',FontSize=20);
% title('Mesh convergence of the three selected cases - Stress');

xticks(commonX);
xticklabels(commonLabels);
ax = gca;
ax.FontSize = 25;      % larger axis numbers
ax.LineWidth = 1;    % slightly thicker axes
ax.TickLabelInterpreter = "latex";

legend('Location','northwest',Interpreter='latex',FontSize=20);

%% ===== LOCAL FUNCTIONS =====

function name = pickVarName(T, candidates)
    namesOriginal = string(T.Properties.VariableNames);
    namesNorm = normalizeNames(namesOriginal);
    candNorm = normalizeNames(string(candidates));

    name = "";
    for i = 1:numel(candNorm)
        idx = find(contains(namesNorm, candNorm(i)), 1, 'first');
        if ~isempty(idx)
            name = namesOriginal(idx);
            return;
        end
    end

    error('Non trovo una colonna compatibile con: %s', strjoin(cellstr(string(candidates)), ', '));
end

function out = normalizeNames(s)
    out = lower(regexprep(string(s), '[^a-zA-Z0-9]', ''));
end

function x = asNumeric(x)
    if isnumeric(x)
        x = double(x);
    elseif iscell(x)
        x = cellfun(@convertOne, x);
    elseif isstring(x)
        x = str2double(x);
    elseif iscategorical(x)
        x = str2double(string(x));
    else
        x = double(x);
    end
end

function y = convertOne(v)
    if isnumeric(v)
        y = double(v);
    elseif isstring(v) || ischar(v)
        y = str2double(string(v));
    else
        y = double(v);
    end
end

function idx = sortMeshLabels(labels)
    labels = lower(strtrim(string(labels)));
    score = nan(size(labels));

    for i = 1:numel(labels)
        s = labels(i);
        if contains(s, "coarse")
            score(i) = 1;
        elseif contains(s, "medium")
            score(i) = 2;
        elseif contains(s, "fine") && ~contains(s, "very")
            score(i) = 3;
        elseif contains(s, "very fine") || contains(s, "veryfine")
            score(i) = 4;
        else
            score(i) = 100 + i;
        end
    end

    [~, idx] = sort(score, 'ascend');
end