function [g_mat,g_stab,V_4,V_8] = portalFrame_2018_08(E,sig_y,L,H,elem,b,h,F_V,F_B,F_H,LV,HH,comment)

%% 2.1. Elastic system stiffness matrix K

% top left (oben links)
K_ol = [E*b(1)*h(1)^3/HH^3+E*b(1)*h(1)^3/(H-HH)^3   0                                   E*b(1)*h(1)^3/(2*HH^2)-E*b(1)*h(1)^3/(2*(H-HH)^2)   -E*b(1)*h(1)^3/(H-HH)^3                 0                                       -E*b(1)*h(1)^3/(2*(H-HH)^2);
        0                                           E*b(1)*h(1)/HH+E*b(1)*h(1)/(H-HH)   0                                                   0                                       -E*b(1)*h(1)/(H-HH)                     0;
        0                                           0                                   E*b(1)*h(1)^3/(3*HH)+E*b(1)*h(1)^3/(3*(H-HH))       E*b(1)*h(1)^3/(2*(H-HH)^2)              0                                       E*b(1)*h(1)^3/(6*(H-HH));
        0                                           0                                   0                                                   E*b(1)*h(1)^3/(H-HH)^3+E*b(3)*h(3)/LV   0                                       E*b(1)*h(1)^3/(2*(H-HH)^2);
        0                                           0                                   0                                                   0                                       E*b(1)*h(1)/(H-HH)+E*b(3)*h(3)^3/LV^3   -E*b(3)*h(3)^3/(2*LV^2);
        0                                           0                                   0                                                   0                                       0                                       E*b(1)*h(1)^3/(3*(H-HH))+E*b(3)*h(3)^3/(3*LV)];
K_ol = K_ol+tril(K_ol',-1);

% top right (oben rechts)
K_or = [0                   0                       0                       0   0   0;
        0                   0                       0                       0   0   0;
        0                   0                       0                       0   0   0;
        -E*b(3)*h(3)/LV     0                       0                       0   0   0;
        0                   -E*b(3)*h(3)^3/LV^3     -E*b(3)*h(3)^3/(2*LV^2) 0   0   0;
        0                   E*b(3)*h(3)^3/(2*LV^2)  E*b(3)*h(3)^3/(6*LV)    0   0   0];

% bottom left (unten links)
K_ul = K_or';

% bottom right (unten rechts)
K_ur = [E*b(3)*h(3)/LV+E*b(3)*h(3)/(L-LV)   0                                           0                                                   -E*b(3)*h(3)/(L-LV)                     0                                       0;
        0                                   E*b(3)*h(3)^3/LV^3+E*b(3)*h(3)^3/(L-LV)^3   E*b(3)*h(3)^3/(2*LV^2)-E*b(3)*h(3)^3/(2*(L-LV)^2)   0                                       -E*b(3)*h(3)^3/(L-LV)^3                 -E*b(3)*h(3)^3/(2*(L-LV)^2);
        0                                   0                                           E*b(3)*h(3)^3/(3*LV)+E*b(3)*h(3)^3/(3*(L-LV))       0                                       E*b(3)*h(3)^3/(2*(L-LV)^2)              E*b(3)*h(3)^3/(6*(L-LV));
        0                                   0                                           0                                                   E*b(3)*h(3)/(L-LV)+E*b(5)*h(5)^3/H^3    0                                       E*b(5)*h(5)^3/(2*H^2);
        0                                   0                                           0                                                   0                                       E*b(3)*h(3)^3/(L-LV)^3+E*b(5)*h(5)/H    E*b(3)*h(3)^3/(2*(L-LV)^2);
        0                                   0                                           0                                                   0                                       0                                       E*b(3)*h(3)^3/(3*(L-LV))+E*b(5)*h(5)^3/(3*H)];
K_ur = K_ur+tril(K_ur',-1);

% in total
K = [K_ol   K_or;
     K_ul   K_ur];

%%  2.2. Load vector P
P = zeros(12,1);
P(1) = F_H;
P(7) = F_B;
P(8) = F_V;

%% 2.3. Calculate global nodal displacements V
V = K\P;
% V = inv(K)*P;

% outputs for data assimilation
V_4 = V(4);
V_8 = V(8);

%% 2.4. Generate local element displacements
v_e = zeros(6,elem);

v_e(:,1) = [0       0       0       -V(2)   V(1)    V(3)]';
v_e(:,2) = [-V(2)   V(1)    V(3)    -V(5)   V(4)    V(6)]';
v_e(:,3) = [V(4)    V(5)    V(6)    V(7)    V(8)    V(9)]';
v_e(:,4) = [V(7)    V(8)    V(9)    V(10)   V(11)   V(12)]';
v_e(:,5) = [V(11)   -V(10)  V(12)   0       0       0]';

%% 2.5. Calculate internal forces of each element

% cross-sectional parameters
A_e = zeros(elem,1);
I_e = zeros(elem,1);
W_e = zeros(elem,1);

for i=1:elem
    A_e(i) = b(i)*h(i);
    I_e(i) = b(i)*h(i)^3/12;
    W_e(i) = b(i)*h(i)^2/6;
end

% element lengths
l_e = zeros(elem,1);
l_e(1) = HH;
l_e(2) = H-HH;
l_e(3) = LV;
l_e(4) = L-LV;
l_e(5) = H;

% elastic element stiffness matrices
k_e = zeros(6,6,elem);
for i=1:elem
    
    k_e(:,:,i) = [  E*A_e(i)/l_e(i) 0                       0                       -E*A_e(i)/l_e(i)    0                       0;
                    0               12*E*I_e(i)/l_e(i)^3    -6*E*I_e(i)/l_e(i)^2    0                   -12*E*I_e(i)/l_e(i)^3   -6*E*I_e(i)/l_e(i)^2;
                    0               0                       4*E*I_e(i)/l_e(i)       0                   6*E*I_e(i)/l_e(i)^2     2*E*I_e(i)/l_e(i);
                    0               0                       0                       E*A_e(i)/l_e(i)     0                       0;
                    0               0                       0                       0                   12*E*I_e(i)/l_e(i)^3    6*E*I_e(i)/l_e(i)^2;
                    0               0                       0                       0                   0                       4*E*I_e(i)/l_e(i)];
    k_e(:,:,i) = k_e(:,:,i)+tril(k_e(:,:,i)',-1);
end

% internal forces
s_e = zeros(6,elem);
for i=1:elem
    s_e(:,i) = k_e(:,:,i)*v_e(:,i);
    
    for j=1:3
        s_e(j,i) = -s_e(j,i);
    end
end

%% 2.6. Calculate maximum stresses

sig_max = zeros(2,elem);

for i=1:elem
    for j=1:2
        sig_max(j,i) = abs(s_e(3*j-2,i)/A_e(i))+abs(s_e(3*j,i)/W_e(i));
    end
end
[max_sig,ort_lr] = max(sig_max);
[max_sig,ort_elem] = max(max_sig);

if comment~=0
    if ort_elem ==1||ort_elem==2
        fprintf('%s','Maxmium stress in left column');
        if ort_elem==1 && ort_lr(ort_elem)==1
            fprintf('%s',' at the bottom');
        elseif ort_elem==2 && ort_lr(ort_elem)==2
            fprintf('%s',' on top');
        else
            fprintf('%s',' at F_H');
        end
    elseif ort_elem ==3||ort_elem==4
        fprintf('%s','Maximum stress in girder');
        if ort_elem==3 && ort_lr(ort_elem)==1
            fprintf('%s',' left');
        elseif ort_elem==4 && ort_lr(ort_elem)==2
            fprintf('%s',' right');
        else
            fprintf('%s',' at F_V and F_B');
        end
    else
        fprintf('%s','Maximum stress in right column');
        if ort_lr(ort_elem)==1
            fprintf('%s',' on top');
        else
            fprintf('%s',' at the bottom');
        end
    end
    fprintf('%s\n\n','.');
end

%% 2.7. Limit state function value (material failure)

g_mat = sig_y/max_sig-1.0;

%% 2.8. Linearized geometric system stiffness matrix K_g

% normal forces
N = zeros(elem,1);
for i=1:elem
    N(i) = s_e(4,i);
end

% top left (oben links)
Kg_ol = [   6*N(1)/(5*HH)+6*N(2)/(5*(H-HH)) 0   N(1)/10-N(2)/10                 -6*N(2)/(5*(H - HH))    0               -N(2)/10;
            0                               0   0                               0                       0               0;
            0                               0   2*N(1)*HH/15+2*N(2)*(H-HH)/15   N(2)/10                 0               -N(2)*(H-HH)/30;
            0                               0   0                               6*N(2)/(5*(H-HH))       0               N(2)/10;
            0                               0   0                               0                       6*N(3)/(5*LV)   -N(3)/10;
            0                               0   0                               0                       0               2*N(2)*(H-HH)/15+2*N(3)*LV/15];
Kg_ol = Kg_ol+tril(Kg_ol',-1);

% top right (oben rechts)
Kg_or = [   0   0               0           0   0   0;
            0   0               0           0   0   0;
            0   0               0           0   0   0;
            0   0               0           0   0   0;
            0   -6*N(3)/(5*LV)  -N(3)/10    0   0   0;
            0   N(3)/10         -N(3)*LV/30 0   0   0];

% bottom left (unten links)
Kg_ul = Kg_or';

% bottom right (unten rechts)
Kg_ur = [   0   0                               0                               0               0                   0;
            0   6*N(3)/(5*LV)+6*N(4)/(5*(L-LV)) N(3)/10-N(4)/10                 0               -6*N(4)/(5*(L-LV))  -N(4)/10;
            0   0                               2*N(3)*LV/15+2*N(4)*(L-LV)/15   0               N(4)/10             -N(4)*(L-LV)/30;
            0   0                               0                               6*N(5)/(5*H)    0                   N(5)/10;
            0   0                               0                               0               6*N(4)/(5*(L-LV))   N(4)/10;
            0   0                               0                               0               0                   2*N(4)*(L-LV)/15+2*N(5)*H/15];
Kg_ur = Kg_ur+tril(Kg_ur',-1);

% in total
Kg = [  Kg_ol   Kg_or;
        Kg_ul   Kg_ur];

%% 2.9. Solving eigenvalue problem

[~,lambda] = eig(K,-Kg);

%% 2.10. Load factor for stability failure

lambda_stab = min(abs(diag(lambda)));

%% 2.11. Limit state function value (stability failure)

g_stab = lambda_stab-1.0;

end
