Baja SAE Chassis Tube Optimization and 3D Frame Analysis

%% baja_full_A_B_C.m
% Combined analytic sweep + Pareto selection + 3D frame FEA overlay + color-mapped deformation
% A+B+C: approximate nodes from PDF + color mapping + overlay all Pareto candidates
clear; clc; close all; warning('off','all');

%% ---------------- User parameters ----------------
in2m = 0.0254; ft2m = 0.3048;
total_tube_length = 25; % m (approx)
E = 205e9; rho = 7850; g = 9.81;
min_Kb = 2620; min_Sb = 374; safety_factor = 1.5;
L_effective = 1.5; % for buckling checks (m)
maxWallRatio = 0.15; % t/D limit in original script
maxParetoToPlot = 6; % max pareto candidates per material to overlay (keeps plot readable)

%% ---------------- Materials (original values; edit cost_per_ft if you want) ----------------
materials = struct( ...
    'name', {'AISI1018', 'AISI4130', 'DocolR8'}, ...
    'Sy', {420e6, 560e6, 830e6}, ...
    'Su', {440e6, 670e6, 900e6}, ...
    'cost_per_ft', {3.50, 8.25, 12.50}, ...
    'weldability', {5,3,4}, ...
    'availability', {5,4,2}, ...
    'color', {'r','b','g'});

%% ---------------- Analytic sweep (original style) ----------------
D_outer = linspace(1.0, 2.5, 40) * in2m; % 1" to 2.5" OD
t_wall = linspace(0.049, 0.134, 25) * in2m; % 0.049" to 0.134" wall

results = {};
for m = 1:length(materials)
    for D = D_outer
        for t = t_wall
            d = D - 2*t;
            if d <= 0 || t/D > maxWallRatio, continue; end
            I = (pi/64) * (D^4 - d^4);
            A = (pi/4) * (D^2 - d^2);
            c = D/2;
            r = sqrt(I/A);
            Kb = E * I;
            Sb = (materials(m).Sy * I) / c;
            P_critical = (pi^2 * E * I) / (L_effective^2);
            J = (pi/32) * (D^4 - d^4);
            Kt = E * J / (2*(1+0.3));
            W = rho * A * g; % N/m
            volume_total = A * total_tube_length;
            weight_total = volume_total * rho; % kg
            cost_total = (total_tube_length / ft2m) * materials(m).cost_per_ft;
            specific_stiffness = Kb / W;
            specific_strength = Sb / W;
            cost_per_stiffness = cost_total / (Kb + eps);
            SF_bending = Sb / min_Sb;
            SF_stiffness = Kb / min_Kb;
            pass_basic = (Kb >= min_Kb) && (Sb >= min_Sb);
            pass_safety = (Kb >= min_Kb) && (Sb >= min_Sb * safety_factor);
            results = [results; {materials(m).name, D/in2m, t/in2m, d/in2m, ...
                A*1e6, I*1e8, Kb, Sb, P_critical/1000, Kt, ...
                W, weight_total, cost_total, specific_stiffness, ...
                specific_strength, cost_per_stiffness, SF_bending, ...
                SF_stiffness, materials(m).weldability, ...
                materials(m).availability, pass_basic, pass_safety, ...
                materials(m).color}];
        end
    end
end

resultsTable = cell2table(results, ...
    'VariableNames', {'Material', 'OD_in', 'Wall_in', 'ID_in', ...
    'Area_mm2', 'I_cm4', 'Kb_Nm2', 'Sb_Nm', 'P_crit_kN', 'Kt_Nm2', ...
    'Weight_Npm', 'TotalWeight_kg', 'TotalCost_USD', 'SpecStiffness', ...
    'SpecStrength', 'CostPerStiff', 'SF_bend', 'SF_stiff', ...
    'Weldability', 'Availability', 'Pass', 'PassSafety', 'Color'});

valid_basic = resultsTable(resultsTable.Pass == true, :);
valid_safety = resultsTable(resultsTable.PassSafety == true, :);

fprintf('\nAnalytic sweep complete: total rows = %d, pass safety = %d\n', height(resultsTable), height(valid_safety));

%% ---------------- Decision & summary (same as original) ----------------
fprintf('\n========= SUMMARY =========\n');
uniqueMats = unique(valid_safety.Material);
summary_data = {};
for i = 1:length(uniqueMats)
    mat = uniqueMats{i};
    idx = strcmp(valid_safety.Material, mat);
    matData = valid_safety(idx,:);
    if ~isempty(matData)
        [~,bestIdx] = min(matData.TotalWeight_kg);
        opt = matData(bestIdx,:);
        summary_data = [summary_data; {mat, opt.OD_in, opt.Wall_in, opt.TotalWeight_kg, opt.TotalCost_USD, opt.SF_bend, opt.Weldability, opt.Availability}];
        fprintf('Material %s: best mass OD=%.3f" wall=%.4f" mass=%.2fkg cost=$%.2f\n', mat, opt.OD_in, opt.Wall_in, opt.TotalWeight_kg, opt.TotalCost_USD);
    end
end

%% ---------------- Plots (original working charts) ----------------
colors = containers.Map({'AISI1018','AISI4130','DocolR8'},{[1 0 0],[0 0 1],[0 0.4 0]});

figure('Position',[200 200 1100 500]);
subplot(1,2,1); hold on; grid on;
title('Cost vs Weight Trade-off (With Safety Factor)');
xlabel('Total Weight (kg)'); ylabel('Total Cost (USD)');
for m = 1:length(materials)
    mat = materials(m).name;
    idx = strcmp(valid_safety.Material, mat);
    if any(idx)
        scatter(valid_safety.TotalWeight_kg(idx), valid_safety.TotalCost_USD(idx), 60, colors(mat), 'filled', 'DisplayName', mat, 'MarkerFaceAlpha',0.6);
    end
end
legend('show','Location','northwest');

subplot(1,2,2); hold on; grid on;
title('Specific Stiffness vs Specific Strength');
xlabel('Specific Stiffness'); ylabel('Specific Strength');
for m = 1:length(materials)
    mat = materials(m).name;
    idx = strcmp(valid_safety.Material, mat);
    if any(idx)
        scatter(valid_safety.SpecStiffness(idx), valid_safety.SpecStrength(idx), 60, colors(mat), 'filled', 'DisplayName', mat, 'MarkerFaceAlpha',0.6);
    end
end
legend('show','Location','northwest');

figure('Position',[200 200 1000 500]);
subplot(1,2,1); hold on; grid on;
title('Valid Tube Geometries'); xlabel('Outer Diameter (in)'); ylabel('Wall Thickness (in)');
for m = 1:length(materials)
    mat = materials(m).name;
    idx = strcmp(valid_safety.Material, mat);
    if any(idx)
        scatter(valid_safety.OD_in(idx), valid_safety.Wall_in(idx), 30, colors(mat), 'filled', 'DisplayName', mat, 'MarkerFaceAlpha',0.6);
    end
end
legend('show','Location','northwest'); xlim([1 2.5]);

subplot(1,2,2); hold on; grid on;
title('Safety Factor Distribution'); xlabel('Bending Safety Factor'); ylabel('Number of Options');
for m = 1:length(materials)
    mat = materials(m).name;
    idx = strcmp(valid_safety.Material, mat);
    if any(idx)
        histogram(valid_safety.SF_bend(idx), 'FaceColor', colors(mat), 'FaceAlpha',0.5, 'DisplayName', mat);
    end
end
yl = ylim; plot([safety_factor safety_factor], yl, 'k--','LineWidth',2,'DisplayName','Required SF');
legend('show','Location','northeast');

%% ---------------- 3D Chassis geometry (approximated from your PDF) ----------------
% Coordinates in meters approximated from the schematic (top, side, roll hoops)
% You can replace this with exact CAD-extracted coordinates later.
nodes = [
    0.00, -0.35, 0.00; %1 front-lower-left
    0.00,  0.35, 0.00; %2 front-lower-right
    1.20, -0.35, 0.00; %3 rear-lower-left
    1.20,  0.35, 0.00; %4 rear-lower-right
    0.20, -0.35, 0.12; %5 front-seat-left
    0.20,  0.35, 0.12; %6 front-seat-right
    0.60, -0.35, 0.12; %7 mid-seat-left
    0.60,  0.35, 0.12; %8 mid-seat-right
    1.05, -0.20, 0.18; %9 roll-lower-left
    1.05,  0.20, 0.18; %10 roll-lower-right
    1.05, -0.20, 0.95; %11 roll-top-left
    1.05,  0.20, 0.95; %12 roll-top-right
    0.05, -0.20, 0.18; %13 front-lower-hoop-left
    0.05,  0.20, 0.18; %14 front-lower-hoop-right
    0.05, -0.20, 0.70; %15 front-top-left
    0.05,  0.20, 0.70; %16 front-top-right
    0.55, 0.00, 0.45;  %17 center-support
    0.90, 0.00, 0.45;  %18 rear-support
    ];

elems = [
    1 2; 3 4; 1 3; 2 4;        % base rectangle
    5 6; 7 8; 5 7; 6 8;        % seat platform
    3 7; 4 8;                  % lower to seat
    7 9; 8 10; 9 11; 10 12; 11 12; % rear roll hoop
    13 14; 15 16; 13 15; 14 16; % front hoop
    1 13; 2 14;                % connect front base to hoop
    5 17; 6 17; 7 17; 8 17;    % center support
    7 18; 8 18; 9 18; 10 18;   % rear support
    17 18; 3 9; 4 10
    ];

elems = unique(sort(elems,2),'rows');

nNodes = size(nodes,1);
nElems = size(elems,1);

%% ---------------- Representative combined load for deformation viz ----------------
F_global = zeros(6*nNodes,1);
driver = 800; seat_nodes = [5 6 7 8];
for n = seat_nodes, F_global((n-1)*6 + 3) = F_global((n-1)*6 + 3) - driver/numel(seat_nodes); end
corner = 2500; % lateral roll load
F_global((11-1)*6 + 2) = -corner*0.6; F_global((12-1)*6 + 2) = -corner*0.4;

% fix base nodes (assume welded to floor or ground)
fixedNodes = [1 2 3 4];
fixedDOFs_master = false(6*nNodes,1);
for fn = fixedNodes, fixedDOFs_master((fn-1)*6 + (1:6)) = true; end

%% ---------------- For each material: compute pareto & run FEA overlays ----------------
figure('Name','3D Deformation Overlay (Pareto candidates)','Position',[100 100 1100 800]);
hold on; axis equal; grid on; view(45,25);
title('Undeformed frame (black) and Pareto candidate deformations (color by displacement)');
% plot undeformed (wireframe)
for e=1:nElems
    n1 = elems(e,1); n2 = elems(e,2);
    plot3([nodes(n1,1) nodes(n2,1)],[nodes(n1,2) nodes(n2,2)],[nodes(n1,3) nodes(n2,3)],'-k','LineWidth',1);
end

colormap(gca, 'jet');
allMaxDisp = 0;
deformationData = {}; % store per plotted overlay for legend

for mi = 1:length(materials)
    matname = materials(mi).name;
    matIdxs = find(strcmp(valid_safety.Material, matname));
    if isempty(matIdxs)
        fprintf('No passing candidates for %s — skipping overlay.\n', matname);
        continue;
    end
    sub = valid_safety(matIdxs,:);
    % produce Pareto front (minimize mass & cost)
    pts = [sub.TotalWeight_kg, sub.TotalCost_USD];
    pfMask = paretofront(pts);
    paretoSet = sub(pfMask,:);
    % If too many Pareto points, reduce to top by variety: choose up to maxParetoToPlot
    if height(paretoSet) > maxParetoToPlot
        % sort by mass*cost normalized and take best few
        normMass = (paretoSet.TotalWeight_kg - min(paretoSet.TotalWeight_kg)) / (max(paretoSet.TotalWeight_kg)-min(paretoSet.TotalWeight_kg)+eps);
        normCost = (paretoSet.TotalCost_USD - min(paretoSet.TotalCost_USD)) / (max(paretoSet.TotalCost_USD)-min(paretoSet.TotalCost_USD)+eps);
        score = normMass + normCost;
        [~, order] = sort(score, 'ascend');
        paretoSet = paretoSet(order(1:maxParetoToPlot), :);
    end
    
    % For each Pareto candidate, run FEA and collect nodal displacements
    for p = 1:height(paretoSet)
        cand = paretoSet(p,:);
        OD = cand.OD_in * in2m; t = cand.Wall_in * in2m;
        d = OD - 2*t;
        if d <= 0, continue; end
        A = pi/4 * (OD^2 - d^2);
        Ix = pi/64 * (OD^4 - d^4);
        Iy = Ix; J = pi/32*(OD^4 - d^4);
        % build section array
        section = struct('A',[],'Ix',[],'Iy',[],'J',[]);
        for e = 1:nElems
            section(e).A = A; section(e).Ix = Ix; section(e).Iy = Iy; section(e).J = J;
        end
        % assemble stiffness
        matStruct.E = E; matStruct.G = E/(2*(1+0.3)); matStruct.Sy = materials(mi).Sy;
        try
            [Kglob, ~] = assembleGlobal3D(nodes, elems, section, matStruct);
        catch ME
            warning('Assembly failed for %s OD=%.3f w=%.4f: %s', matname, cand.OD_in, cand.Wall_in, ME.message);
            continue;
        end
        % apply BCs & solve
        fixedDOFs = fixedDOFs_master;
        [Kff,Kfp,freeDOFs,Ufull] = applyBCandSolve(Kglob, F_global, fixedDOFs);
        if isempty(Ufull)
            warning('FEA solve failed for %s OD=%.3f w=%.4f', matname, cand.OD_in, cand.Wall_in);
            continue;
        end
        % nodal translational displacement
        Utrans = zeros(nNodes,3);
        for n=1:nNodes
            Utrans(n,:) = Ufull((n-1)*6 + (1:3))';
        end
        dispMag = sqrt(sum(Utrans.^2,2));
        maxDisp = max(dispMag);
        allMaxDisp = max(allMaxDisp, maxDisp); %#ok<SAGROW>
        deformationData{end+1}.mat = matname;
        deformationData{end}.candidate = cand;
        deformationData{end}.Utrans = Utrans;
        deformationData{end}.dispMag = dispMag;
        deformationData{end}.maxDisp = maxDisp;
        deformationData{end}.color = materials(mi).color; % basic color for legend grouping
    end
end

if isempty(deformationData)
    warning('No Pareto candidates produced displacement results — nothing to overlay.');
else
    % choose scale for visualization: scale so largest displacement ~ 8% of chassis diagonal
    chassisDiag = norm(max(nodes) - min(nodes));
    allMaxDispVals = cellfun(@(d)d.maxDisp, deformationData);
    allMaxDisp = max(allMaxDispVals);
    if allMaxDisp > 0
        scale = 0.08 * chassisDiag / allMaxDisp;
    else
        scale = 1.0;
    end

    % prepare color mapping based on displacement magnitude (global)
    cmap = jet(256);
    globalMin = inf; globalMax = -inf;

    % compute element displacements for color scaling
    for k = 1:numel(deformationData)
        def = deformationData{k};
        def.elemDisp = zeros(nElems,1);
        for e = 1:nElems
            n1 = elems(e,1); n2 = elems(e,2);
            midDisp = 0.5*(def.Utrans(n1,:) + def.Utrans(n2,:));
            def.elemDisp(e) = norm(midDisp);
        end
        deformationData{k} = def;
        globalMin = min(globalMin, min(def.elemDisp));
        globalMax = max(globalMax, max(def.elemDisp));
    end
    if globalMin == globalMax, globalMin = 0; end

    % plot overlays
    for k = 1:numel(deformationData)
        def = deformationData{k};
        Utrans = def.Utrans;
        elemDisp = def.elemDisp;
        cand = def.candidate;
        matcol = def.color;
        for e = 1:nElems
            n1 = elems(e,1); n2 = elems(e,2);
            p1 = nodes(n1,:) + scale * Utrans(n1,:);
            p2 = nodes(n2,:) + scale * Utrans(n2,:);
            % map color to displacement magnitude
            if globalMax > 0
                tval = (elemDisp(e)-globalMin)/(globalMax-globalMin+eps);
            else
                tval = 0;
            end
            idxc = max(1, min(256, round(1 + tval*255)));
            col = cmap(idxc,:);
            plot3([p1(1) p2(1)], [p1(2) p2(2)], [p1(3) p2(3)], '-', ...
                'Color', col, 'LineWidth', 2.0);
        end
    end

    % Add colorbar reflecting displacement magnitude
    caxis([globalMin globalMax]);
    hcb = colorbar;
    ylabel(hcb, 'Element midpoint displacement (m)');

    % Determine unique material names for legend
    matNamesUnique = unique(cellfun(@(d)d.mat, deformationData, 'UniformOutput', false));
    legendEntries = [{'Undeformed'}, matNamesUnique];
    legend(legendEntries, 'Location', 'bestoutside');

    % Annotate candidates for each material
    yl = ylim; xl = xlim; zl = zlim;
    textPos = [xl(2) + 0.05*(xl(2)-xl(1)), yl(2), zl(1)];
    ypos = yl(2);
    for k = 1:numel(deformationData)
        cand = deformationData{k}.candidate;
        matn = deformationData{k}.mat;
        lab = sprintf('%s: OD=%.3f" w=%.4f" maxDisp=%.3f mm', ...
            matn, cand.OD_in, cand.Wall_in, deformationData{k}.maxDisp*1000);
        ypos = ypos - 0.03*(yl(2)-yl(1));
        text(textPos(1), ypos, textPos(3), lab, 'FontSize', 9, 'Interpreter', 'none');
    end
end

xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)');
title('Undeformed (black) and Pareto candidate deformations (color by displacement)');

%% ---------------- Export results ----------------
writetable(valid_safety, 'baja_tube_analysis_results.csv');
fprintf('\nSaved results to baja_tube_analysis_results.csv\n');

%% ---------------- Helper functions ----------------

function is_p = paretofront(X)
% Simple pareto front (minimize columns)
n = size(X,1);
is_p = true(n,1);
for i=1:n
    if is_p(i)
        dominated = all(X <= X(i,:),2) & any(X < X(i,:),2);
        is_p(dominated) = false;
    end
end
end

function [Kglob, elemBases] = assembleGlobal3D(nodes, elems, section, mat)
% Assemble global stiffness matrix for 3D frame elements (12x12 local -> global)
nNodes = size(nodes,1); nElems = size(elems,1);
ndof = 6*nNodes;
Kglob = zeros(ndof, ndof);
elemBases = cell(nElems,1);
for e=1:nElems
    ni = elems(e,1); nj = elems(e,2);
    x1 = nodes(ni,:); x2 = nodes(nj,:);
    v = x2 - x1; L = norm(v); ex = v / L;
    arbitrary = [0 0 1];
    if abs(dot(arbitrary, ex)) > 0.9, arbitrary = [0 1 0]; end
    ez = cross(ex, arbitrary); ez = ez / norm(ez);
    ey = cross(ez, ex); ey = ey / norm(ey);
    T = [ex; ey; ez]; % 3x3 rotation (local axes in global coords)
    A = section(e).A; Ix = section(e).Ix; Iy = section(e).Iy; J = section(e).J;
    E = mat.E; G = mat.G;
    k_local = frameElementStiffness3D(E,G,A,Iy,Ix,J,L); % 12x12 in local coordinates
    % build 12x12 Lambda transform (block diag of T for each 3-DOF group)
    Lambda = zeros(12,12);
    Lambda(1:3,1:3) = T; Lambda(4:6,4:6) = T; Lambda(7:9,7:9) = T; Lambda(10:12,10:12) = T;
    kglob = Lambda' * k_local * Lambda;
    idxi = (ni-1)*6 + (1:6); idxj = (nj-1)*6 + (1:6);
    Kglob(idxi, idxi) = Kglob(idxi, idxi) + kglob(1:6,1:6);
    Kglob(idxi, idxj) = Kglob(idxi, idxj) + kglob(1:6,7:12);
    Kglob(idxj, idxi) = Kglob(idxj, idxi) + kglob(7:12,1:6);
    Kglob(idxj, idxj) = Kglob(idxj, idxj) + kglob(7:12,7:12);
    elemBases{e}.L = L; elemBases{e}.T = T; elemBases{e}.k_local = k_local; elemBases{e}.nodes = [ni,nj];
end
end

function Kloc = frameElementStiffness3D(E,G,A,Iy,Ix,J,L)
% 12x12 Euler-Bernoulli frame stiffness in element local x-axis coordinates
% DOF order local: [u v w rx ry rz] per node
kax = E*A / L;
k_b_y = 12*E*Iy / L^3;
k_b_z = 12*E*Ix / L^3;
k_b_y_m = 6*E*Iy / L^2;
k_b_z_m = 6*E*Ix / L^2;
k_rot_y = 4*E*Iy / L;
k_rot_z = 4*E*Ix / L;
k_rot_y_c = 2*E*Iy / L;
k_rot_z_c = 2*E*Ix / L;
ktors = G * J / L;
Kloc = zeros(12,12);
% axial
Kloc(1,1)=kax; Kloc(1,7)=-kax; Kloc(7,1)=-kax; Kloc(7,7)=kax;
% bending v (dof2) with rot rz (dof6)
Kloc(2,2)=k_b_y; Kloc(2,6)=k_b_y_m; Kloc(2,8)=-k_b_y; Kloc(2,12)=k_b_y_m;
Kloc(6,2)=k_b_y_m; Kloc(6,6)=k_rot_y; Kloc(6,8)=-k_b_y_m; Kloc(6,12)=k_rot_y_c;
Kloc(8,2)=-k_b_y; Kloc(8,6)=-k_b_y_m; Kloc(8,8)=k_b_y; Kloc(8,12)=-k_b_y_m;
Kloc(12,2)=k_b_y_m; Kloc(12,6)=k_rot_y_c; Kloc(12,8)=-k_b_y_m; Kloc(12,12)=k_rot_y;
% bending w (dof3) with rot ry (dof5)
Kloc(3,3)=k_b_z; Kloc(3,5)=-k_b_z_m; Kloc(3,9)=-k_b_z; Kloc(3,11)=-k_b_z_m;
Kloc(5,3)=-k_b_z_m; Kloc(5,5)=k_rot_z; Kloc(5,9)=k_b_z_m; Kloc(5,11)=k_rot_z_c;
Kloc(9,3)=-k_b_z; Kloc(9,5)=k_b_z_m; Kloc(9,9)=k_b_z; Kloc(9,11)=k_b_z_m;
Kloc(11,3)=-k_b_z_m; Kloc(11,5)=k_rot_z_c; Kloc(11,9)=k_b_z_m; Kloc(11,11)=k_rot_z;
% torsion rx
Kloc(4,4)=ktors; Kloc(4,10)=-ktors; Kloc(10,4)=-ktors; Kloc(10,10)=ktors;
% ensure symmetric
Kloc = (Kloc + Kloc') - diag(diag(Kloc));
end

function [Kff,Kfp,freeDOFs,Ufull] = applyBCandSolve(Kglob,F,fixedDOFs)
ndof = length(F);
freeDOFs = find(~fixedDOFs);
fixDOFs = find(fixedDOFs);
if isempty(freeDOFs)
    Ufull = zeros(ndof,1); Kff=[]; Kfp=[];
    return;
end
Kff = Kglob(freeDOFs, freeDOFs);
Kfp = Kglob(freeDOFs, fixDOFs);
Ff = F(freeDOFs);
% Solve
tol = 1e-12;
% Add tiny regularization if singular
[Ufree,flag] = solveLinear(Kff, Ff);
if flag ~= 0
    Ufull = [];
    return;
end
Ufull = zeros(ndof,1);
Ufull(freeDOFs) = Ufree;
Ufull(fixDOFs) = 0;
end

function [x,flag] = solveLinear(A,b)
% Solve linear system robustly with fallback to pinv
flag = 0;
x = [];
try
    x = A\b;
catch
    % try least squares
    x = pinv(A)*b;
    flag = 0;
end
if any(isnan(x)) || any(isinf(x))
    flag = 1;
end
end

% End of script