% CÁLCULO ESTÁTICO DE PÓRTICOS PLANOS

clear;
fprintf('%s\n','Cálculo Estático de Pórticos Planos');

%  Leitura de Dados
fprintf('%s\n','Leitura dos dados...');
[arqfinal,patharq]=uigetfile('*.txt','Arquivo de Leitura');
%  Abre o arquivo
fid = fopen(arqfinal,'rt');
%  Lê a primeira linha do arquivo
dados = fscanf(fid,'%d',[5]);
%  Obtenção do número de nós:
nnos=dados(1);
%  Cálculo no número de graus de liberdade:
ngl=nnos*3;
%  Obtenção do número de elementos:
nelm=dados(2);
%  Obtenção do número de nós restringidos:
nnr=dados(3);
%  Obtenção do número de nós carregados:
nnc=dados(4);
%  Obtenção do número de vigas carregadas:
nvc=dados(5);

coord = fscanf(fid,'%10e',[2,nnos]);
%  Definição das coordenadas x dos nós: x(i)=abcissa do nó i
x=coord(1,:);
%  Definição das coordenadas y dos nós: y(i)=ordenada do nó i
y=coord(2,:);
prop = fscanf(fid,'%10e',[5,nelm]);
%  Definição do vetor nó1: no1(i)=1º nó que define o elemento i
no1=prop(1,:);
%  Definição do vetor nó2: no2(i)=2º nó que define o elemento i
no2=prop(2,:);
%  Definição do vetor área: area(i)=área do elemento i
area=prop(3,:);
%  Definição do vetor mom inércia
mom=prop(4,:);
%  Definição do vetor módulo de elasticidade: mod(i)=mód. elast. do elemento i
mod=prop(5,:);
restr = fscanf(fid,'%10e',[4,nnr]);
%  Definição do vetor restrições IC:
%  se IC(i)=0 o GL i é livre; se IC(i)=1 o GL i está restringido
IC(nnos*3)=0;
for i=1:nnr,
   nn=restr(1,i);
   IC((nn-1)*3+1)=restr(2,i);
   IC((nn-1)*3+2)=restr(3,i);
   IC((nn-1)*3+3)=restr(4,i);
end
forcas = fscanf(fid,'%10e',[4,nnc]);
%  Definição do vetor de forças globais F:
F(ngl)=0.0;
for i=1:nnc,
   n=forcas(1,i);
   F((n-1)*3+1)=forcas(2,i);
   F((n-1)*3+2)=forcas(3,i);
   F((n-1)*3+3)=forcas(4,i);
end
%  Definição das vigas carregadas:
carreg_vig = fscanf(fid,'%10e',[2,nvc]);
%  fecha o arquivo de dados:
st = fclose(fid);

%  Plota a estrutura (configuração inicial)
fprintf('%s\n','Desenhando a estrutura...');
figure(1);
clf;
for i=1:nelm,
   line([x(no1(i)) x(no2(i))],[y(no1(i)) y(no2(i))]);
end
axis equal;
axis off;
title('Estrutura em Estudo');
for i=1:nelm,
   text((x(no1(i))+x(no2(i)))/2,(y(no1(i))+y(no2(i)))/2,num2str(i));
end
for i=1:nnos,
   texto=text(x(i),y(i),num2str(i));
   set(texto,'Color','red')
end

%  Calcula a matriz de rigidez global pelo processo da rigidez direta
fprintf('%s\n','Calculando a Matriz de Rigidez Global...');
Kg(ngl,ngl)=0.0;
Kg(:,:)=0.0;
S0(nelm,6)=0.0;
Sl(:,:)=0.0;

% Início do Loop  
for nel=1:nelm,
   kk(1)=no1(nel);
   kk(2)=no2(nel);
   % Cálculo da matriz de rotação
   dx=x(kk(2))-x(kk(1));
   dy=y(kk(2))-y(kk(1));
   % Cálculo do comprimento do elemento
   L=sqrt(dx^2+dy^2);
   % Cálculo dos cossenos diretores
   cx=dx/L;
   cy=dy/L;
   % Definição da Matriz de Rotação
   R=[ cx  cy  0  0  0  0;
      -cy  cx  0  0  0  0;
      0    0   1  0  0  0
      0    0   0  cx cy 0;
      0    0   0 -cy cx 0;
      0    0   0  0  0  1];
   % Definição da Matriz de Rigidez do elemento de treliça Ke
   for i=1:nvc,
   ne=carreg_vig(1,i);
   if ne==nel
       q=carreg_vig(2,i)
       S0(ne,:)=R*[0 q*dx/2 q*dx^2/12 0 q*dx/2 -q*dx^2/12]';
       %S0(ne,1)=0;
       %S0(ne,2)=q*L/2;
       %S0(ne,3)=q*L^2/12;
       %S0(ne,4)=0;
       %S0(ne,5)=q*L/2;
       %S0(ne,6)=-q*L^2/12;
       L2=L*abs(cx);
       F((no1(ne)-1)*3+2)=F((no1(ne)-1)*3+2)-q*L2/2;
       F((no1(ne)-1)*3+3)=F((no1(ne)-1)*3+3)-q*L2^2/12;
       F((no2(ne)-1)*3+2)=F((no2(ne)-1)*3+2)-q*L2/2;
       F((no2(ne)-1)*3+3)=F((no2(ne)-1)*3+3)+q*L2^2/12;
    end
   end


   E=mod(nel);
   A=area(nel);   
   J=mom(nel);
   Ke=[E*A/L  0          0          -E*A/L 0           0;
      0      12*E*J/L^3  6*E*J/L^2  0      -12*E*J/L^3 6*E*J/L^2;
      0      6*E*J/L^2   4*E*J/L    0      -6*E*J/L^2  2*E*J/L;
      -E*A/L 0           0          E*A/L  0           0;
      0      -12*E*J/L^3 -6*E*J/L^2 0      12*E*J/L^3 -6*E*J/L^2;
      0      6*E*J/L^2   2*E*J/L    0      -6*E*J/L^2  4*E*J/L];
   % Cálculo da Matriz rotacionada
   RKeR=R'*Ke*R;
   
   % Acumula na matriz de rigidez global Kg
   for i=1:2,
      % Cálculo do número de GLs existentes antes dos GLs relativos ao iºnó
      m=3*(kk(i)-1);
      for j=1:2,
      	% Cálculo do número de GLs existentes antes dos GLs relativos ao iºnó
         n=3*(kk(j)-1);
         k1=0;
         for k=3*i-2:3*i,
            % GL global k1 referente à coor local k
            k1=k1+1;
            l1=0;
            for l=3*j-2:3*j,
               % GL global l1 referente à coor local l
               l1=l1+1;
               Kg(m+k1,n+l1)=Kg(m+k1,n+l1)+RKeR(k,l);
            end
         end
      end
   end
   
% Fim do Loop de montagem da matriz de rigidez global
end
Kgl=Kg;
% Imposição das restrições de apoio na matriz de rigidez global
fprintf('%s\n','Impondo-se as restrições de apoio...');
for n=1:nnos,
	for j=1:3,
      i=3*(n-1)+j;
      if (IC(i)==1)
         F(i)=0.0;
         % Zerando linhas e colunas dos GLs restritos
			Kg(:,i)=0.0;
			Kg(i,:)=0.0;
			Kg(i,i)=1.0;
      end
   end
end

% Resolução da Equação de Equilíbrio:

U(ngl)=0.0;
U=inv(Kg)*F';
      
% Atualizar Geometria

x1(nnos)=0.0;
y1(nnos)=0.0;
for i=1:nnos,
   x1(i)=x(i)+U(3*i-2);
   y1(i)=y(i)+U(3*i-1);
end

%  Plota a estrutura deformada
xd(nnos)=0.0;
yd(nnos)=0.0;
escala=input('Digite a escala a ser utilizada(100):');
if (isempty(escala)==1)
   escala=100;
end
for i=1:nnos,
   xd(i)=x(i)+escala*U(3*i-2);
   yd(i)=y(i)+escala*U(3*i-1);
end
fprintf('%s\n','Desenhando a estrutura deformada...');
figure(2);
clf;
for i=1:nelm,
   linha=line([x(no1(i)) x(no2(i))],[y(no1(i)) y(no2(i))]);
   set(linha,'Color','blue')
   set(linha,'LineStyle','--')
   linha=line([xd(no1(i)) xd(no2(i))],[yd(no1(i)) yd(no2(i))]);
   set(linha,'Color','red')
end
axis equal;
axis off;
title('Estrutura Deformada');

% Obtenção dos esforços
% Def. vetor esforços segundo as coordenadas globais, por elemento
Sg(6)=0.0;
Sg(:)=0.0;
Sgl(nel,6)=0.0;
% Def. vetor esforços por elementos, segundo suas coord. locais 
Sl(nelm,6)=0.0;
Sl(:,:)=0.0;
% Def. vetor reações 
reac(ngl)=0.0;
reac(:)=0.0;
% Início do Loop 
for nel=1:nelm,
   kk(1)=no1(nel);
   kk(2)=no2(nel);
   % Cálculo da matriz de rotação
   dx=x(kk(2))-x(kk(1));
   dy=y(kk(2))-y(kk(1));
   % Cálculo do comprimento do elemento
   L=sqrt(dx^2+dy^2);
   % Cálculo dos cossenos diretores
   cx=dx/L;
   cy=dy/L;
   % Definição da Matriz de Rotação
   R=[ cx  cy  0  0  0  0;
      -cy  cx  0  0  0  0;
      0    0   1  0  0  0
      0    0   0  cx cy 0;
      0    0   0 -cy cx 0;
      0    0   0  0  0  1];
   % Definição da Matriz de Rigidez do elemento de treliça Ke
   E=mod(nel);
   A=area(nel);
   J=mom(nel);
   Kl=[E*A/L  0          0          -E*A/L 0           0;
      0      12*E*J/L^3  6*E*J/L^2  0      -12*E*J/L^3 6*E*J/L^2;
      0      6*E*J/L^2   4*E*J/L    0      -6*E*J/L^2  2*E*J/L;
      -E*A/L 0           0          E*A/L  0           0;
      0      -12*E*J/L^3 -6*E*J/L^2 0      12*E*J/L^3 -6*E*J/L^2;
      0      6*E*J/L^2   2*E*J/L    0      -6*E*J/L^2  4*E*J/L];
   % Cálculo do vetor dos delocamentos locais segundo as coordenadas globais ug
   ug(6)=0.0;
   ug(:)=0.0;
   for m=1:2,
      k=3*(kk(m)-1);
      l=0;
      for n=3*m-2:3*m,
         l=l+1;
         ug(n)=U(k+l);
      end
   end
   % Cálculo do vetor dos delocamentos locais segundo as coordenadas locais ul
   ul=R*ug';
   % Cálculo dos esforços segundo as coordenadas locais ul
   Sl(nel,:)=S0(nel,:)+(Kl*ul)';
   % Cálculo dos esforços segundo as coordenadas globais
   Sg(:)=0.0;
   Sg=(R'*Sl(nel,:)')';
   Sgl(nel,:)=Sg;
   % Cálculo das reações de apoio
   for i=1:2,
      for j=1:3,
         if (IC(3*(kk(i)-1)+j)~=0)
            % Cálculo da coord global m
            m=3*(kk(i)-1)+j;
            n=3*(i-1)+j;
            % Cálculo das reações de apoio
            reac(m)=reac(m)+Sg(n)';
         end
      end
   end
end

% Plota o Diagrama de esforços normais e reações de apoio
fprintf('%s\n','Desenhando o DEN...');
figure(3);
clf;
% Desenha os esforços normais
for i=1:nelm,
   linha=line([x(no1(i)) x(no2(i))],[y(no1(i)) y(no2(i))]);
   set(linha,'Color','blue')
   texto=text((x(no1(i))+x(no2(i)))/2,(y(no1(i))+y(no2(i)))/2,num2str(-Sl(i,1),4));
   set(texto,'Color','red')
end
for i=1:nnr,
   nn=restr(1,i);
   valor=strcat('(',num2str(reac((nn-1)*3+1),3),',',num2str(reac((nn-1)*3+2),3),',',num2str(reac((nn-1)*3+3),3),')');
   texto=text(x(nn),y(nn),valor);
   set(texto,'Color','y')
end
% Desenha as reações de apoio
axis equal;
axis off;
title('Esforços Normais e Reações de Apoio');


fprintf('%s\n','Desenhando o DMF...');
figure(4);
clf;
% Desenha o DMF
mommax=max(abs(Sl(:,3)));
xmax=max(x)/8.0;
prompt=strcat('Digite a escala a ser utilizada(',num2str(xmax),')');
escala=input(prompt);
if (isempty(escala)==0)
   xmax=escala;
end

for nel=1:nelm,
   kk(1)=no1(nel);
   kk(2)=no2(nel);
   % Cálculo da matriz de rotação
   dx=x(kk(2))-x(kk(1));
   dy=y(kk(2))-y(kk(1));
   % Cálculo do comprimento do elemento
   L=sqrt(dx^2+dy^2);
   % Cálculo dos cossenos diretores
   cx=dx/L;
   cy=dy/L;
   X1=x(kk(1));
   X4=x(kk(2));
   Y1=y(kk(1));
   Y4=y(kk(2));
   M1=Sl(nel,3)/mommax*xmax;
   M2=Sl(nel,6)/mommax*xmax;
   X2=X1-M1*cy;
   Y2=Y1+M1*cx;
   X3=X4+M2*cy;
   Y3=Y4-M2*cx;
   linha=line([X1 X4],[Y1 Y4]);
   set(linha,'Color','black')
   aviso1=0;
   for i=1:nvc,
       clear vetorX vetorY vetorX1 vetorY1;
   ne=carreg_vig(1,i);
   if ne==nel
       q=carreg_vig(2,i);
       aviso1=1;
       L2=L*abs(cx);
       p = polyfit([0 L/2 L],[M1 ((M1-M2)/2-q*dx^2/8/mommax*xmax) -M2],2);
       vetorX=0:L/20:L;
       vetorY=polyval(p,vetorX);
       vetorX1=X1*ones(1,length(vetorX))+vetorX*cx;
       vetorX1=vetorX1-vetorY*cy;
       if dy==0,
            vetorY1=Y1*ones(1,21);
       else
            vetorY1=Y1:dy/20:Y4;
       end 
       vetorY1=vetorY1+vetorY*cx;
       grafo=line([X1 vetorX1 X4],[Y1 vetorY1 Y4]);
   end
   end
   if aviso1==0, grafo=line([X1 X2 X3 X4],[Y1 Y2 Y3 Y4]);end
   set(grafo,'Color','blue')
%   texto1=text(X2,Y2,num2str(Sl(nel,3),4));
%   set(texto1,'Color','red')
%   texto2=text(X3,Y3,num2str(Sl(nel,6),4));
%   set(texto2,'Color','red')
end
axis equal;
axis off;
title('Diagrama de Momentos Fletores');
