%
% Exemple de FEM per elements lineals: LINK
% (2-Dim: despla\c{c}aments en x i y)
%
% (c) Toni Susin, 2012
%-----------------------------------------------------------


%------------------------------------
% Carreguem la geometria: nodes i elements
%------------------------------------
nod=load('ArmaduraHoweNodes.txt','-ascii');
[num_nod,ndim]=size(nod);
%Matriu de connectivitat
B=load('ArmaduraHoweConnectivitat.txt','-ascii');
elem=B; %son la mateixa cosa
[num_elem,ndim_elem]=size(elem);

%------------------------------------
%Pintem els nodes i els elements
%------------------------------------
figure(1)
pintemNodElem(nod,elem);
%------------------------------------
%Propietats del material
%------------------------------------
% m\'{o}dul de Young
%Y=50;
Y=200; %GPa=KN/mm^2
E=Y*ones(1, num_elem); %tots els elements del mateix material
%E(3)=100;
%------------------------------------
%Constants reals: Area de secci\'{o} de la barra
%------------------------------------
%Area=2;

Area=4800; %mm^2
A=Area*ones(1,num_elem);
%A(2)=1;
%A(3)=2*sqrt(2);
%------------------------------------
%Pintem la matriu de rigidesa
%------------------------------------
figure(2)
K=zeros(num_nod*ndim);
for e=1:num_elem
    Ke=PlaneLinkBar_Mrigidesa(nod,elem,e,E,A)
    files=[elem(e,1)*2-1; elem(e,1)*2; elem(e,2)*2-1; elem(e,2)*2];
    colums=files;
    K(files,colums)=K(files,colums)+Ke;
    spy(K)
    pause(0.1);
end
%matriu de rigidesa
K
%------------------------------------
% SOLVE: Apliquem condicions de contorn
%------------------------------------
%Inicialitzem forces
f=zeros(num_nod*ndim,1);
fixed=zeros(1,num_nod*ndim);
val=zeros(1,num_nod*ndim);
fapp=zeros(num_nod*ndim,1);

% Fixem els punts: 1 all DOF, 3
% codifiquem amb un 1 els punts en els que fixem despla�aments
% exemp: fixed=[1,1,0,1,0,0];

%fixed(1*ndim-1)=1; %(u1_x=0)
%fixed(1*ndim)=1; %(u1_y=0)
%fixed(2*ndim)=1; %(u2_y=0)

fixed(1*ndim-1)=1; %(u1_x=0)
fixed(1*ndim)=1;   %(u1_y=0)
fixed(8*ndim)=1;   %(u8_y=0=

% codifiquem les forces i els valors dels despla\c{c}aments
% exemp: val=[0,0,0,0,2,1];
% forces al node 3 : (fx=2, fy=1)

%nodf=3; %node en el que apliquem la for\c{c}a
%val(nodf*ndim-1)=2;
%val(nodf*ndim)=1;
%fapp(nodf*ndim-1)=2;
%fapp(nodf*ndim)=1;

nodf=1;
val(nodf*ndim-1)=0; %u1_x=0;
val(nodf*ndim)=0;   %u1_y=0;
nodf=2;
val(nodf*ndim-1)=30; %KN
val(nodf*ndim)=-50;  %kN
nodf=4;
val(nodf*ndim-1)=0;  %KN
val(nodf*ndim)=-40;  %KN
nodf=6;
val(nodf*ndim-1)=0;  %KN
val(nodf*ndim)=-50;  %KN
nodf=8;
val(nodf*ndim)=0;   %u8_y=0;

nodf=2;
fapp(nodf*ndim-1)=30; %KN
fapp(nodf*ndim)=-50;  %KN
nodf=4;
fapp(nodf*ndim-1)=0;  %KN
fapp(nodf*ndim)=-40;  %KN
nodf=6;
fapp(nodf*ndim-1)=0;  %KN
fapp(nodf*ndim)=-50;  %KN

[Km,fm]=ImposemCC(K,f,fixed,val);
Km
fm
%------------------------------------
% SOLVE: Resolem el sistema i calculem els despla�aments
%------------------------------------

u=Km\fm

%------------------------------------
% POSTPROCES: Pintem despla\c{c}aments i calculem forces i tensions
%------------------------------------

% Representem la soluci\'{o} despla\c{c}ada
figure(3)
esc=100; %escala per representar els despla\c{c}aments
pintemNodElemDespl(nod,elem,u,esc);
%Forces de Reacci\'{o}
fr=K*u-fapp
%Tensions Axials
[T, F]=calculTensionsAxials(nod,elem,u,E,A)

%---------------------------------------------------------------
% POSTPROCES: Sortida de les solucions
%---------------------------------------------------------------
%fout=1;
fout=fopen('Howe_sol.txt','w');
fprintf(fout,'\nSOLUTION:\n');
%fprintf(fout,'UNITS: UX,UY(mm), FRX,FRY(N), F(N), T(N/mm^2)');
fprintf(fout,'====================================================\n');
fprintf(fout,'NODAL SOLUTION PER NODE:\n');
fprintf(fout,'%4s %23s %23s\n','NODE','UX (mm)','UY (mm)');
for i=2:ndim:length(u);
    fprintf(fout,'%4d %23.15E %23.15E\n',i/ndim,u(i-1),u(i));
end
fprintf(fout,'\nREACTION SOLUTIONS PER NODE:\n');
fprintf(fout,'%4s %23s %23s\n','NODE','FX (kN)','FY (kN)');
for i=2:ndim:length(fr);
    fprintf(fout,'%4d %23.15E %23.15E\n',i/ndim,fr(i-1),fr(i));
end
fprintf(fout,'\nAXIAL STRESS:\n');
fprintf(fout,'%4s %23s\n','ELEM','TAXL (GPa)');
for i=1:length(F)
    fprintf(fout,'%4d %23.15E\n',i,T(i));
end
fprintf(fout,'====================================================\n');
%Despla\c{c}aments en els punts F i H
numNodeE=5; %num.node E = 5
numNodeH=8; %num.node H = 8
fprintf(fout,'\nDISPLACEMENTS AT POINTS E AND H (mm):\n');
fprintf(fout,'%30s %23s\n','E','H');
fprintf(fout,'%6s %23.15E %23.15E\n','U(x)',u(numNodeE*ndim-1),u(numNodeH*ndim-1));
fprintf(fout,'%6s %23.15E %23.15E\n','U(y)',u(numNodeE*ndim),u(numNodeH*ndim));

%Forces de reaccio en els punts A i H
numNodeA=1; %num.node A = 1
numNodeH=8; %num.node H = 8
fprintf(fout,'\nREACTION FORCES AT POINTS A AND H (kN):\n');
fprintf(fout,'%30s %23s\n','A','H');
fprintf(fout,'%6s %23.15E %23.15E\n','FX',fr(numNodeA*ndim-1),fr(numNodeH*ndim-1));
fprintf(fout,'%6s %23.15E %23.15E\n','FY',fr(numNodeA*ndim),fr(numNodeH*ndim));

%Tensions axials als elements AC i FH
numElemAC=2;  %num.element E = 2
numElemFH=13; %num.element H = 13
fprintf(fout,'\nAXIAL STRESS ON ELEMENTS AC AND FH (GPa):\n');
fprintf(fout,'%30s %23s\n','AC','FH');
fprintf(fout,'%6s %23.15E %23.15E\n','TAXL',T(numElemAC),T(numElemFH));

% Matriu de Rigidesa de l'element BE
numElemBE=5; % L'element BE és l'element 5 segons la nostra numeració
fprintf(fout,'\nSTIFFNESS MATRIX OF ELEMENT BE:\n')
Ke=PlaneLinkBar_Mrigidesa(nod,elem,numElemBE,E,A);
for i=1:4
    for j=1:4
        fprintf(fout,'%25.15e',Ke(i,j));
    end
    fprintf(fout,'\n');
end
fprintf(fout,'\nTRACE OF THE STIFFNESS MATRIX OF ELEMENT BE:\n  %23.15E',trace(Ke));
fclose(fout);
