Untitled

ui = input('Is this plane stress (Type 1) or plane strain (Type 2)');

E = 10 * 10^9;
v = 0.3; % Poissons ratio 
L2 = 0.1*cos(63 * (pi/180));
L1 = 0.14 - L2; % m
L3 = 0.1; % m

% Create the proper elastic modulus tensor
if ui == 1
    
    C = (E/(1 - v^2)) * [1, v, 0; v, 1, 0; 0, 0, (0.5*(1-v))];
    % Plane stress elastic mod
elseif ui == 2

    C = (E/((1+v)*(1-2*v))) * [1-v, v, 0; v, 1-v, 0; 0, 0, (0.5*(1-2*v))];
    % Plane strain elastic mod
end

% Return the stiffness matrix for a CST (Constant strain triangle) element
% given inputs L1 L2 L3 E v

D = [1/(L1 + L2), 0, -1/(L1 + L2), 0, 0, 0;
    0, L2 / ((L1 + L2)*L3), 0, -1 / ((L1 + L2)*L3), 0, 1/L3;
    L2 / ((L1 + L2)*L3), 1/(L2 + L3), -1 / ((L1 + L2)*L3), -1 / ((L1 + L2)*L3), 1/L3, 0];

K = (L1 + L2)*L3*0.5 * D * D.' * C * D