%
% 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('nodesTorre.txt','-ascii');
[num_nod,ndim]=size(nod);
%Matriu de connectivitat
B=load('connectivitatTorre.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�dul de Young
Y=1100; %N/mm^2
E=Y*ones(1, num_elem); %tots els elements del mateix material
E(1:5)=350; %N/mm^2
%------------------------------------
%Constants reals: Area de secci\'o de la barra
%------------------------------------
Area=9000; %mm^2
A=Area*ones(1,num_elem);
A(1:5)=4500; %mm^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 i 2, all DOF 
% codifiquem amb un 1 els punts en els que fixem despla\c{c}aments
% exemp: fixed=[1,1,0,1,0,0];
fixed(1*ndim-1)=1; %Node 1 en la direcci\'o x
fixed(1*ndim)=1;   %Node 1 en ladirecci\'o y
fixed(2*ndim-1)=1; %Node 2 en la direcci\'o x
fixed(2*ndim)=1;   %Node 2 en la direcci\'o y
% 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)

% Posem primer **totes** les forces aplicades al vector val 

nodf=13 % forces aplicades al node 13:
val(nodf*ndim-1)= 1112 %N en la direcci\'o x
val(nodf*ndim)= 0; %N en la direcci\'o y 

nodf=14 % forces aplicades al node 14:
val(nodf*ndim-1)= 0 %N en la direcci\'o x
val(nodf*ndim)= -4448; %N en la direcci\'o y 

nodf=17 % forces aplicades al node 17:
val(nodf*ndim-1)= -1112; %N en la direcci\'o x
val(nodf*ndim)= 0; % kN en la direcci\'o y  

nodf=18 % forces aplicades al node 18:
val(nodf*ndim-1)= 0; %N en la direcci\'o x
val(nodf*ndim) = -4448 %N en la direcci\o y

% Copiem les forces aplicades al vector fapp:
fapp=val';

% Fixem el valors dels despla\c{c}aments fixats al vector val
nodfix=1;
val(nodfix*ndim-1)=0; %despla\c{c}ament en x del node 1 fixat a 0;
val(nodfix*ndim)=0; %despla\c{c}ament en y del node 1 fixat a 0;

nodfix=2;
val(nodfix*ndim-1)=0; %despla\c{c}ament en x del node 2 fixat a 0;
val(nodfix*ndim)=0; %despla\c{c}ament en y del node 2 fixat a 0;

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

u=Km\fm

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

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

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sortida:
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fout = fopen('Solucions_Pont.txt','w');
fout=1;
fprintf(fout,'------------------------------------------------------\n');
fprintf(fout,'DESPLAZAMIENTOS:\n');
fprintf(fout,'%4s %24s %24s\n','NODO','UX','UY');
for i=1:num_nod
    fprintf(fout,'%4d %24.16e %24.16e\n',i,u(i*ndim-1),u(i*ndim));
end
fprintf(fout,'\nTENSIONES AXIALES Y FUERZAS SOBRE CADA ELEMENTO:\n');
fprintf(fout,'%4s %24s %24s\n','elem','T','F');
for i=1:num_elem
    fprintf(fout,'%4d %24.16e %24.16e\n',i,T(i),F(i));
end
%fclose(fout);
%Maximo desplazamiento en Y, caso modificado: 
maxy=max(abs(u(2:2:end)));
fprintf('\n\nDESPLAZAMIENTO MAXIMO EN Y, CASO MODIFICADO: %24.16e\n',maxy);   

