clear all
close all

E=25e6;
ni=0.16;
Nodes = 4;
Nel = 2;
t = 0.2;
Ndof = 2*Nodes;   

Coord = [0 1.5; 0 0; 2 0.5; 2 1.5];
bc = [1 0; 2 0; 3 0; 4 0];
Inc = [1 2 3; 1 3 4];

for i=1:Nel
    n1 = Inc(i,1);
    n2 = Inc(i,2);
    n3 = Inc(i,3);
    Ex(i,1:3)=[Coord(n1,1),Coord(n2,1),Coord(n3,1)];
    Ey(i,1:3)=[Coord(n1,2),Coord(n2,2),Coord(n3,2)];
    Edof(i,1:7)=[i, n1*2-1, n1*2,   n2*2-1, n2*2,  n3*2-1, n3*2];
end

figure(1);
eldraw2(Ex,Ey,[1 1 1],Edof(:,1));
hold on;
grid on;
axis equal;

KG = zeros(Ndof,Ndof);
F = zeros(Ndof,1);

D=E/(1-ni^2)*[1 ni 0; ni 1 0; 0 0 (1-ni)/2];  
for i=1:Nel
     Kel=plante(Ex(i,:),Ey(i,:),[1 t],D);
     KG=assem(Edof(i,:),KG,Kel);
end
F(2) = -25;
F(8) = -50;

[u,r]=solveq(KG,F,bc)

Ed = extract(Edof,u);
sfac=scalfact2(Ex,Ey,Ed,0.2);
eldisp2(Ex,Ey,Ed,[2 2 1],sfac);
hold on
grid on
title('Przemieszczenie [m]');

Sig = zeros(Nel,3);
Eps = zeros(Nel,3);
for i=1:Nel
     [s,e]=plants(Ex(i,:),Ey(i,:),[1 t],D,Ed(i,:))
     Sig(i,:) = s;
     hmh(i) = (1/sqrt(2)) * sqrt( (s(1) - s(2))^2 + 6 * s(3)^2 );
     Eps(i,:) = e;
 end
 maxs = max(abs(Sig));
 maxe = max(abs(Eps));
 figure(2);
 naz1 = ['Sxx';'Syy';'Sxy'];
 naz2 = ['Exx';'Eyy';'Exy'];
 for i=1:Nel
     for j=1:3
         subplot(2,3,j);
         fill(Ex(i,:),Ey(i,:),(1-abs(Sig(i,j)/maxs(j)))*[1 1 1]);
         hold on
         grid on
         title([naz1(j,:) '[kN/m^2]']);
         subplot(2,3,3+j);
         fill(Ex(i,:),Ey(i,:),(1-abs(Eps(i,j)/maxe(j)))*[1 1 1]);
         hold on
         grid on
         title([naz2(j,:) ' [-]']);
     end
end

figure(3);
maxh = max(abs(hmh));
for i=1:Nel
    fill(Ex(i,:),Ey(i,:),(1-abs(hmh(i)/maxh))*[1 1 1]);
    hold on
    grid on
    title('Huber-Mises-Hencky (von Mises) [kN/m^2]');
end
    

