https://www.youtube.com/watch?v=b9XgSlBBiNw&ab_channel=NicolasMoraBowen

https://www.youtube.com/watch?v=djnzOTUDTqQ

https://www.youtube.com/watch?v=9TJUb1ZFQ4o&ab_channel=NicolasMoraBowen

https://www.youtube.com/watch?v=IDsK5L63b5I&ab_channel=NicolasMoraBowen

https://www.youtube.com/watch?v=rxb0QhewWAA&ab_channel=NicolasMoraBowen

https://www.youtube.com/watch?v=E8Og3AZvC_4&ab_channel=NicolasMoraBowen

https://www.youtube.com/watch?v=KmXWcCqiDnE&ab_channel=NicolasMoraBowen

https://www.youtube.com/watch?v=nYYqGoKWrFQ&ab_channel=NicolasMoraBowen

classdef nudo < Preferencias
    %UNTITLED2 Summary of this class goes here
    %   Detailed explanation goes here

    properties
        numero
        coordenadas
        condiciones_borde
        load
        index
        Dj_existentes=[0 0 0];
        Dj
    end

    methods
        % Funcion para instanciar el objeto
        function obj = nudo(numero,coordenadas,condiciones_borde,load)
            if nargin>0
                obj.numero=numero;
                obj.coordenadas=coordenadas;
                obj.condiciones_borde=condiciones_borde;
                obj.load=load;
                obj.index=obj.indices();
            end
        end
        function [idx]=indices(obj)
            idx=zeros(1,obj.nDof);
            for i=1:obj.nDof
                j=obj.nDof-i;
                idx(1,i)=(obj.nDof*obj.numero)-j;
            end
        end
    end

end
classdef barra < Preferencias
    %UNTITLED5 Summary of this class goes here
    %   Detailed explanation goes here

    properties
        nudo_j
        nudo_k
        Material
        Seccion
        L
        alpha
        sm
        T
        Sm
        index
        Dm
        dm
        Am_EQ
        am_EQ
        x_local
        a_internas_EQ
        A_internas_EQ
        T_coordenadas
        am_FE=[0;0;0;0;0;0];
        Am_FE=[0;0;0;0;0;0];
        a_internas_FE
        A_internas_FE
        a_internas
        A_internas
        am
        Am
        % a para articulado, na para no articulado
        releases=["na","na"]
        check_releases
    end

    properties (Hidden)
        numero_divisiones=10;
        sm_base
    end

    methods (Static)
        function phi=Calculo_Matriz_de_Equilibrio(x)
            phi=[-1  0  0
                  0 -1  0
                  0  x -1];
        end
    end

    methods
        function obj = barra(nudo_j,nudo_k,Material,Seccion,releases)
            %BARRA Instancia el objeto utilizando
            %   barra(nudo_j,nudo_k,Material,Seccion,releases)
            if nargin>0
                % Asignamos las propiedades del objeto
                obj.nudo_j = nudo_j;
                obj.nudo_k = nudo_k;
                obj.Material = Material;
                obj.Seccion = Seccion;
                obj.releases=releases;
                
                % Calculamos las propiedades del objeto
                obj.check_releases=obj.Validacion_Articulaciones();
                [obj.L,obj.alpha]=obj.Calculo_Geometria();
                obj.index=obj.Calculo_indices();
                [obj.x_local,obj.a_internas_FE]=obj.Calculo_Subdivision_X();
                % Calculo de matrices de transformacion
                obj.T=obj.Calculo_Matriz_Transformacion();
                obj.T_coordenadas=obj.Calculo_Matriz_Transformacion_Coordenadas();
                % Calculos de matrices de rigidez
                obj.sm_base=obj.Calculo_Matriz_Rigidez_Local();
                obj.sm=obj.Calculo_Matriz_Rigidez_Local_Articulada();
                obj.Sm=obj.Calculo_Matriz_Rigidez_Global();
                   
            end
        end

        function obj=Asignar_Carga_Puntual(obj,Carga,a,angulo)
            % Asignar_Carga_Puntual Asigna las cargas a_FE, A_FE, y
            % a_internas_FE
            % Cargas
            % a
            % alpha
            [obj.am_FE,obj.Am_FE,obj.a_internas_FE]=obj.Carga_Puntual(Carga,a,angulo);
        end
        function obj=Asignar_Momento_Puntual(obj,Carga,a)
            % Asignar_Carga_Puntual Asigna las cargas a_FE, A_FE, y
            % a_internas_FE
            % Cargas
            % a
            % alpha
            [obj.am_FE,obj.Am_FE,obj.a_internas_FE]=obj.Momento_Puntual(Carga,a);
        end
        function obj=Asignar_Carga_Distribuida(obj,wj,wk,angulo)
            % Asignar_Carga_Puntual Asigna las cargas a_FE, A_FE, y
            % a_internas_FE
            % Cargas
            % a
            % alpha
            [obj.am_FE,obj.Am_FE,obj.a_internas_FE]=obj.Carga_Distribuida(wj,wk,angulo);
        end

        % Funciones para incluir cargas entre nudos
        function [a_FE,A_FE,a_internas_FE]=Carga_Puntual(obj,A,a,angulo)
            % A = Carga →+ ↑+ 
            % angulo = angulo de la carga respecto a la barra 
            % a = Distancia del extremo j al punto de aplicacion de la carga 
            b=obj.L-a;
            Ax=A*cosd(angulo);
            Ay=A*sind(angulo);
            u1=(12*obj.Material.E*obj.Seccion.Inercia)/(obj.Material.G*obj.Seccion.Area_cortante*a^2);
            u2=(12*obj.Material.E*obj.Seccion.Inercia)/(obj.Material.G*obj.Seccion.Area_cortante*b^2);
            if a<obj.L
                a_FE_new=[-Ax*(b/(a+b))
                         -(Ay*(u2+1)*b^3+3*Ay*a*b^2)/(3*a*b^2+3*a^2*b+a^3*u1+a^3+b^3*(u2+1))
                          ((Ay*b^2)/(2*(a+b)))-((Ay*b^2*(a+b)*(3*a+b+b*u2))/(2*(3*a*b^2+3*a^2*b+b^3*u2+b^3+a^3*(u1+1))))
                          -Ax*(a/(a+b))
                          -(Ay*(u1+1)*a^3+3*Ay*b*a^2)/(3*a*b^2+3*a^2*b+a^3*u1+a^3+b^3*(u2+1))
                          ((Ay*a^2*(a+b)*(a+3*b+a*u1))/(2*(3*a*b^2+3*a^2*b+a^3*u1+a^3+b^3*(u2+1))))-((Ay*a^2)/(2*(a+b)))];
                % Sin Deformaciones por cortante
%                 a_FE_new=[-A*cosd(angulo)*(1-a/obj.L)
%                          (-A*sind(angulo)*(obj.L+2*a)*(b)^2)/obj.L^3
%                          (-A*sind(angulo)*a*b^2)/obj.L^2
%                           -A*cosd(angulo)*(a/obj.L)
%                          (-A*sind(angulo)*(3*obj.L-2*a)*(a)^2)/obj.L^3
%                          ( A*sind(angulo)*a^2*b)/obj.L^2];
                if obj.check_releases==true
                    index_articulacion=find(obj.releases=="a");
                    index_0=index_articulacion*obj.nDof;
                    index_c=1:1:obj.nDof*2;
                    index_c(index_0)=[];
                    
                    sff=obj.sm_base(index_0,index_0);
                    srf=obj.sm_base(index_c,index_0);
                    af_FE=a_FE_new(index_0);
                    ar_FE=a_FE_new(index_c);
                    
                    df=-sff^-1*af_FE;
                    ar=srf*df+ar_FE;
                    
                    a_FE_new=zeros(1,obj.nDof*2)';
                    a_FE_new(index_c)=ar;
                end
                a_FE=obj.am_FE+a_FE_new;
                A_FE=obj.T'*a_FE;
            end

            % Calculamos las acciones internas por el empotramiento

            a_internas_FE_new=zeros(3,width(obj.x_local));
            a_nudo_j=a_FE_new(1:3,1);
                for i=1:width(obj.x_local)
                    phi_x=obj.Calculo_Matriz_de_Equilibrio(obj.x_local(i));
                    if obj.x_local(i)<a
                        An=[0;0;0];
                    else
                        An=[-A*cosd(angulo)
                            -A*sind(angulo)
                             A*sind(angulo)*(obj.x_local(i)-a)];
                    end
                    a_internas_FE_new(:,i)=phi_x*a_nudo_j+An;
                end
                % Grabamos sobre el vector original para poder usar
                % multiples cargas entre nudos
                a_internas_FE=obj.a_internas_FE;
                a_internas_FE.ai(2,2:end-1)=a_internas_FE.ai(2,2:end-1)+a_internas_FE_new(1,:);
                a_internas_FE.vi(2,2:end-1)=a_internas_FE.vi(2,2:end-1)+a_internas_FE_new(2,:);
                a_internas_FE.mi(2,2:end-1)=a_internas_FE.mi(2,2:end-1)+a_internas_FE_new(3,:);
        end
        function [a_FE,A_FE,a_internas_FE]=Momento_Puntual(obj,A,a)
            % A = Carga + antihorario
            % a = Distancia del extremo j al punto de aplicacion de la carga 
            b=obj.L-a;
            u1=(12*obj.Material.E*obj.Seccion.Inercia)/(obj.Material.G*obj.Seccion.Area_cortante*a^2);
            u2=(12*obj.Material.E*obj.Seccion.Inercia)/(obj.Material.G*obj.Seccion.Area_cortante*b^2);
            if a<obj.L
                a_FE_new=[0
                          (6*A*a*b)/(3*a*b^2+3*a^2*b+a^3*u1+a^3+b^3*(u2+1))
                          (3*A*a*b*(a+b))/(3*a*b^2+3*a^2*b+a^3*u1+a^3+b^3*(u2+1))-(A*b)/(a+b)
                          0
                          -(6*A*a*b)/(3*a*b^2+3*a^2*b+a^3*u1+a^3+b^3*(u2+1))
                          (3*A*a*b*(a+b))/(3*a*b^2+3*a^2*b+a^3*u1+a^3+b^3*(u2+1))-(A*a)/(a+b)];
                % Sin deformaciones por cortante
%                 a_FE_new=[0
%                          (6*A*a*b)/obj.L^3
%                          (A*b*(2*a-b))/obj.L^2
%                           0
%                          (-6*A*a*b)/obj.L^3
%                          (A*a*(2*b-a))/obj.L^2];
                if obj.check_releases==true
                    index_articulacion=find(obj.releases=="a");
                    index_0=index_articulacion*obj.nDof;
                    index_c=1:1:obj.nDof*2;
                    index_c(index_0)=[];
                    
                    sff=obj.sm_base(index_0,index_0);
                    srf=obj.sm_base(index_c,index_0);
                    af_FE=a_FE_new(index_0);
                    ar_FE=a_FE_new(index_c);
                    
                    df=-sff^-1*af_FE;
                    ar=srf*df+ar_FE;
                    
                    a_FE_new=zeros(1,obj.nDof*2)';
                    a_FE_new(index_c)=ar;
                end                

                a_FE=obj.am_FE+a_FE_new;
                A_FE=obj.T'*a_FE;
            end

            % Calculamos las acciones internas por el empotramiento

            a_internas_FE_new=zeros(3,width(obj.x_local));
            a_nudo_j=a_FE_new(1:3,1);
                for i=1:width(obj.x_local)
                    phi_x=obj.Calculo_Matriz_de_Equilibrio(obj.x_local(i));
                    if obj.x_local(i)<a
                        An=[0;0;0];
                    else
                        An=[0
                            0
                           -A];
                    end
                    a_internas_FE_new(:,i)=phi_x*a_nudo_j+An;
                end
                % Grabamos sobre el vector original para poder usar
                % multiples cargas entre nudos
                a_internas_FE=obj.a_internas_FE;
                a_internas_FE.ai(2,2:end-1)=a_internas_FE.ai(2,2:end-1)+a_internas_FE_new(1,:);
                a_internas_FE.vi(2,2:end-1)=a_internas_FE.vi(2,2:end-1)+a_internas_FE_new(2,:);
                a_internas_FE.mi(2,2:end-1)=a_internas_FE.mi(2,2:end-1)+a_internas_FE_new(3,:); 

        end
        function [a_FE,A_FE,a_internas_FE]=Carga_Distribuida(obj,wj,wk,angulo)
            % wj = Carga →+ ↑+ 
            % wk = Carga →+ ↑+
            % angulo = angulo de la carga respecto a la barra
            w2=wj-wk; % Porcion de la carga triangular
            w1=wk; % Porcion rectangular
            a_FE_new=[-w1*cosd(angulo)*obj.L/2-w2*cosd(angulo)*obj.L/3
                      -w1*sind(angulo)*obj.L/2-7*w2*sind(angulo)*obj.L/20
                      -w1*sind(angulo)*obj.L^2/12-w2*sind(angulo)*obj.L^2/20
                      -w1*cosd(angulo)*obj.L/2-w2*cosd(angulo)*obj.L/6
                      -w1*sind(angulo)*obj.L/2-3*w2*sind(angulo)*obj.L/20
                       w1*sind(angulo)*obj.L^2/12+w2*sind(angulo)*obj.L^2/30];

            % Calculamos la correccion por articulaciones
   
            if obj.check_releases==true
                index_articulacion=find(obj.releases=="a");
                index_0=index_articulacion*obj.nDof;
                index_c=1:1:obj.nDof*2;
                index_c(index_0)=[];
                
                sff=obj.sm_base(index_0,index_0);
                srf=obj.sm_base(index_c,index_0);
                af_FE=a_FE_new(index_0);
                ar_FE=a_FE_new(index_c);
                
                df=-sff^-1*af_FE;
                ar=srf*df+ar_FE;
                
                a_FE_new=zeros(1,obj.nDof*2)';
                a_FE_new(index_c)=ar;
            end

            a_FE=obj.am_FE+a_FE_new;
            A_FE=obj.T'*a_FE;

            % Calculamos las acciones internas por el empotramiento

            a_internas_FE_new=zeros(3,width(obj.x_local));
            a_nudo_j=a_FE_new(1:3,1);
                for i=1:width(obj.x_local)
                    phi_x=obj.Calculo_Matriz_de_Equilibrio(obj.x_local(i));
                    w_x =((wk-wj)/obj.L)*obj.x_local(i)+wj;
                    An=[-cosd(angulo)*((wj+w_x)/2)*obj.x_local(i)
                        -sind(angulo)*((wj+w_x)/2)*obj.x_local(i)
                         sind(angulo)*((2*wj+w_x)/6)*obj.x_local(i)^2];
                    a_internas_FE_new(:,i)=phi_x*a_nudo_j+An;
                end
                % Grabamos sobre el vector original para poder usar
                % multiples cargas entre nudos
                a_internas_FE=obj.a_internas_FE;
                a_internas_FE.ai(2,2:end-1)=a_internas_FE.ai(2,2:end-1)+a_internas_FE_new(1,:);
                a_internas_FE.vi(2,2:end-1)=a_internas_FE.vi(2,2:end-1)+a_internas_FE_new(2,:);
                a_internas_FE.mi(2,2:end-1)=a_internas_FE.mi(2,2:end-1)+a_internas_FE_new(3,:);

        end
        function Graficos_FE(obj)
          
            figure
            hold on
            
            subplot(3,1,1)
            pgon=polyshape(obj.a_internas_FE.mi(1,:),-obj.a_internas_FE.mi(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama de Momentos: Cordenadas Locales')
            xlabel('Distancia')
            ylabel('Momentos: Locales')
            grid on
            subplot(3,1,2)
            pgon=polyshape(obj.a_internas_FE.vi(1,:),obj.a_internas_FE.vi(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama de Cortante: Cordenadas Locales')
            xlabel('Distancia')
            ylabel('Cortantes: Locales')
            grid on
            subplot(3,1,3)
            pgon=polyshape(obj.a_internas_FE.ai(1,:),obj.a_internas_FE.ai(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama Axial: Cordenadas Locales')
            xlabel('Distancia')
            ylabel('Axiales: Locales')
            grid on

            set(gcf,'Visible','on','position',[0 0 1920 1080])
        end
        
        % Funciones propiedades del objeto
        function Check=Validacion_Articulaciones (obj)
            temp=find(obj.releases=="a",1);
            if isempty(temp)
                Check=false;
            else
                Check=true;
            end
        end
        function [L,alpha]=Calculo_Geometria(obj)
            Delta=obj.nudo_k.coordenadas-obj.nudo_j.coordenadas;
            L=((Delta(1,1))^2+(Delta(1,2))^2)^0.5;
            alpha=atan2d(Delta(1,2),Delta(1,1));
        end
        function [idx]=Calculo_indices(obj)
            idx=[obj.nudo_j.index,obj.nudo_k.index];
        end
        function [x_locales,a_internas_FE]=Calculo_Subdivision_X(obj)
             x_locales=0:obj.L/obj.numero_divisiones:obj.L;
             % Creamos un vector de ceros para poder grabar varias cargas entre nudo
             a_internas_FE.ai=zeros(2,width(x_locales)+2);
             a_internas_FE.vi=zeros(2,width(x_locales)+2);
             a_internas_FE.mi=zeros(2,width(x_locales)+2);
             a_internas_FE.ai(1,2:end)=[x_locales,obj.L];
             a_internas_FE.vi(1,2:end)=[x_locales,obj.L];
             a_internas_FE.mi(1,2:end)=[x_locales,obj.L];                  
        end

        % Funciones de matrices de transformacion
        function [T]=Calculo_Matriz_Transformacion(obj)

            c=cosd(obj.alpha);
            s=sind(obj.alpha);

            T=[c  s  0  0  0  0
                -s  c  0  0  0  0
                0  0  1  0  0  0
                0  0  0  c  s  0
                0  0  0 -s  c  0
                0  0  0  0  0  1];

        end
        function T_coordenadas=Calculo_Matriz_Transformacion_Coordenadas(obj)
            c=cosd(obj.alpha);
            s=sind(obj.alpha);
            T_coordenadas=[c -s
                           s  c];
        end

        % Funciones para el calculo de matrices de rigidez
        function [sm]=Calculo_Matriz_Rigidez_Local(obj)

            u=(12*obj.Material.E*obj.Seccion.Inercia)/(obj.Material.G*obj.L^2*obj.Seccion.Area_cortante);
            %u=0;

            a=obj.Seccion.Area*obj.Material.E/obj.L;
            b=12*obj.Material.E*obj.Seccion.Inercia/((1+u)*obj.L^3);
            c=6*obj.Material.E*obj.Seccion.Inercia/((1+u)*obj.L^2);
            d=(4+u)*(obj.Material.E*obj.Seccion.Inercia)/((1+u)*obj.L);
            e=(2-u)*(obj.Material.E*obj.Seccion.Inercia)/((1+u)*obj.L);

            sm=[ a  0  0 -a  0  0
                 0  b  c  0 -b  c
                 0  c  d  0 -c  e
                -a  0  0  a  0  0
                 0 -b -c  0  b -c
                 0  c  e  0 -c  d];

        end
        function [Sm]=Calculo_Matriz_Rigidez_Global(obj)
            Sm=obj.T'*obj.sm*obj.T;
        end
        function [sm_mod]=Calculo_Matriz_Rigidez_Local_Articulada(obj)
            
            if obj.check_releases==false
                sm_mod=obj.sm_base;
            else
            index_articulacion=find(obj.releases=="a");
            index_0=index_articulacion*obj.nDof;
            index_c=1:1:obj.nDof*2;
            index_c(index_0)=[];
            
            scc=obj.sm_base(index_c,index_c);
            sc0=obj.sm_base(index_c,index_0);
            s00=obj.sm_base(index_0,index_0);
            s0c=obj.sm_base(index_0,index_c);
            
            sm_calc=scc-sc0*s00^-1*s0c;
            sm_mod=zeros(2*obj.nDof,2*obj.nDof);
            sm_mod(index_c,index_c)=sm_calc;
            end

        end
        
        
        function Matriz_Transformada=Transformacion_Local_a_Global(obj,local)
            % Para un objeto tomamos un struct que contienen las matrices
            % locales, se lo escala, rota, y desplaza; y devolvemos un
            % nuevo struct
            % Determinamos las matrices de accione internas con los factores de escala
            a_escalado=(local.ai).*[1;obj.factor_escala_axial];
            v_escalado=(local.vi).*[1;obj.factor_escala_corte];
            m_escalado=(local.mi).*[1;obj.factor_escala_momento];
            % Generar las matrices de acciones rotadas
            a_rot=obj.T_coordenadas*(a_escalado);
            v_rot=obj.T_coordenadas*(v_escalado);
            m_rot=obj.T_coordenadas*(m_escalado);
            % Desplazamos el objeto en el espacio
            Matriz_Transformada.Ai=a_rot+[obj.nudo_j.coordenadas(1,1);obj.nudo_j.coordenadas(1,2)];
            Matriz_Transformada.Vi=v_rot+[obj.nudo_j.coordenadas(1,1);obj.nudo_j.coordenadas(1,2)];
            Matriz_Transformada.Mi=m_rot+[obj.nudo_j.coordenadas(1,1);obj.nudo_j.coordenadas(1,2)];
        end
        function Diagramas_Acciones_Internas(obj)
          
            figure
            hold on
            
            subplot(3,3,1)
            pgon=polyshape(obj.a_internas.mi(1,:),-obj.a_internas.mi(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama de Momentos')
            xlabel('Distancia')
            ylabel('Momentos: Locales')
            grid on
            subplot(3,3,2)
            pgon=polyshape(obj.a_internas_EQ.mi(1,:),-obj.a_internas_EQ.mi(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama de Momento: Equivalente')
            xlabel('Distancia')
            ylabel('Momento: Locales')
            grid on
            subplot(3,3,3)
            pgon=polyshape(obj.a_internas_FE.mi(1,:),-obj.a_internas_FE.mi(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama Momento: Empotrado')
            xlabel('Distancia')
            ylabel('Momento: Locales')
            grid on
            subplot(3,3,4)
            pgon=polyshape(obj.a_internas.vi(1,:),obj.a_internas.vi(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama de Cortante')
            xlabel('Distancia')
            ylabel('Momentos: Locales')
            grid on
            subplot(3,3,5)
            pgon=polyshape(obj.a_internas_EQ.vi(1,:),obj.a_internas_EQ.vi(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama de Cortante: Equivalente')
            xlabel('Distancia')
            ylabel('Cortante: Locales')
            grid on
            subplot(3,3,6)
            pgon=polyshape(obj.a_internas_FE.vi(1,:),obj.a_internas_FE.vi(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama Cortante: Empotrado')
            xlabel('Distancia')
            ylabel('Cortante: Locales')
            grid on
            subplot(3,3,7)
            pgon=polyshape(obj.a_internas.ai(1,:),obj.a_internas.ai(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama de Axial')
            xlabel('Distancia')
            ylabel('Axial: Locales')
            grid on
            subplot(3,3,8)
            pgon=polyshape(obj.a_internas_EQ.ai(1,:),obj.a_internas_EQ.ai(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama de Axial: Equivalente')
            xlabel('Distancia')
            ylabel('Axial: Locales')
            grid on
            subplot(3,3,9)
            pgon=polyshape(obj.a_internas_FE.ai(1,:),obj.a_internas_FE.ai(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama Axial: Empotrado')
            xlabel('Distancia')
            ylabel('Axiales: Locales')
            grid on

            set(gcf,'Visible','on','position',[0 0 1920 1080])
        end

        % Funciones utilizadas posterior al calculo y almacenamiento de Df
        
        
        % Funciones utilizadas para solucion
        function obj=Calculo_Desplazamientos_Locales(obj)
            obj.dm=obj.T*obj.Dm;
        end
        function obj=Calculo_Acciones(obj)
            obj.Am_EQ=obj.Sm*obj.Dm;
            obj.am_EQ=obj.sm*obj.dm;
            obj.am=obj.am_EQ+obj.am_FE;
            obj.Am=obj.Am_EQ+obj.Am_FE;
        end
        function obj=Calculo_Acciones_Internas(obj)
            % Separamos las fuerzas del nudo j
            a_nudo_j=obj.am_EQ(1:3,1);
            resumen=zeros(3,width(obj.x_local));
                for i=1:width(obj.x_local)
                    phi_x=obj.Calculo_Matriz_de_Equilibrio(obj.x_local(i));
                    resumen(:,i)=phi_x*a_nudo_j;
                end
             obj.a_internas_EQ.ai=[0,obj.x_local,obj.L;0,resumen(1,:),0];
             obj.a_internas_EQ.vi=[0,obj.x_local,obj.L;0,resumen(2,:),0];
             obj.a_internas_EQ.mi=[0,obj.x_local,obj.L;0,resumen(3,:),0];
             
             obj.a_internas.ai=obj.a_internas_EQ.ai+[zeros(1,width(obj.a_internas_EQ.ai));obj.a_internas_FE.ai(2,:)];
             obj.a_internas.vi=obj.a_internas_EQ.vi+[zeros(1,width(obj.a_internas_EQ.vi));obj.a_internas_FE.vi(2,:)];
             obj.a_internas.mi=obj.a_internas_EQ.mi+[zeros(1,width(obj.a_internas_EQ.mi));obj.a_internas_FE.mi(2,:)];
        end
        function obj=Calculo_Acciones_Internas_Globales(obj)
            obj.A_internas_EQ=obj.Transformacion_Local_a_Global(obj.a_internas_EQ);
            obj.A_internas_FE=obj.Transformacion_Local_a_Global(obj.a_internas_FE);
            obj.A_internas=obj.Transformacion_Local_a_Global(obj.a_internas);
        end
        
        function Graficos_Fuerzas_Internas(obj)
            figure
            hold on
            
            subplot(3,1,1)
            pgon=polyshape(obj.a_internas_EQ.mi(1,:),-obj.a_internas_EQ.mi(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama de Momentos: Cordenadas Locales')
            xlabel('Distancia')
            ylabel('Momentos: Locales')
            grid on
            subplot(3,1,2)
            pgon=polyshape(obj.a_internas_EQ.vi(1,:),obj.a_internas_EQ.vi(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama de Cortante: Cordenadas Locales')
            xlabel('Distancia')
            ylabel('Cortantes: Locales')
            grid on
            subplot(3,1,3)
            pgon=polyshape(obj.a_internas_EQ.ai(1,:),obj.a_internas_EQ.ai(2,:),'Simplify',false);
            plot(pgon)
            title('Diagrama Axial: Cordenadas Locales')
            xlabel('Distancia')
            ylabel('Axiales: Locales')
            grid on

            set(gcf,'Visible','on','position',[0 0 1920 1080])
        end
        

    end

end
classdef Material
    %MATERIAL Clase que contienen las propiedades del material
    %   E = Modulos de elasticidad

    properties
        E
        G
        u
    end

    methods
        function obj = Material(E,u)
            % Funcion para instanciar la clase con los datos del material
            %   obj = Material(E)
            if nargin>0
                obj.E = E;
                obj.u=u;
                obj.G=obj.E/(2*(1+u));
            end
        end

    end
end
classdef Seccion
    %Clase de la Seccion
    %   Detailed explanation goes here

    properties
        Area
        Inercia
        Area_cortante
    end

    methods
        function obj = Seccion(A,I,Av)
            % Funcion que instancia el obejeto de la seccion
            %   Detailed explanation goes here
            if nargin>0
                obj.Area = A;
                obj.Inercia = I;
                obj.Area_cortante=Av;
            end
        end

    end
end
classdef Preferencias
    %UNTITLED6 Summary of this class goes here
    %   Detailed explanation goes here

    properties
        nDof=3;
        factor_escala_momento=-1/2000;
        factor_escala_corte=1/10;
        factor_escala_axial=1/10;
    end

end
classdef Modelo < Preferencias
    %UNTITLED4 Summary of this class goes here
    %   Detailed explanation goes here

    properties
        Nudos
        Barras
        Dj
        Aj
        Sj
    end

    methods
        function obj = Modelo(Nudos,Barras)
            if nargin>0
                obj.Nudos=Nudos;
                obj.Barras=Barras;
            end

        end

        function obj=Calculos(obj)
            % Calcular la matriz de rigidez del sistema
                nDof=obj.nDof;
                nudos=obj.Nudos;
                barras=obj.Barras;

                numero_nudos=height(nudos);
                
                Sj=zeros(nDof*numero_nudos);
                fuerzas_FE=zeros(1,nDof*numero_nudos);
                
                for i=1:height(barras)
                    % Ensamblamos la matriz de rigidez global del sistema
                    Sj(barras(i,1).index,barras(i,1).index)=Sj(barras(i,1).index,barras(i,1).index)+barras(i,1).Sm;
                    % Ensamblamos el vector de fuerzas de empotramiento perfecto del
                    % sistema
                    fuerzas_FE(barras(i,1).index)=fuerzas_FE(barras(i,1).index)+barras(i,1).Am_FE';
                end
                
                
                indices_sistema=zeros(1,nDof*numero_nudos);
                desplazamientos_iniciales=zeros(1,nDof*numero_nudos);
                fuerzas_nudos=zeros(1,nDof*numero_nudos);
                restricciones_sistema=repmat("",1,nDof*numero_nudos);
                
                % Recopilamos la informacion del sistema
                for i=1:height(nudos)
                    % Determinamos el vector con los indices del sistema
                    indices_sistema(nudos(i,1).index)=nudos(i,1).index;
                    % Determinamos un vector con las condiciones de borde del sistema
                    restricciones_sistema(nudos(i,1).index)=nudos(i,1).condiciones_borde;
                    % Determinamos el vector de desplazamientos iniciales del sistema
                    desplazamientos_iniciales(nudos(i,1).index)=nudos(i,1).Dj_existentes;
                    % Determinamos el vector de fuerzas aplicadas a los nudos
                    fuerzas_nudos(nudos(i,1).index)=nudos(i,1).load;
                end
                
                % Determinamos un vector con los indices en funcion de las condiciones de
                % borde
                index_free=find(restricciones_sistema=="f");
                index_restringido=find(restricciones_sistema=="r");
                
                % Determinar las matrices Sff, Sfr, Srf, Srr
                Sff=Sj(index_free,index_free);
                Sfr=Sj(index_free,index_restringido);
                Srf=Sj(index_restringido,index_free);
                Srr=Sj(index_restringido,index_restringido);
                
                % Resolver [Df]=[Sff]^-1([Af]-[Af_FE]-[Sfr][Dr])
                Dr=desplazamientos_iniciales(index_restringido)';
                Af_nudo=fuerzas_nudos(index_free)';
                Ar_nudo=fuerzas_nudos(index_restringido)';
                Af_FE=fuerzas_FE(index_free)';
                Ar_FE=fuerzas_FE(index_restringido)';
                
                Df=Sff^-1*(Af_nudo-Af_FE-Sfr*Dr);
                
                % Calculo de Desplazamientos del Sistema
                
                Dj=zeros(nDof*numero_nudos,1);
                Dj(index_free)=Df;
                Dj(index_restringido)=Dr;
                
                % Calculo de las reacciones
                % Resolver [Ar]+[Ar_nudo]-[Ar_FE]=[Srf][Df]+[Srr][Dr]
                
                Ar=Srf*Df+Srr*Dr-Ar_nudo+Ar_FE;
                
                % Ensamblamos el vector de fuerzas del nudo del sistema
                
                Aj=zeros(nDof*numero_nudos,1);
                Aj(index_free)=Af_nudo;
                Aj(index_restringido)=Ar;
                
                obj.Dj=Dj;
                obj.Aj=Aj;
                obj.Sj=Sj;
        end
        function obj=Resultados_Nudos(obj)
            nudos=obj.Nudos;
            for i=1:height(obj.Nudos)
                nudos(i,1).Dj=obj.Dj(nudos(i,1).index);
            end
            obj.Nudos=nudos;
        end
        function obj=Resultados_Barras(obj)
            barras=obj.Barras;
            for i=1:height(barras)
                % Grabamos los desplazamientos calculados en la barra
                barras(i,1).Dm=obj.Dj(barras(i,1).index);
                % Corremos la funcion interna para el calculo en la barra
                barras(i,1)=barras(i,1).Calculo_Desplazamientos_Locales();
                barras(i,1)=barras(i,1).Calculo_Acciones();
                barras(i,1)=barras(i,1).Calculo_Acciones_Internas();
                barras(i,1)=barras(i,1).Calculo_Acciones_Internas_Globales();
            end
            obj.Barras=barras;
        end
        function Esquema_Geometrico(obj)
            nudos=obj.Nudos;
            barras=obj.Barras;
            figure
            daspect([1 1 1]);
            hold on
            daspect([1 1 1]);
            for i=1:height(nudos)
                plot(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),'.','MarkerSize',15,'Color','k')
                text(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),num2str(nudos(i,1).numero))
            end
            
            for i=1:height(barras)
                XY=[barras(i,1).nudo_j.coordenadas(1,1) barras(i,1).nudo_k.coordenadas(1,1)
                    barras(i,1).nudo_j.coordenadas(1,2) barras(i,1).nudo_k.coordenadas(1,2)];
                line(XY(1,:),XY(2,:),'Color','k','LineWidth',0.75)
            end
            
            title('Esquema Geometrico')
            
            grid on
        end
        function Esquema_Nudos(obj)
            nudos=obj.Nudos;
            figure
            daspect([1 1 1]);
            hold on
            daspect([1 1 1]);
            for i=1:height(nudos)
                plot(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),'.','MarkerSize',15,'Color','k')
                text(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),num2str(nudos(i,1).numero))
            end
            
            
            title('Esquema Geometrico')
            
            grid on
        end
        function Diagramas(obj)
            nudos=obj.Nudos;
            barras=obj.Barras;
            
            figure

            subplot(4,3,2)
            daspect([1 1 1]);
            hold on
            daspect([1 1 1]);
            
            for i=1:height(nudos)
                plot(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),'.','MarkerSize',15,'Color','k')
            end
            
            for i=1:height(barras)
                XY=[barras(i,1).nudo_j.coordenadas(1,1) barras(i,1).nudo_k.coordenadas(1,1)
                    barras(i,1).nudo_j.coordenadas(1,2) barras(i,1).nudo_k.coordenadas(1,2)];
                line(XY(1,:),XY(2,:),'Color','k','LineWidth',0.75)
            end
            
            title('Esquema Geometrico')
            
            grid on
            
            subplot(4,3,5)
            daspect([1 1 1]);
            hold on
            
            for i=1:height(nudos)
                plot(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),'.','MarkerSize',15,'Color','k')
            end
            
            for i=1:height(barras)
                XY=[barras(i,1).nudo_j.coordenadas(1,1) barras(i,1).nudo_k.coordenadas(1,1)
                    barras(i,1).nudo_j.coordenadas(1,2) barras(i,1).nudo_k.coordenadas(1,2)];
                line(XY(1,:),XY(2,:),'Color','k','LineWidth',0.75)    
                pgon=polyshape(barras(i,1).A_internas_EQ.Mi(1,:),barras(i,1).A_internas_EQ.Mi(2,:),'Simplify',false);
                plot(pgon)
            end
            
            grid on
            title('Diagrama de Momentos Equivalentes')
            
            subplot(4,3,8)
            daspect([1 1 1]);
            hold on
            
            for i=1:height(nudos)
                plot(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),'.','MarkerSize',15,'Color','k')
            end
            
            for i=1:height(barras)
                XY=[barras(i,1).nudo_j.coordenadas(1,1) barras(i,1).nudo_k.coordenadas(1,1)
                    barras(i,1).nudo_j.coordenadas(1,2) barras(i,1).nudo_k.coordenadas(1,2)];
                line(XY(1,:),XY(2,:),'Color','k','LineWidth',0.75)    
                pgon=polyshape(barras(i,1).A_internas_EQ.Vi(1,:),barras(i,1).A_internas_EQ.Vi(2,:),'Simplify',false);
                plot(pgon)
            end
            
            grid on
            title('Diagrama de Cortantes Equivalentes')
            
            subplot(4,3,11)
            daspect([1 1 1]);
            hold on
            
            for i=1:height(nudos)
                plot(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),'.','MarkerSize',15,'Color','k')
            end
            
            for i=1:height(barras)
                XY=[barras(i,1).nudo_j.coordenadas(1,1) barras(i,1).nudo_k.coordenadas(1,1)
                    barras(i,1).nudo_j.coordenadas(1,2) barras(i,1).nudo_k.coordenadas(1,2)];
                line(XY(1,:),XY(2,:),'Color','k','LineWidth',0.75)    
                pgon=polyshape(barras(i,1).A_internas_EQ.Ai(1,:),barras(i,1).A_internas_EQ.Ai(2,:),'Simplify',false);
                plot(pgon)
            end
            
            grid on
            title('Diagrama de Axiales Equivalentes')
            
            subplot(4,3,6)
            daspect([1 1 1]);
            hold on
            
            for i=1:height(nudos)
                plot(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),'.','MarkerSize',15,'Color','k')
            end
            
            for i=1:height(barras)
                XY=[barras(i,1).nudo_j.coordenadas(1,1) barras(i,1).nudo_k.coordenadas(1,1)
                    barras(i,1).nudo_j.coordenadas(1,2) barras(i,1).nudo_k.coordenadas(1,2)];
                line(XY(1,:),XY(2,:),'Color','k','LineWidth',0.75)    
                pgon=polyshape(barras(i,1).A_internas_FE.Mi(1,:),barras(i,1).A_internas_FE.Mi(2,:),'Simplify',false);
                plot(pgon)
            end
            
            grid on
            title('Diagrama de Momentos Empotrados')
            
            subplot(4,3,9)
            daspect([1 1 1]);
            hold on
            
            for i=1:height(nudos)
                plot(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),'.','MarkerSize',15,'Color','k')
            end
            
            for i=1:height(barras)
                XY=[barras(i,1).nudo_j.coordenadas(1,1) barras(i,1).nudo_k.coordenadas(1,1)
                    barras(i,1).nudo_j.coordenadas(1,2) barras(i,1).nudo_k.coordenadas(1,2)];
                line(XY(1,:),XY(2,:),'Color','k','LineWidth',0.75)    
                pgon=polyshape(barras(i,1).A_internas_FE.Vi(1,:),barras(i,1).A_internas_FE.Vi(2,:),'Simplify',false);
                plot(pgon)
            end
            
            grid on
            title('Diagrama de Corte Empotrados')
            
            subplot(4,3,12)
            daspect([1 1 1]);
            hold on
            
            for i=1:height(nudos)
                plot(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),'.','MarkerSize',15,'Color','k')
            end
            
            for i=1:height(barras)
                XY=[barras(i,1).nudo_j.coordenadas(1,1) barras(i,1).nudo_k.coordenadas(1,1)
                    barras(i,1).nudo_j.coordenadas(1,2) barras(i,1).nudo_k.coordenadas(1,2)];
                line(XY(1,:),XY(2,:),'Color','k','LineWidth',0.75)    
                pgon=polyshape(barras(i,1).A_internas_FE.Ai(1,:),barras(i,1).A_internas_FE.Ai(2,:),'Simplify',false);
                plot(pgon)
            end
            
            grid on
            title('Diagrama de Axiales Empotrados')
            
            subplot(4,3,4)
            daspect([1 1 1]);
            hold on
            
            for i=1:height(nudos)
                plot(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),'.','MarkerSize',15,'Color','k')
            end
            
            for i=1:height(barras)
                XY=[barras(i,1).nudo_j.coordenadas(1,1) barras(i,1).nudo_k.coordenadas(1,1)
                    barras(i,1).nudo_j.coordenadas(1,2) barras(i,1).nudo_k.coordenadas(1,2)];
                line(XY(1,:),XY(2,:),'Color','k','LineWidth',0.75)    
                pgon=polyshape(barras(i,1).A_internas.Mi(1,:),barras(i,1).A_internas.Mi(2,:),'Simplify',false);
                plot(pgon)
            end
            
            grid on
            title('Diagrama de Momentos')
            
            subplot(4,3,7)
            daspect([1 1 1]);
            hold on
            
            for i=1:height(nudos)
                plot(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),'.','MarkerSize',15,'Color','k')
            end
            
            for i=1:height(barras)
                XY=[barras(i,1).nudo_j.coordenadas(1,1) barras(i,1).nudo_k.coordenadas(1,1)
                    barras(i,1).nudo_j.coordenadas(1,2) barras(i,1).nudo_k.coordenadas(1,2)];
                line(XY(1,:),XY(2,:),'Color','k','LineWidth',0.75)    
                pgon=polyshape(barras(i,1).A_internas.Vi(1,:),barras(i,1).A_internas.Vi(2,:),'Simplify',false);
                plot(pgon)
            end
            
            grid on
            title('Diagrama de Cortantes')
            
            subplot(4,3,10)
            daspect([1 1 1]);
            hold on
            
            for i=1:height(nudos)
                plot(nudos(i,1).coordenadas(1,1),nudos(i,1).coordenadas(1,2),'.','MarkerSize',15,'Color','k')
            end
            
            for i=1:height(barras)
                XY=[barras(i,1).nudo_j.coordenadas(1,1) barras(i,1).nudo_k.coordenadas(1,1)
                    barras(i,1).nudo_j.coordenadas(1,2) barras(i,1).nudo_k.coordenadas(1,2)];
                line(XY(1,:),XY(2,:),'Color','k','LineWidth',0.75)    
                pgon=polyshape(barras(i,1).A_internas.Ai(1,:),barras(i,1).A_internas.Ai(2,:),'Simplify',false);
                plot(pgon)
            end
            
            grid on
            title('Diagrama de Axiales')
        end
        
    end
end
% Calculo de la rigidez lateral de un portico
% Variacion del Peralte de vigas
clc
clear variables

format longG

% Sistema de unidades [MKS]

% Generamos los nudos y barras del problema
L=800; % Longitud del marco
H=300; % Altura del marco

% Parametros del material
HA=Material(15100*(240)^0.5,0.2);

% Parametros de las secciones
H_viga=0.01:0.1:150;

Column=Seccion(40*40,40*40^3/12,(5/6)*40*40);

Beam=repmat(Seccion,1,width(H_viga));
for i=1:width(H_viga)
    % Asumimos que la viga es infinitamente rigida
    Beam(i)=Seccion(10^5,30*H_viga(i)^3/12,(5/6)*30*H_viga(i));
end

% Creamos los nudos
nudos=[
       nudo(1,[0,0],["r" "r" "r"],[0 0 0])
       nudo(2,[L,0],["r" "r" "r"],[0 0 0])
       nudo(3,[0,H],["f" "f" "f"],[1000 0 0])
       nudo(4,[L,H],["f" "f" "f"],[0 0 0])
       ];

% Creamos las barras
barras=repmat(barra,3,width(H_viga));
for i=1:width(H_viga)
    barras(:,i)=[
                 barra(nudos(1,1),nudos(3,1),HA,Column,["na" "na"])
                 barra(nudos(2,1),nudos(4,1),HA,Column,["na" "na"])
                 barra(nudos(3,1),nudos(4,1),HA,Beam(i),["na" "na"])
                   ];
end

% Generamos el modelo

modelos=repmat(Modelo,1,width(H_viga));

for i=1:width(H_viga)
    modelos(i)=Modelo(nudos,barras(:,i));
    modelos(i)=modelos(i).Calculos();
    modelos(i)=modelos(i).Resultados_Nudos();
    modelos(i)=modelos(i).Resultados_Barras();
end

% Calculamos los resultados
delta=zeros(1,width(H_viga));
for i=1:width(H_viga)
    delta(i)=modelos(i).Nudos(4,1).Dj(1,1);
end

figure
hold on
plot(H_viga,delta,LineWidth=1.5)

grid on
title('Flexibilidad del Marco versus Peralte de Viga')
xlabel('Peralte de Viga [cm]')
ylabel('Deformacion [cm]')

filename='Flexibilidad';
saveas(gcf,filename,'jpg');

figure
hold on
plot(H_viga,delta/min(delta),LineWidth=1.5)

grid on
title('Flexibilidad del Marco versus Peralte de Viga')
xlabel('Peralte de Viga [cm]')
ylabel('Deformacion/Deformacion Minima')
ylim([0,4])

filename='Flexibilidad Relativa';
saveas(gcf,filename,'jpg');

figure
hold on
plot(H_viga,1./delta,LineWidth=1.5)

grid on
title('Rigidez del Marco versus Peralte de Viga')
xlabel('Peralte de Viga [cm]')
ylabel('Rigidez [kgf]')

filename='Rigidez';
saveas(gcf,filename,'jpg');

figure
hold on
plot(H_viga,1./(delta/max(delta)),LineWidth=1.5)

grid on
title('Rigidez del Marco versus Peralte de Viga')
xlabel('Peralte de Viga [cm]')
ylabel('Rigidez/Rigidez Minima')
ylim([0,4])

filename='Rigidez Relativa';
saveas(gcf,filename,'jpg');

Inversion de la Matriz