% ------------------------------------------------------------------------------------------------------------------ #
%
%   DYNARDO example (c) Dynardo 2012
%
%   ten bar truss (Haftka & Guerdal 1992)
%
%   Minimize the mass
%   0.1 <= a_i <= 20 in^2
%   Maximum stress is 25000
%
% ------------------------------------------------------------------------------------------------------------------ #

% Check if MATLAB is called outside optiSLang
if ~exist('OSL_REGULAR_EXECUTION','var')
	clear all;
	close all;
    OSL_REGULAR_EXECUTION = false;
end

% Define cross section areas
if OSL_REGULAR_EXECUTION == false 
    area01 = 10.00000; % Optimal value = 7.94
    area02 = 10.00000; % Optimal value = 0.10
    area03 = 10.00000; % Optimal value = 8.06
    area04 = 10.00000; % Optimal value = 3.94
    area05 = 10.00000; % Optimal value = 0.10
    area06 = 10.00000; % Optimal value = 0.10
    area07 = 10.00000; % Optimal value = 5.74
    area08 = 10.00000; % Optimal value = 5.57
    area09 = 10.00000; % Optimal value = 5.57
    area10 = 10.00000; % Optimal value = 0.10
end       

% Material properties
Youngs  = 10000000; 
density = 0.100000; 

% Loads
F1 = -100000.000; 
F2 = -100000.000; 

% Length
L = 360;

% Calculate mass
mass=L*density*(area01+area02+area03+area04+area05+area06+sqrt(2)*(area07+area08+area09+area10))

% Define global vectors, dofs, incidences
areas=[area01,area02,area03,area04,area05,area06,area07,area08,area09,area10];
node_coordinates=[2*L,L; 2*L,0; L,L; L,0; 0,L; 0,0];
element_incidences=[3,5; 1,3; 4,6; 2,4; 3,4; 1,2; 4,5; 3,6; 2,3; 1,4];
num_nodes = 6;
num_elements = 10;
available_dof = [1,2,3,4,5,6,7,8];
num_available_dof = 8;
global_load = [0,0, 0,F2, 0,0, 0,F1, 0,0, 0,0];

% Remove node 1 from dof list, if unimportant trusses are set to zero
if area02 == 0 && area06 == 0 && area10 == 0
    available_dof = [3,4,5,6,7,8];
    num_available_dof = 6;
end

% Build up global stiffness
global_stiffness = zeros(2*num_nodes,2*num_nodes);
for i=1:num_elements
    node1 = element_incidences(i,1);
    node2 = element_incidences(i,2);
    coor_node1 = node_coordinates(node1,:);
    coor_node2 = node_coordinates(node2,:);
    area = areas(i);

    % Local element stiffness
    element_len = sqrt((coor_node1(1)-coor_node2(1))^2+(coor_node1(2)-coor_node2(2))^2);
    local_stiffness = Youngs*area/element_len*[1,-1;-1,1];

    % Transform local stiffness to global coordinates
    angle = atan2(coor_node2(2)-coor_node1(2),coor_node2(1)-coor_node1(1));
    trafo_matrix = [cos(angle),sin(angle),0,0; 0,0,cos(angle),sin(angle)];
    element_stiffness=trafo_matrix'*local_stiffness*trafo_matrix;

    % Assemble element matrix in global matrix
    nodes=[node1,node2];
    for k=1:2
        for l=1:2
            pos1=2*nodes(k)-1;
            pos2=2*nodes(l)-1;
            global_stiffness(pos1,  pos2  ) = global_stiffness(pos1,  pos2  )+element_stiffness(2*k-1,2*l-1);
            global_stiffness(pos1,  pos2+1) = global_stiffness(pos1,  pos2+1)+element_stiffness(2*k-1,2*l  );
            global_stiffness(pos1+1,pos2  ) = global_stiffness(pos1+1,pos2  )+element_stiffness(2*k  ,2*l-1);
            global_stiffness(pos1+1,pos2+1) = global_stiffness(pos1+1,pos2+1)+element_stiffness(2*k  ,2*l  );
        end
    end
end

% Reduced global stiffness and loads according to free dofs
global_stiffness_help = global_stiffness(:,available_dof);
global_stiffness_reduced = global_stiffness_help(available_dof,:);
global_load_reduced = global_load(available_dof);

% Solve linear equation, extract global displacements
global_disp_reduced = global_load_reduced/global_stiffness_reduced;
global_disp = zeros(1,num_nodes*2);
global_disp(available_dof)=global_disp_reduced

% Calculate stresses of elements
for i=1:num_elements
    node1 = element_incidences(i,1);
    node2 = element_incidences(i,2);
    coor_node1 = node_coordinates(node1,:);
    coor_node2 = node_coordinates(node2,:);
    area = areas(i);

    % Transform global disp in element coordinates
    element_len = sqrt((coor_node1(1)-coor_node2(1))^2+(coor_node1(2)-coor_node2(2))^2);
    angle = atan2(coor_node2(2)-coor_node1(2),coor_node2(1)-coor_node1(1));
    trafo_matrix = [cos(angle),sin(angle),0,0; 0,0,cos(angle),sin(angle)];
    disp_matrix=[global_disp(2*node1-1), global_disp(2*node1), global_disp(2*node2-1), global_disp(2*node2)];
    disp_local = trafo_matrix*disp_matrix';

    % Calculate local stresses
    stress = Youngs*(disp_local(2)-disp_local(1))/element_len;
    element_stresses(i) = stress;
    if area == 0
        element_stresses(i) = 0;
    end

end
element_stresses

max_disp=max(abs(global_disp))
max_stress=max(abs(element_stresses))
