classdef Barra2D
    % Clase creadas para barras 2D
    %   Detailed explanation goes here

    properties
        Name
        E
        A
        node_j
        node_k
    end

    properties
        Sm
        Sj
        T
    end
    
    properties (Hidden)
        vector % Para calcular longitud y angulo
        nDof=2;
    end

    properties
        L
        alpha
        index
    end

    properties
        Dj
        Aj
        Dm
        Am
    end

    methods
        % Funcion para instanciar las propiedades de las barras (Geometria y rigidez)
        function obj=Barra2D(Name,E,A,node_j,node_k)
            if nargin>0
                obj.Name=Name;
                obj.E=E;
                obj.A=A;
                obj.node_j=node_j;
                obj.node_k=node_k;
                [obj.L,obj.alpha,obj.index]=obj.Geometria;
                [obj.Sm,obj.Sj,obj.T]=obj.Propiedades;
            end
        end
    end

    methods
        % Funciones de apoyo para el calculo de la geometria de la barra
        function [L,alpha,index]=Geometria(obj)
            % Calculamos los datos geometricos de los elementos
            n_j=obj.node_j;
            n_k=obj.node_k;
            obj.vector=n_k.coordenada-n_j.coordenada;
            alpha=atan2d(obj.vector(1,2),obj.vector(1,1));
            L=(obj.vector(1,2)^2+obj.vector(1,1)^2)^0.5;

            % Realizamos un algoritmo para definir los indices de la barra
            i_j=n_j.index;
            i_k=n_k.index;
            index=[i_j,i_k];
        end

        %Funciones de apoyo para el calculo de la rigidez de la barra
        function [Sm,Sj,T] = Propiedades(obj)
            % Calculamos la matriz de rigidez local y global del elemento 
            Sm=obj.E*obj.A/obj.L*[1 -1
                                 -1  1];
            c=cosd(obj.alpha);
            s=sind(obj.alpha);
            Sj=obj.E*obj.A/obj.L*[c^2  c*s   -c^2 -c*s
                                  c*s  s^2   -c*s -s^2
                                 -c^2 -c*s    c^2  c*s
                                 -c*s -s^2    c*s  s^2];
            T=[c s 0 0
               0 0 c s];
        end
 
        function [Dj,Dm,Aj,Am] = Calculos_locales(Dj)
            Dj=Dj(obj.index);
            % Calculo de desplazamientos locales
            Dm=obj.T*obj.Dj;
            % Calculo de fuerzas locales
            Am=obj.Sm*obj.Dm;
            % Calculo de fuerzas de nudos
            Aj=obj.T'*obj.Am;
        end
    end

end
classdef Node2D
    %UNTITLED Summary of this class goes here
    %   numero - Se refiere al codigo individual de cada nudo, los nudos
    %   debes ser enteros, positivos, y tener un valor propio que no se
    %   repita
    %   coordenada - son las coordenadas cartesianas definidas como un vector [i j k]
    %   restrain - corresponde a un vector que indica los grados de
    %   libertad libre y los restringidos [f/r f/r f/r] - f es libres - r
    %   es restringidos
    %   desplazamiento - corresponde a los desplazamientos del nudo

    properties
        node_number
        coordenada
        load
        restrain
        restrain_displacement=[0,0];
        desplazamiento
    end

    properties
        nDof=2;
        index
        boundary_condition;
    end

    methods
        % Funcion para instanciar el objeto
        function obj = Node2D(numero,coordenada,load,restrain)
            % Instanciamos las propiedades del objeto
            %   Detailed explanation goes here
            if nargin>0
                obj.node_number = numero;
                obj.coordenada=coordenada;
                obj.load=load;
                obj.restrain=restrain;
                % Creamos el vector del nudo con los indices
                obj.index=[obj.nDof*obj.node_number-1,obj.nDof*obj.node_number];
                % Creamos el vector del nudo con los indices que relacionan
                % los nudos libres - 1=restringido - 0=libre
                obj.boundary_condition=obj.restrain;
                obj.boundary_condition(obj.boundary_condition~="r")=0;
                obj.boundary_condition(obj.boundary_condition=="r")=1;
                obj.boundary_condition=str2double(obj.boundary_condition);
            end 
        end
    end

end
% Calculo general de celocias
clc
clear variables

% Determinamos datos generales del modelo
E=200*10^9; %[Pa] - Acero
A=(8/100*8/100); %[m^2]
nDof_nudo=2; % Numero de grados de libertad por nudo

% Deteminamos los nudos
% Node2D(numero,coordenada,load,restrain)
nodes=[Node2D(1,[0,0],[0,0],["r","r"])
       Node2D(2,[6,2],[0,-500*10^3],["f","f"])
       Node2D(3,[6,5],[100*10^3,0],["f","f"])
       Node2D(4,[10,2],[0,-500*10^3],["f","f"])
       Node2D(5,[10,5],[0,0],["f","f"])
       Node2D(6,[18,2],[0,0],["r","r"])
       Node2D(7,[18,5],[0,0],["f","f"])];

% Determinamos los elementos
% Barra2D(Name,E,A,node_j,node_k)

elements=[Barra2D("Barra 1",E,A,nodes(1,1),nodes(2,1))
          Barra2D("Barra 2",E,A,nodes(1,1),nodes(3,1))
          Barra2D("Barra 3",E,A,nodes(2,1),nodes(4,1))
          Barra2D("Barra 4",E,A,nodes(4,1),nodes(6,1))
          Barra2D("Barra 5",E,A,nodes(3,1),nodes(5,1))
          Barra2D("Barra 6",E,A,nodes(5,1),nodes(7,1))
          Barra2D("Barra 7",E,A,nodes(6,1),nodes(7,1))
          Barra2D("Barra 8",E,A,nodes(4,1),nodes(5,1))
          Barra2D("Barra 9",E,A,nodes(2,1),nodes(3,1))
          Barra2D("Barra 10",E,A,nodes(3,1),nodes(4,1))
          Barra2D("Barra 11",E,A,nodes(4,1),nodes(7,1))];

factor_escala=50; % Factor de escala para el dibujo
delta=1; % Offset para el dibujo

clear E A

Calculos automaticos
% Determinamos parametros globales del analisis
numero_nudos=height(nodes);
numero_barras=height(elements);
nDof_sistema=numero_nudos*nDof_nudo;

% Ensamblamos la matriz de rigidez global de la estructra
Sj_sistema=zeros(nDof_sistema);

for i=1:numero_barras
    Sj_temp=elements(i,1).Sj;
    index=elements(i,1).index;
    Sj_sistema(index,index)=Sj_sistema(index,index)+Sj_temp;
end

clear i index Sj_temp

% Determinamos los indices requeridos para identificar los grados de
% libertad libres y los grados de libertada restringidos
% Adicionalmente ensamblamos el vector de cargas para la solucion de los
% desplazamientos
% Ensamblamos el vector de desplazamiento de los apoyos para el calculo de
% los desplazamientos del sistema

index_global=zeros(1,nDof_sistema);
index_restrain=zeros(1,nDof_sistema);
A_join_load=zeros(1,nDof_sistema);
D_restrain=zeros(1,nDof_sistema);

for i=1:numero_nudos
    index=nodes(i,1).index;
    boundary_condition=nodes(i,1).boundary_condition;
    index_global(index)=index;
    index_restrain(index)=boundary_condition;
    A_join_load(index)=nodes(i,1).load;
    D_restrain(index)=nodes(i,1).restrain_displacement;
end

A_join_load=A_join_load';
D_restrain=D_restrain';

clear node_temp index boundary_condition i

% Extraemos los indices correspondientes a los nudos libres
[~,idx]=find(index_restrain==0);
index_f=index_global(idx);

% Extraemos los indices correspondeinetes a los nudos restringidos
[~,idx]=find(index_restrain==1);
index_r=index_global(idx);

clear idx index_global

% Particionamos la matriz de rigidez del sistema
% Determinamos la matriz Sff
Sff=Sj_sistema(index_f,index_f);
Sfr=Sj_sistema(index_f,index_r);
Srf=Sj_sistema(index_r,index_f);
Srr=Sj_sistema(index_r,index_r);

% Determinamos el vector de carga de los nudos libres
Af=A_join_load(index_f');

% Determinamos el vector de desplazamientos para los nudos restringidos
Dr=D_restrain(index_r');

% Resolvemos el sistema de ecuaciones para los delsplazamientos
Df=Sff^-1*(Af-Sfr*Dr);

% Ensamblamos el vector de desplazamiento de los nudos
Dj=zeros(nDof_sistema,1);
Dj(index_f')=Df;
Dj(index_r')=D_restrain(index_r');

% Calculamos las acciones de nudos
Aj=Sj_sistema*Dj;

Rutina para almacenar los resultados a las barras
% Grabamos los resultados obtenidos en el objeto
for i=1:numero_barras
    elements(i,1).Dj=Dj(elements(i,1).index);
    % Calculo de desplazamientos locales
    elements(i,1).Dm=elements(i,1).T*elements(i,1).Dj;
    % Calculo de fuerzas locales
    elements(i,1).Am=elements(i,1).Sm*elements(i,1).Dm;
    % Calculo de fuerzas de nudos
    elements(i,1).Aj=elements(i,1).T'*elements(i,1).Am;
end

Rutina para almacenar los desplazamientos de nudos
% Grabamos los deplazamientos de los nudos en el objeto
for i=1:numero_nudos
    nodes(i,1).desplazamiento=Dj(nodes(i,1).index');
end

Rutina para graficar
% Ensamblamos un vector de coordenadas X y Y para dibujar

j=1;
X=zeros(1,numero_barras*2);
X_deformada=zeros(1,numero_barras*2);
Y=zeros(1,numero_barras*2);
Y_deformada=zeros(1,numero_barras*2);

figure
hold on
for i=1:numero_barras
    plot([elements(i,1).node_j.coordenada(1,1),elements(i,1).node_k.coordenada(1,1)],[elements(i,1).node_j.coordenada(1,2),elements(i,1).node_k.coordenada(1,2)],'.r','MarkerSize',15,'Color','k','LineWidth',1,'LineStyle','--')
    plot([elements(i,1).node_j.coordenada(1,1)+elements(i,1).Dj(1,1)*factor_escala,elements(i,1).node_k.coordenada(1,1)+elements(i,1).Dj(3,1)*factor_escala],[elements(i,1).node_j.coordenada(1,2)++elements(i,1).Dj(2,1)*factor_escala,elements(i,1).node_k.coordenada(1,2)+elements(i,1).Dj(4,1)*factor_escala],'.r','MarkerSize',15,'Color','k','LineWidth',1.5,'LineStyle','-')
    % Grabamos las coordenadas en una lista para obtener los maximos y
    % minimos para el dibujo
    X(1,j)=elements(i,1).node_j.coordenada(1,1);
    Y(1,j)=elements(i,1).node_j.coordenada(1,2);
    X_deformada(1,j)=elements(i,1).node_j.coordenada(1,1)+elements(i,1).Dj(1,1)*factor_escala;
    Y_deformada(1,j)=elements(i,1).node_j.coordenada(1,2)+elements(i,1).Dj(2,1)*factor_escala;
    j=j+1;
    X(1,j)=elements(i,1).node_k.coordenada(1,1);
    Y(1,j)=elements(i,1).node_k.coordenada(1,2);
    X_deformada(1,j)=elements(i,1).node_k.coordenada(1,1)+elements(i,1).Dj(3,1)*factor_escala;
    Y_deformada(1,j)=elements(i,1).node_k.coordenada(1,2)+elements(i,1).Dj(4,1)*factor_escala;
    j=j+1;
end

clear j i

grid on
axis([min(X)-delta,max(X)+delta,min(Y)-delta,max(Y)+delta])
title("Estructura")
daspect([1 1 1]) 
set(gcf,'Visible','on','position',[0 0 1920 1080])

clear X Y X_deformada Y_deformada