%MATLAB M-file pro prvky s linearni aproximaci posunu
%lokalizace pro bi-linearni prac. diagram
%mekci prvek je uprostred tyce
%23.10.2005 vit.smilauer@fsv.cvut.cz


%|-------------------------   ---->F

%|--> x

clear all; 
echo off;

%input data=======================================================

E = 210.e+6;  %ocel [kPa]
E_red_coef = 0.001; %reduction of Young's modulus
epsilon_zlom = 0.001;
l = 10.; %delka tyce
A = 0.001; %plocha rezu [m2]
P = 210; %max. koncova sila v kN
pocet_kroku = 1000; %pocet rovnomernych kroku pritizeni do max. sily P
redukce = 0.99; %redukce plochy prutu pro lokalizaci
pocet_prvku = 3; %nejlepe lichy pocet, pak je mekci prvek uprostred

%end of input data================================================

L = l/pocet_prvku; %delka prvku

%linearni aproximace
k = E*A/L*...
[1 -1;
 -1 1];

%iniciace
r_g=zeros(pocet_prvku,1); %redukovany vektor posunu
F_g=zeros(pocet_prvku,1); %redukovany vektor zatizeni

for cyklus=1:1:pocet_kroku
  K_g=zeros(pocet_prvku+1);
  delta_F_g=zeros(pocet_prvku,1);
  
  %lokalizace matic tuhosti (tecna matice tuhosti)
  for i=1:1:(pocet_prvku)
    if i==1
      epsilon = r_g(1,1)/L;
    else
      epsilon = (r_g(i,1)-r_g(i-1,1))/L;
    end
    epsilon   
        
    if epsilon < epsilon_zlom
       k_help = k;
    else
       k_help = E_red_coef * k;
    end
    
    if i == floor((pocet_prvku+1)/2);%prostredni prvek - redukce plochy
       k_help = redukce * k_help;
    end
  
    K_g(i:i+1,i:i+1)=K_g(i:i+1,i:i+1)+k_help(1:2,1:2);
  end;

  %redukce globalni matice tuhosti vlivem vetknuti (zrusit prvni radek a
  %sloupec)
  K_red=zeros(pocet_prvku);
  K_red(1:pocet_prvku,1:pocet_prvku) = K_g(2:pocet_prvku+1,2:pocet_prvku+1);

  delta_F_g(pocet_prvku,1)=P/pocet_kroku;
  %reseni soustavy
  delta_r_g = inv(K_red) * delta_F_g;
  
  %tisk posunu koncu
  delta_r_g
  
  %pricteni prirustku
  r_g = r_g + delta_r_g
  F_g = F_g + delta_F_g
end;

echo on

