Contents

CÀLCUL D'INESTABILITAT PER AL PROGRAMA CME2.

MANUEL GIL RAMOS 2.013-2.042 GRAU EN ENGINYERIA DE LA CONSTRUCCIÓ UNIVERSITAT POLITÈCNICA DE CATALUNYA

Versió data 10/06/2.014

Entrada: n = sencer nombre de nusos de l'estructura. N = sencer nombre de barres de l'estructura. Ncoord = matriu n*3 amb les coordenades x, y i z dels nusos. Ncond = condicions de contorn dels nusos: nus lliure =0; rodet inclinat =1; recolzament simple =2; encastament =3; desplaçaments elàstics i girs lliures =4; girs elàstics i desplaçaments nuls =5; desplaçaments i girs elàstics =6; desplaçaments imposats i girs lliures =7; girs imposats i desplaçaments lliures =8; desplaçaments i girs imposats =9. Nextfor = vector n*6 de forces aplicades als nusos. Lnods = matriu N*2 amb nombre de nus inferior i superior de cada barra. Lartyp = matriu N*2 amb tipus d'extrem de barra inferior i superior, rígid =0; articulat =1. Lmecp = matriu N*9 de propietats mecaniques de cada barra: E, G, A, Iz, Iy, zx, zy, zz, J. barra_carrega_rep = matriu N*6 amb la càrrega repartida sobre les barres, mòdul, direcció x, y i z, distancia a extrem inferior i superior. barra_carrega_rep_proj = matriu N*6 amb la càrrega repartida projectada sobre les barres (longitut de la càrrega = longitut de la barra), mòdul, direcció x, y i z, distància a extrem inferior i superior. barra_carrega_punt = matriu N*5 amb la càrrega puntual sobre les barres, mòdul, direcció x, y i z i distància al nus inferior. barra_carrega_mom = matriu N*5 amb el moment aplicat sobre les barres, mòdul, direcció x, y i z i distància al nus inferior. iteracions = nombre d'iteracions a realitzar pel mètode iteratiu de resolució de 2on ordre alpha1 = primera variable del mètode de resolució Runge-Kutta alpha2 = segona variable del mètode de resolució Runge-Kutta mu = tercera variable del mètode de resolució Runge-Kutta metode = mètode de resolució a utilitzar pel càlcul de 2on ordre

Sortida: coef = vector amb els coeficients de vinclament ordenats per valor absolut. mode_vincl = matriu amb els moviments del nusos segons el mode de vinclament associat. TL2G = matriu 3*3*N amb les matrius de canvi de base local a global de cada barra. ierr = comptador d'errors que retorna un numero entre el 0 i el 4 segons el tipus d'error

function [coef,mode_vincl,TL2G,ierr] = calculinestabilitat(n,N,Ncoord,...
    Ncond,Nextfor,Lnods,Lartyp,Lmecp,barra_carrega_rep,...
    barra_carrega_rep_proj,barra_carrega_punt,barra_carrega_mom,...
    iteracions,alpha1,alpha2,mu,metode)

coef=zeros(1);

% Si tots els nusos tenen la mateixa coordenada x dimx=0.
dimx=any(Ncoord(:,1)~=Ncoord(1,1));

% Si tots els nusos tenen la mateixa coordenada y dimy=0.
dimy=any(Ncoord(:,2)~=Ncoord(1,2));

% Si tots els nusos tenen la mateixa coordenada z dimz=0.
dimz=any(Ncoord(:,3)~=Ncoord(1,3));

% Si existeix una barra col·locada en una direcció...
if (dimx&~dimy&~dimz)|(~dimx&dimy&~dimz)|(~dimx&~dimy&dimz)

    % Si la barra té les mateixes coordenades x i z (barra en eix y)
    if ~any(Ncoord(:,1)~=Ncoord(1,1)) && ~any(Ncoord(:,3)~=Ncoord(1,3))

        Ncoord(n+1,1)=Ncoord(1,1)+1; % Creació d'un punt imaginari fora de
                                     % l'eix de la barra per evitar
                                     % eliminar certs moviments

    % Si la barra té les mateixes coordenades y i z (barra en eix x)
    elseif ~any(Ncoord(:,2)~=Ncoord(1,2)) && ~any(Ncoord(:,3)~=Ncoord(1,3))

        Ncoord(n+1,2)=Ncoord(1,2)+1; % Creació d'un punt imaginari fora de
                                     % l'eix de la barra per evitar
                                     % eliminar certs moviments

    % Si la barra té les mateixes coordenades z (barra en eix x o y)
    elseif ~any(Ncoord(:,3)~=Ncoord(1,3))
        ierr=5; % Com les estructures bidimensionals han d'estar als
                % eixos x-y, retorna un error

        TL2G=[]; % Dimensiona la matriu de canvi de base buida

        mode_vincl=[]; % Dimensiona la matriu amb els modes de vinclament
                       % buida
    end
end

% Si existeixen càrregues en barra (repartides, moments, puntuals...)
if any(barra_carrega_rep(:,1))|any(barra_carrega_rep_proj(:,1))|...
        any(barra_carrega_punt(:,1))|any(barra_carrega_mom(:,1))

    ierr=3; % Error que indica càrregues en barra
    TL2G=[]; % Dimensiona la matriu de canvi de base buida
    mode_vincl=[]; % Dimensiona la matriu amb els modes de vinclament
                   % buida
    return; % Retorna l'execució a la funció cme2.m
end

% Cridem a la funció del càlcul estàtic amb els arguments de sortida del
% nostre interès
[~,~,~,TL2G,ierr,K_Glob,Vector_esforcos,Vector_esforcos_prescrits,...
    sistema,ns,nn] = cme2_calcul_estatic(n,N,Ncoord,Ncond,Nextfor,...
    Lnods,Lartyp,Lmecp,barra_carrega_rep,barra_carrega_rep_proj,...
    barra_carrega_punt,barra_carrega_mom);

K_Geom_Glob_cte=zeros(ns,ns); % Dimensionem la matriu geomètrica global
                              % perquè el programa s'executi mès ràpid

if ierr==1 % Si l'execució del càlcul estàtic ens retorna un error...
    mode_vincl=[]; % Dimensiona la matriu amb els modes de vinclament
                   % buida
    return % Retorna l'execució a la funció cme2.m
end

kriggeom_cte_tot=zeros(12,12,N); % Dimensiona la hipermatriu de matrius
                                 % geomètriques locals en eixos locals
kriggeom_cte=zeros(12,12,N); % Dimensiona la hipermatriu de matrius
                             % geomètriques locals en eixos locals

for i=1:N % Per a cada barra...
    n1=Lnods(i,1); % node inferior
    n2=Lnods(i,2); % node superior
    extrems=Lartyp(i,1)*2+Lartyp(i,2); % rígid-rígid =0
                                       % rígid-articulat =1
                                       % articulat-rígid =2
                                       % articulat-articulat =3

    vector_barra=Ncoord(n2,1:3)-Ncoord(n1,1:3); % Vector de la barra.
    L=norm(vector_barra); % Longitut de la barra.
    T=[TL2G(:,:,i),zeros(3);zeros(3),TL2G(:,:,i)]; % Matriu canvi de base
                                                   % completa 6x6
    [kriggeom_cte(:,:,i)]=calcul_submatrius_geometriques(Lmecp(i,:),...
        extrems,L); % Càlcul de les submatrius geomètriques en eixos
                    % locals.
    kriggeom_cte_tot(:,:,i)=[T*kriggeom_cte(1:6,1:6,i)*T',T*...
        kriggeom_cte(1:6,7:12,i)*T';T*kriggeom_cte(7:12,1:6,i)*T'...
        ,T*kriggeom_cte(7:12,7:12,i)*T']; % Submatriu geomètrica en eixos
                                          % globals.
    kriggeom_cte_tot(:,:,i)=contorn_rodet_inclinat(kriggeom_cte_tot...
        (:,:,i),n1,n2,Ncond); % Aplicació de la condició de contorn de
                              % rodet inclinat a la submatriu geometrica,
                              % si s'escau.
end

Vector_moviments=zeros(ns,1); % Dimensiona el vector de moviments per
                              % optimitzar l'execució
incr=iteracions; % Li dona a la variable incr el valor de les iteracions
                 % introduides a l'interfície gràfica

switch metode % Switch per escollir mètode
    case 1 % Cas del mètode 1 (Euler)
        [K_Geom_Glob_cte,errmet]=Euler_iterative...
            (Vector_esforcos_prescrits,K_Glob,TL2G,K_Geom_Glob_cte,...
            sistema,Vector_moviments,Vector_esforcos,n,nn,ns,Ncond,N,...
            barra_carrega_rep,barra_carrega_rep_proj,barra_carrega_punt,...
            barra_carrega_mom,kriggeom_cte_tot,kriggeom_cte,Lnods,...
            Lartyp,Lmecp,Ncoord,incr);
    case 2 % Cas del mètode 2 (Runge-Kutta)
        [K_Geom_Glob_cte,errmet]=Runge_Kutta(Vector_esforcos_prescrits,...
            K_Glob,TL2G,K_Geom_Glob_cte,sistema,Vector_moviments,...
            Vector_esforcos,n,nn,ns,Ncond,N,barra_carrega_rep,...
            barra_carrega_rep_proj,barra_carrega_punt,barra_carrega_mom,...
            kriggeom_cte_tot,kriggeom_cte,Lnods,Lartyp,Lmecp,Ncoord,...
            incr,alpha1,alpha2,mu);
end

if errmet==1 % Si l'execució dels mètodes ens ha retornat l'error de canvi
             % de signe

    ierr=4; % Dona a la variable ierr el valor pertinent
    TL2G=[]; % Dimensiona la matriu de canvi de base buida
    mode_vincl=[]; % Dimensiona la matriu amb els modes de vinclament buida
    return; % Retorna l'execució a la funció cme2.m
end

B=K_Geom_Glob_cte(1:nn,1:nn); % Dimensiona la variable B amb la part de
                             % la matriu de rigidesa que ens interessa
A=K_Glob(1:nn,1:nn); % Dimensiona la variable A amb la part de la matriu
                     % de rigidesa global que ens interessa

[vect,mat_coefcoef]=eig(A,B); % Resolució del problema d'autovalors amb
                              % diverses variables de retorn pels
                              % coeficients i els modes de vinclament

dimvec=size(vect); % Es troba el nombre de moviments determinats pel mode
                   % de vinclament

if dimvec(1)<6*n % Si el nombre de moviments determinats pel mode de
                 % vinclament és inferior al total...

    vect((dimvec(1)+1):(6*n),:)=0; % La resta de moviments són igual a 0

end

dimvec=size(vect); % Es troba el nombre de moviments totals

vect_reord=zeros(dimvec(1),1);

for i=1:nn % Bucle que reordena els vectors de la matriu amb els modes
           % de vinclament. Col·loca primer el desplaçament en x del node
           % 1, el desplaçament en y, etc.

    vect_reord(:,i)=reordenacio_vec(vect(:,i),sistema,n,Ncond);

end

i=1;

while i<=nn % Bucle per extreure la diagonal de la matriu de coeficients

    vect_coef(i,1)=mat_coefcoef(i,i);
    i=i+1;

end

vect_coef_abs=abs(vect_coef); % Vector amb els valors absoluts
[~,ord]=sort(vect_coef_abs); % Busquem l'ordre del vector absolut

for i=1:max(ord) % Bucle que ordena els coeficients i modes de vinclament
                 % segons el valor absolut dels coeficients

    coeford(i,1)=vect_coef(ord(i),1);
    vect_vincl(:,i)=vect_reord(:,ord(i));

end

if all(isinf(coeford))|all(isnan(coeford)) % Si tots els coeficients
                                           % són infinit...
    ierr=2; % Dona a la variable ierr el valor pertinent
    mode_vincl=[]; % Dimensiona la matriu amb els modes de vinclament buida
    return; % Retorna l'execució a la funció cme2.m
end

k=1;

for i=1:max(ord) % Bucle que elimina els coeficients que són infinit
                 % (no són d'interès)
    if ~isinf(coeford(i))

        coef(k)=coeford(i);
        mode_vincl(:,k)=vect_vincl(:,i);
        k=k+1;

    end
end

end
Input argument "Ncoord" is undefined.

Error in ==> calculinestabilitat at 58
dimx=any(Ncoord(:,1)~=Ncoord(1,1)); 

SUBMATRIUS GEOMETRIQUES EN EIXOS LOCALS.

% Funció que calcula les matrius geomètriques locals segons els tipus
% d'extrems de les barres.

function [kriggeom_cte]=calcul_submatrius_geometriques(Lmecp,extrems,l)

     switch extrems
         case 0   %si la barra es rigida en tots dos extrems
            kriggeom_cte(1,1)=0;             %k11
            kriggeom_cte(2,2)=36;
            kriggeom_cte(2,6)=3*l;
            kriggeom_cte(3,3)=36;
            kriggeom_cte(3,5)=-3*l;
            kriggeom_cte(4,4)=0;
            kriggeom_cte(5,3)=-3*l;
            kriggeom_cte(5,5)=4*l^2;
            kriggeom_cte(6,2)=3*l;
            kriggeom_cte(6,6)=4*l^2;

            kriggeom_cte(1,7)=0;            %k12
            kriggeom_cte(2,8)=-36;
            kriggeom_cte(2,12)=3*l;
            kriggeom_cte(3,9)=-36;
            kriggeom_cte(3,11)=-3*l;
            kriggeom_cte(4,10)=0;
            kriggeom_cte(5,9)=-3*l;
            kriggeom_cte(5,11)=-l^2;
            kriggeom_cte(6,8)=-3*l;
            kriggeom_cte(6,12)=-l^2;

            kriggeom_cte(7,1)=0;             %k21
            kriggeom_cte(8,2)=-36;
            kriggeom_cte(8,6)=-3*l;
            kriggeom_cte(9,3)=-36;
            kriggeom_cte(9,5)=3*l;
            kriggeom_cte(10,4)=0;
            kriggeom_cte(11,3)=-3*l;
            kriggeom_cte(11,5)=-l^2;
            kriggeom_cte(12,2)=3*l;
            kriggeom_cte(12,6)=-l^2;

            kriggeom_cte(7,7)=0;             %k22
            kriggeom_cte(8,8)=36;
            kriggeom_cte(8,12)=-3*l;
            kriggeom_cte(9,9)=36;
            kriggeom_cte(9,11)=3*l;
            kriggeom_cte(10,10)=0;
            kriggeom_cte(11,9)=3*l;
            kriggeom_cte(11,11)=4*l^2;
            kriggeom_cte(12,8)=-3*l;
            kriggeom_cte(12,12)=4*l^2;


        case 1 %si la barra es articulada en l'extrem inferior

            kriggeom_cte(1,1)=0;
            kriggeom_cte(2,2)=36;
            kriggeom_cte(2,6)=0;
            kriggeom_cte(3,3)=36;
            kriggeom_cte(3,5)=0;
            kriggeom_cte(4,4)=0;
            kriggeom_cte(5,3)=-0;
            kriggeom_cte(5,5)=0;
            kriggeom_cte(6,2)=0;
            kriggeom_cte(6,6)=0;

            kriggeom_cte(1,7)=0;
            kriggeom_cte(2,8)=-36;
            kriggeom_cte(2,12)=6*l;
            kriggeom_cte(3,9)=-36;
            kriggeom_cte(3,11)=-6*l;
            kriggeom_cte(4,10)=-0;
            kriggeom_cte(5,9)=0;
            kriggeom_cte(5,11)=0;
            kriggeom_cte(6,8)=-0;
            kriggeom_cte(6,12)=0;

            kriggeom_cte(7,1)=0;
            kriggeom_cte(8,2)=-36;
            kriggeom_cte(8,6)=0;
            kriggeom_cte(9,3)=-36;
            kriggeom_cte(9,5)=0;
            kriggeom_cte(10,4)=-0;
            kriggeom_cte(11,3)=-6*l;
            kriggeom_cte(11,5)=0;
            kriggeom_cte(12,2)=6*l;
            kriggeom_cte(12,6)=0;

            kriggeom_cte(7,7)=0;
            kriggeom_cte(8,8)=36;
            kriggeom_cte(8,12)=-6*l;
            kriggeom_cte(9,9)=36;
            kriggeom_cte(9,11)=6*l;
            kriggeom_cte(10,10)=0;
            kriggeom_cte(11,9)=6*l;
            kriggeom_cte(11,11)=6*l^2;
            kriggeom_cte(12,8)=-6*l;
            kriggeom_cte(12,12)=6*l^2;
        case 2     %si la barra es articulada en l'extrem superior
            kriggeom_cte(1,1)=0;
            kriggeom_cte(2,2)=36;
            kriggeom_cte(2,6)=6*l;
            kriggeom_cte(3,3)=36;
            kriggeom_cte(3,5)=-6*l;
            kriggeom_cte(4,4)=0;
            kriggeom_cte(5,3)=-6*l;
            kriggeom_cte(5,5)=6*l^2;
            kriggeom_cte(6,2)=6*l;
            kriggeom_cte(6,6)=6*l^2;

            kriggeom_cte(1,7)=0;
            kriggeom_cte(2,8)=-36;
            kriggeom_cte(2,12)=0;
            kriggeom_cte(3,9)=-36;
            kriggeom_cte(3,11)=-0;
            kriggeom_cte(4,10)=-0;
            kriggeom_cte(5,9)=6*l;
            kriggeom_cte(5,11)=0;
            kriggeom_cte(6,8)=-6*l;
            kriggeom_cte(6,12)=0;

            kriggeom_cte(7,1)=0;
            kriggeom_cte(8,2)=-36;
            kriggeom_cte(8,6)=-6*l;
            kriggeom_cte(9,3)=-36;
            kriggeom_cte(9,5)=-6*l;
            kriggeom_cte(10,4)=-0;
            kriggeom_cte(11,3)=-0;
            kriggeom_cte(11,5)=0;
            kriggeom_cte(12,2)=0;
            kriggeom_cte(12,6)=0;

            kriggeom_cte(7,7)=0;
            kriggeom_cte(8,8)=36;
            kriggeom_cte(8,12)=-0;
            kriggeom_cte(9,9)=36;
            kriggeom_cte(9,11)=0;
            kriggeom_cte(10,10)=0;
            kriggeom_cte(11,9)=0;
            kriggeom_cte(11,11)=0;
            kriggeom_cte(12,8)=-0;
            kriggeom_cte(12,12)=0;

        case 3   %si la barra es articulada en tots dos extrems
            kriggeom_cte(1,1)=0;
            kriggeom_cte(2,2)=30;
            kriggeom_cte(2,6)=0;
            kriggeom_cte(3,3)=30;
            kriggeom_cte(3,5)=0;
            kriggeom_cte(4,4)=0;
            kriggeom_cte(5,3)=0;
            kriggeom_cte(5,5)=0;
            kriggeom_cte(6,2)=0;
            kriggeom_cte(6,6)=0;

            kriggeom_cte(1,7)=0;
            kriggeom_cte(2,8)=-30;
            kriggeom_cte(2,12)=0;
            kriggeom_cte(3,9)=-30;
            kriggeom_cte(3,11)=-0;
            kriggeom_cte(4,10)=-0;
            kriggeom_cte(5,9)=0;
            kriggeom_cte(5,11)=0;
            kriggeom_cte(6,8)=0;
            kriggeom_cte(6,12)=0;

            kriggeom_cte(7,1)=0;
            kriggeom_cte(8,2)=-30;
            kriggeom_cte(8,6)=0;
            kriggeom_cte(9,3)=-30;
            kriggeom_cte(9,5)=0;
            kriggeom_cte(10,4)=-0;
            kriggeom_cte(11,3)=-0;
            kriggeom_cte(11,5)=0;
            kriggeom_cte(12,2)=0;
            kriggeom_cte(12,6)=0;

            kriggeom_cte(7,7)=0;
            kriggeom_cte(8,8)=30;
            kriggeom_cte(8,12)=-0;
            kriggeom_cte(9,9)=30;
            kriggeom_cte(9,11)=0;
            kriggeom_cte(10,10)=0;
            kriggeom_cte(11,9)=0;
            kriggeom_cte(11,11)=0;
            kriggeom_cte(12,8)=-0;
            kriggeom_cte(12,12)=0;

    end

end

SUBMATRIUS DE REGIDESA EN EIXOS LOCALS.

% Funció que calcula les matrius de rigidesa en eixos locals. Al cridar a
% la funció de càlcul estàtic nomès obtenim la matriu de rigidesa global.

function KL=calcul_submatrius_rigidesa(Lmecp,extrems,L)
E=Lmecp(1); % Mòdul de Young
Iz=Lmecp(4); % Inèrcia eix z
if any(Lmecp(6:8)) % Si existeix vector direcció d'Iz
    Iy=Lmecp(5); % Inèrcia eix y
else
    Iy=Iz;
end
GItL=Lmecp(9)/L; % Inèrcia a torsió
A=Lmecp(3); % Àrea de la secció
EAL=E*A/L;
EIzL=E*Iz/L;
EIyL=E*Iy/L;
EIzL2=E*Iz/L^2;
EIyL2=E*Iy/L^2;
EIzL3=E*Iz/L^3;
EIyL3=E*Iy/L^3;
% Sub-matrius de rigidesa segons els extrems siguin rígids o articulats.
switch extrems
    case 0 % rígid-rígid
        KL(1:6,1:6)=[EAL,0,0,0,0,0;0,12*EIzL3,0,0,0,6*EIzL2;0,0,...
            12*EIyL3,0,-6*EIyL2,0;0,0,0,GItL,0,0;0,0,-6*EIyL2,0,4*...
            EIyL,0;0,6*EIzL2,0,0,0,4*EIzL];
        KL(1:6,7:12)=[-EAL,0,0,0,0,0;0,-12*EIzL3,0,0,0,6*EIzL2;0,0,...
            -12*EIyL3,0,-6*EIyL2,0;0,0,0,-GItL,0,0;0,0,6*EIyL2,0,2*EIyL,...
            0;0,-6*EIzL2,0,0,0,2*EIzL];
        KL(7:12,1:6)=KL(1:6,7:12)';
        KL(7:12,7:12)=[EAL,0,0,0,0,0;0,12*EIzL3,0,0,0,-6*EIzL2;0,0,...
            12*EIyL3,0,6*EIyL2,0;0,0,0,GItL,0,0;0,0,6*EIyL2,0,4*EIyL,0;...
            0,-6*EIzL2,0,0,0,4*EIzL];
     case 1 % rígid-articulat
        KL(1:6,1:6)=[EAL,0,0,0,0,0;0,3*EIzL3,0,0,0,3*EIzL2;0,0,3*EIyL3,...
            0,-3*EIyL2,0;0,0,0,0,0,0;0,0,-3*EIyL2,0,3*EIyL,0;0,3*EIzL2,...
            0,0,0,3*EIzL];
        KL(1:6,7:12)=[[-EAL,0,0;0,-3*EIzL3,0;0,0,-3*EIyL3;0,0,0;0,0,3*...
            EIyL2;0,-3*EIzL2,0],zeros(6,3)];
        KL(7:12,1:6)=KL(1:6,7:12)';
        KL(7:12,7:12)=[[EAL,0,0;0,3*EIzL3,0;0,0,3*EIyL3],zeros(3,3);...
            zeros(3,6)];
    case 2 % articulat-rígid
        KL(1:6,1:6)=[[EAL,0,0;0,3*EIzL3,0;0,0,3*EIyL3],zeros(3);...
            zeros(3,6)];
        KL(1:6,7:12)=[[-EAL,0,0,0,0,0;0,-3*EIzL3,0,0,0,3*EIzL2;0,0,...
            -3*EIyL3,0,-3*EIyL2,0];zeros(3,6)];
        KL(7:12,1:6)=KL(1:6,7:12)';
        KL(7:12,7:12)=[EAL,0,0,0,0,0;0,3*EIzL3,0,0,0,-3*EIzL2;0,0,...
            3*EIyL3,0,3*EIyL2,0;0,0,0,0,0,0;0,0,3*EIyL2,0,3*EIyL,0;0,...
            -3*EIzL2,0,0,0,3*EIzL];
    case 3 % articulat-articulat
        KL(1:6,1:6)=[EAL,0,0,0,0,0;zeros(5,6)];
        KL(1:6,7:12)=[-EAL,0,0,0,0,0;zeros(5,6)];
        KL(7:12,1:6)=KL(1:6,7:12);
        KL(7:12,7:12)=KL(1:6,1:6);
end
end

CÀLCUL DE LES SOL·LICITACIONS ALS NUSOS PER CÀRREGUES A LES BARRES.

%  Per a cada tipus d'extrem de barra, rígid o articulat, en calcula els
%  esforços d'empotrament perfecte, per a càrregues repartides projectades
%  (longitut de la càrrega = longitut de la barra) o no sobre la barra,
%  càrregues puntuals i moments aplicats a la bara. Per a càrregues
%  repartides permet definir distàncies al nus inferior i superior per
%  delimitar la zona de la barra que rep la càrrega. Igualment amb la
%  distància al nus inferior es situa la càrrega puntual i el moment a la
%  barra. Totes les càrregues s'introdueixen, i els resultats es retornen,
%  en eixos globals, però es calculen en eixos locals (eix de la barra
%  coincideix amb l'eix X local).


function reaccions=calcul_carregues_barra_local(carrega_rep,...
    carrega_rep_proj,carrega_punt,carrega_mom,TL2G,L,extrems)
reaccions=[0,0,0,0,0,0,0,0,0,0,0,0];
if carrega_rep_proj(1) % Càlcul de les forces equivalents als nusos a la
                       % càrrega repartida projectada sobre la barra.
    vbarra=[L,0,0]; % Vector barra en eixos locals;
    carrega=TL2G'*carrega_rep_proj(2:4)'; % Vector direcció de la càrrega
                                          % en eixos locals.
    carrega=carrega_rep_proj(1)*carrega/norm(carrega); % Vector càrrega en
                                                       % eixos locals.
    vperp=cross(vbarra,carrega);
    vperp=cross(vperp,carrega);
    vperp=vperp/norm(vperp);
    Lc=abs(dot(vperp,vbarra)); % Longitut de la càrrega.
    reaccions=-calcul_repartida(carrega,carrega_rep_proj(5),...
        carrega_rep_proj(6),L,Lc,extrems);
end
if carrega_rep(1) % Càlcul de les forces equivalents als nusos a la
                  % càrrega repartida a la barra.
    carrega=TL2G'*carrega_rep(2:4)'; % Vector direcció de la càrrega
                                     % repartida en eixos locals.
    carrega=carrega_rep(1)*carrega'/norm(carrega); % Vector càrrega
    reaccions=reaccions-calcul_repartida(carrega,carrega_rep(5),...
        carrega_rep(6),L,L,extrems);
end
if carrega_punt(1) % Càlcul de les forces equivalents als nusos a la
                   % càrrega puntual a la barra.
    carrega=TL2G'*carrega_punt(2:4)'; % Vector direcció de la carrega
                                      % puntual
    carrega=carrega'*carrega_punt(1)/norm(carrega); % Vector carrega.
    reaccions=reaccions-calcul_puntual(carrega,carrega_punt(5),L,extrems);
end
if carrega_mom(1) % Càlcul de les forces equivalents als nusos al moment
                  % aplicat a la barra.
    carrega=TL2G'*carrega_mom(2:4)'; % Vector direcció del moment.
    carrega=carrega'*carrega_mom(1)/norm(carrega); % Vector moment.
    reaccions=reaccions-calcul_moment(carrega,carrega_mom(5),L,extrems);
end
end

APLICACIÓ CONDICIÓ DE CONTORN RODET INCLINAT A LES MATRIUS DE RIGIDESA.

% Condició de contorn rodet inclinat.
% Si algún dels extrems connecta amb un rodet inclinat, multiplica les
% submatrius per la matriu de canvi de base als eixos del rodet, els
% moviments i esforços del nus corresponent queden en eixos locals del
% rodet.

function KG=contorn_rodet_inclinat(KG,n1,n2,Ncond)
if Ncond(n1,1)==1 % Si el nus inferior és un rodet inclinat
    Tr=calcul_canvi_base(Ncond(n1,2:4)); % Calcula matriu canvi de base
                                         % del rodet.
    T=[Tr,zeros(3);zeros(3),Tr];
    KG(1:6,1:6)=T'*KG(1:6,1:6)*T; % Sub-matriu en la base del rodet.
    KG(1:6,7:12)=T'*KG(1:6,7:12); % Sub-matriu en la base del rodet.
    KG(7:12,1:6)=KG(7:12,1:6)*T; % Sub-matriu en la base del rodet.
end
if Ncond(n2,1)==1 % Si el nus superior és un rodet inclinat
    Tr=calcul_canvi_base(Ncond(n2,2:4)); % Calcula matriu canvi de base
                                         % del rodet.
    T=[Tr,zeros(3);zeros(3),Tr];
    KG(1:6,7:12)=KG(1:6,7:12)*T; % Sub-matriu en la base del rodet.
    KG(7:12,1:6)=T'*KG(7:12,1:6); % Sub-matriu en la base del rodet.
    KG(7:12,7:12)=T'*KG(7:12,7:12)*T; % Sub-matriu en la base del rodet.
end
end

RESOLUCIÓ DEL PROBLEMA NO LINEAL MITJANÇANT EL MÈTODE D'EULER

% Funció que retorna la matriu geomètrica global calculada pel mètode
% d'Euler

function [K_Geom_Glob_cte,errmet]=Euler_iterative...
    (Vector_esforcos_prescrits,K_Glob,TL2G,K_Geom_Glob_cte,...
    sistema,Vector_moviments,Vector_esforcos,n,nn,ns,Ncond,N,...
    barra_carrega_rep,barra_carrega_rep_proj,barra_carrega_punt,...
    barra_carrega_mom,kriggeom_cte_tot,kriggeom_cte,Lnods,Lartyp,Lmecp,...
    Ncoord,incr)

errmet=0;
dland=1/incr; % La inversa dels increments ens servirà per trobar el
              % diferencial de càrrega
Ax=zeros(N,1); % Dimensionament d'una matriu que emmagatzemarà les dàdes
               % dels axils

dp=dland*Vector_esforcos_prescrits; % Diferencial de càrrega (utilitzant
                                    % la inversa dels increments)

% Dimensiona a zero totes les dades que utlitzarem per optimitzar
% l'execució
K_Geom_Glob_cte=K_Geom_Glob_cte*0; % Matriu geomètrica global
sum_dif_mov=Vector_moviments*0; % Moviments totals
d_vector_moviments=Vector_moviments*0; % Moviments diferencials
kriggeom_loc_ev=kriggeom_cte*0; % Matriu gemètrica local

for j=1:incr;

    % Troba el diferencial de moviment d'aquest pas
    d_vector_moviments(1:nn,1)=(K_Glob(1:nn,1:nn)-K_Geom_Glob_cte...
        (1:nn,1:nn))\dp(1:nn,1);

    while det(K_Glob(1:nn,1:nn)-K_Geom_Glob_cte(1:nn,1:nn))<0 % Si el
                                      % determinant de la matriu geomètrica
                                      % global es superior al determinant
                                      % de la matriu de rigidesa global

        K_Geom_Glob_cte=[]; % Dimensionament de la matriu goemètrica global
                            % buida
        errmet=1; % Valor pertinent per la variable d'error del mètode
        return % Retorna l'execució a la funció calculinestabilitat.m
    end

    sum_dif_mov=sum_dif_mov+d_vector_moviments; % Sumatori de moviments

    % A continuació e calcularàn els esforços en les barres per trobar
    % l'axil de casascuna d'elles i introduir-lo a les matrius geomètriques
    % locals

    % Calcula les reaccions als recolzaments de l'estructura.
    Vector_reaccions=K_Glob*sum_dif_mov-Vector_esforcos;

    % Retorn de reaccions i moviments en l'ordre original dels nusos i
    % eixos globals.
    [Retorn_reaccions,Retorn_moviments]=ordre_inicial_nusos...
        (sistema,Vector_reaccions,sum_dif_mov,n,Ncond); %

    Retorn_esforcos_barres=zeros(N,12); % Càlcul dels esforços en extrems
                                        % de barra.
    K_Geom_Glob_cte=zeros(ns); % Dimensiona a zero la matriu geomètrica
                               % global per evitar sumar l'actual a
                               % l'anterior.

for i=1:N
    n1=Lnods(i,1); % node inferior
    n2=Lnods(i,2); % node superior
    extrems=Lartyp(i,1)*2+Lartyp(i,2); % rígid-rígid =0
                                       % rígid-articulat =1
                                       % articulat-rígid =2
                                       % articulat-articulat =3

    vector_barra=Ncoord(n2,1:3)-Ncoord(n1,1:3); % Vector de la barra.
    L=norm(vector_barra); % Longitut de la barra.
    T=[TL2G(:,:,i),zeros(3);zeros(3),TL2G(:,:,i)]; % Matriu canvi de base
                                                   % completa 6x6

    % Càlcul de la submatriu de rigidesa en eixos locals.
    KL=calcul_submatrius_rigidesa(Lmecp(i,:),extrems,L);
    vector_mov_local=[T'*Retorn_moviments(n1*6-5:n1*6);T'*...
        Retorn_moviments(n2*6-5:n2*6)]; % Vector moviments en eixos locals.
    Nepfor(i,:)=calcul_carregues_barra_local(barra_carrega_rep...
        (i,:),barra_carrega_rep_proj(i,:),barra_carrega_punt(i,:),...
        barra_carrega_mom(i,:),TL2G(:,:,i),L,extrems);
    Retorn_esforcos_barres(i,:)=(KL-kriggeom_loc_ev(:,:,i))*...
        vector_mov_local-Nepfor(i,:)'; % Esforços en extrems de cada barra
                                       % en eixos locals.
    Ax(i,j+1)=Retorn_esforcos_barres(i,1); % Emmagatzema els axils de barra

    % A coninuació evaula la matriu geomètrica local i l'acobla a la matriu
    % geomètrica global. Aquesta matriu geomètrica serà la que utilitzarem
    % al pas següent.

    kriggeom_cte_ev(:,:,i)=(Ax(i,j+1)/(30*L))*kriggeom_cte_tot(:,:,i);
    kriggeom_loc_ev(:,:,i)=(Ax(i,j+1)/(30*L))*kriggeom_cte(:,:,i);
    K_Geom_Glob_cte=K_Geom_Glob_cte+ensamblatge_matriu_rigidesa...
        (kriggeom_cte_ev(:,:,i),ns,sistema,n1,n2);
end

end

    % Un cop s'han realitzat tots els increments trobarem la matriu
    % geomètrica global a partir dels desplaçaments finals dels nussos.


    % Calcula les reaccions als recolzaments de l'estructura.
    Vector_reaccions=K_Glob*sum_dif_mov-Vector_esforcos;

    % Retorn de reaccions i moviments en l'ordre original dels nusos i
    % eixos globals.
    [Retorn_reaccions,Retorn_moviments]=ordre_inicial_nusos(sistema,...
        Vector_reaccions,sum_dif_mov,n,Ncond);

    Retorn_esforcos_barres=zeros(N,12); % Càlcul dels esforços en extrems
                                        % de barra.
    K_Geom_Glob_cte=zeros(ns);

    for i=1:N
    n1=Lnods(i,1); % node inferior
    n2=Lnods(i,2); % node superior
    extrems=Lartyp(i,1)*2+Lartyp(i,2); % rígid-rígid =0
                                       % rígid-articulat =1
                                       % articulat-rígid =2
                                       % articulat-articulat =3
    vector_barra=Ncoord(n2,1:3)-Ncoord(n1,1:3); % Vector de la barra.
    L=norm(vector_barra); % Longitut de la barra.
    T=[TL2G(:,:,i),zeros(3);zeros(3),TL2G(:,:,i)]; % Matriu canvi de base
                                                   % completa 6x6

    % Càlcul de la submatriu de rigidesa en eixos locals.
    KL=calcul_submatrius_rigidesa(Lmecp(i,:),extrems,L);
    vector_mov_local=[T'*Retorn_moviments(n1*6-5:n1*6);T'*...
        Retorn_moviments(n2*6-5:n2*6)]; % Vector moviments en eixos locals.
    Nepfor(i,:)=calcul_carregues_barra_local(barra_carrega_rep(i,:),...
        barra_carrega_rep_proj(i,:),barra_carrega_punt(i,:),...
        barra_carrega_mom(i,:),TL2G(:,:,i),L,extrems);
    Retorn_esforcos_barres(i,:)=(KL-kriggeom_loc_ev(:,:,i))*...
        vector_mov_local-Nepfor(i,:)'; % Esforços en extrems de cada barra
                                       % en eixos locals.
    Ax_fin(i,1)=Retorn_esforcos_barres(i,1);

    % A coninuació evaula la matriu geomètrica local i l'acobla a la matriu
    % geomètrica global. Aquesta matriu geomètrica serà la que la funció
    % retornarà.
    kriggeom_cte_ev(:,:,i)=(Ax_fin(i,1)/(30*L))*kriggeom_cte_tot(:,:,i);
    kriggeom_loc_ev(:,:,i)=(Ax_fin(i,1)/(30*L))*kriggeom_cte(:,:,i);
    K_Geom_Glob_cte=K_Geom_Glob_cte+ensamblatge_matriu_rigidesa...
        (kriggeom_cte_ev(:,:,i),ns,sistema,n1,n2);

end

end

RESOLUCIÓ DEL PROBLEMA NO LINEAL MITJANÇANT EL MÈTODE DE RUNGE_KUTTA

% Funció que retorna la matriu geomètrica global calculada pel mètode
% de Runge-Kutta

function [K_Geom_Glob_cte,errmet]=Runge_Kutta(Vector_esforcos_prescrits,...
    K_Glob,TL2G,K_Geom_Glob_cte,sistema,Vector_moviments,...
    Vector_esforcos,n,nn,ns,Ncond,N,barra_carrega_rep,...
    barra_carrega_rep_proj,barra_carrega_punt,barra_carrega_mom,...
    kriggeom_cte_tot,kriggeom_cte,Lnods,Lartyp,Lmecp,Ncoord,incr,alpha1,...
    alpha2,mu)

errmet=0;
dland=1/incr; % La inversa dels increments ens servirà per trobar el
              % diferencial de càrrega
Ax=zeros(N,1); % Dimensionament d'una matriu que emmagatzemarà les dàdes
               % dels axils
Ax_prim=zeros(N,1); % Dimensionament d'una matriu que emmagatzemarà les
                    % dàdes dels axils dels passos intermedis
dp=dland*Vector_esforcos_prescrits; % Diferencial de càrrega (utilitzant
                                    % la inversa dels increments)

% Dimensiona a zero totes les dades que utlitzarem per optimitzar
% l'execució
K_Geom_Glob_cte=K_Geom_Glob_cte*0;
K_Geom_Glob_cte_prim=K_Geom_Glob_cte*0;
sum_dif_mov=Vector_moviments*0;
d_vector_moviments=Vector_moviments*0;
d_vector_moviments_prim=Vector_moviments*0;
kriggeom_loc_ev=kriggeom_cte*0;


for j=1:incr;

    % Busca la matriu de rigidesa tangent a l'inici del pas a partir
    % dels moviments totals dels nussos al pas anteior. Si és el primer
    % pas, la matriu geomètrica serà igual a 0 i per tant la matriu de
    % rigidesa tangent serà igual a la matriu de rigidesa global.

    % Calcula les reaccions als recolzaments de l'estructura.
    Vector_reaccions=K_Glob*sum_dif_mov-Vector_esforcos;

    % Retorn de reaccions i moviments en l'ordre original dels nusos i
    % eixos globals.
    [Retorn_reaccions,Retorn_moviments]=ordre_inicial_nusos...
        (sistema,Vector_reaccions,sum_dif_mov,n,Ncond);

    Retorn_esforcos_barres=zeros(N,12); % Càlcul dels esforços en extrems
                                        % de barra.
    K_Geom_Glob_cte=zeros(ns);

for i=1:N
    n1=Lnods(i,1); % node inferior
    n2=Lnods(i,2); % node superior
    extrems=Lartyp(i,1)*2+Lartyp(i,2); % rígid-rígid =0
                                       % rígid-articulat =1
                                       % articulat-rígid =2
                                       % articulat-articulat =3
    vector_barra=Ncoord(n2,1:3)-Ncoord(n1,1:3); % Vector de la barra.
    L=norm(vector_barra); % Longitut de la barra.
    T=[TL2G(:,:,i),zeros(3);zeros(3),TL2G(:,:,i)]; % Matriu canvi de base
                                                   % completa 6x6

    % Càlcul de la submatriu de rigidesa en eixos locals.
    KL=calcul_submatrius_rigidesa(Lmecp(i,:),extrems,L);
    vector_mov_local=[T'*Retorn_moviments(n1*6-5:n1*6);T'*...
        Retorn_moviments(n2*6-5:n2*6)]; % Vector moviments en eixos locals.
    Nepfor(i,:)=calcul_carregues_barra_local(barra_carrega_rep(i,:),...
        barra_carrega_rep_proj(i,:),barra_carrega_punt(i,:),...
        barra_carrega_mom(i,:),TL2G(:,:,i),L,extrems);
    Retorn_esforcos_barres(i,:)=(KL-kriggeom_loc_ev(:,:,i))*...
        vector_mov_local-Nepfor(i,:)'; % Esforços en extrems de cada barra
                                       % en eixos locals.
    Ax(i,j+1)=Retorn_esforcos_barres(i,1);

    % A coninuació evaula la matriu geomètrica local i l'acobla a la matriu
    % geomètrica global. Aquesta matriu geomètrica serà la que el mètode
    % utilitzarà pel pas intermedi (K1).
    kriggeom_cte_ev(:,:,i)=(Ax(i,j+1)/(30*L))*kriggeom_cte_tot(:,:,i);
    kriggeom_loc_ev(:,:,i)=(Ax(i,j+1)/(30*L))*kriggeom_cte(:,:,i);
    K_Geom_Glob_cte=K_Geom_Glob_cte+ensamblatge_matriu_rigidesa...
        (kriggeom_cte_ev(:,:,i),ns,sistema,n1,n2);

end

    % Diferencial dels moviments al pas intermedi. Aquest moviments
    % s'hauràn de sumar al sumatori total de desplaçaments.
    d_vector_moviments_prim(1:nn,1)=(K_Glob(1:nn,1:nn)-...
        K_Geom_Glob_cte(1:nn,1:nn))\(mu*dp(1:nn,1));
    sum_dif_mov_prim=sum_dif_mov+d_vector_moviments_prim;

    while det(K_Glob(1:nn,1:nn)-K_Geom_Glob_cte(1:nn,1:nn))<0 % Si el
                                      % determinant de la matriu geomètrica
                                      % global es superior al determinant
                                      % de la matriu de rigidesa global
        K_Geom_Glob_cte=[]; % Dimensionament de la matriu goemètrica global
                            % buida
        errmet=1; % Valor pertinent per la variable d'error del mètode
        return % Retorna l'execució a la funció calculinestabilitat.m
    end

    % A continuació e calcularàn els esforços en les barres per trobar
    % l'axil de casascuna d'elles i introduir-lo a les matrius geomètriques
    % locals

    % Calcula les reaccions als recolzaments de l'estructura.
    Vector_reaccions=K_Glob*sum_dif_mov_prim-Vector_esforcos;

    % Retorn de reaccions i moviments en l'ordre original dels nusos i
    % eixos globals.
    [Retorn_reaccions,Retorn_moviments]=ordre_inicial_nusos(sistema,...
        Vector_reaccions,sum_dif_mov_prim,n,Ncond);
    Retorn_esforcos_barres=zeros(N,12); % Càlcul dels esforços en extrems
                                        % de barra.
    K_Geom_Glob_cte_prim=zeros(ns);

for i=1:N
    n1=Lnods(i,1); % node inferior
    n2=Lnods(i,2); % node superior
    extrems=Lartyp(i,1)*2+Lartyp(i,2); % rígid-rígid =0
                                       % rígid-articulat =1
                                       % articulat-rígid =2
                                       % articulat-articulat =3
    vector_barra=Ncoord(n2,1:3)-Ncoord(n1,1:3); % Vector de la barra.
    L=norm(vector_barra); % Longitut de la barra.
    T=[TL2G(:,:,i),zeros(3);zeros(3),TL2G(:,:,i)]; % Matriu canvi de base
                                                   % completa 6x6

    % Càlcul de la submatriu de rigidesa en eixos locals.
    KL=calcul_submatrius_rigidesa(Lmecp(i,:),extrems,L);
    vector_mov_local=[T'*Retorn_moviments(n1*6-5:n1*6);T'*...
        Retorn_moviments(n2*6-5:n2*6)]; % Vector moviments en eixos locals.
    Nepfor(i,:)=calcul_carregues_barra_local(barra_carrega_rep(i,:),...
        barra_carrega_rep_proj(i,:),barra_carrega_punt(i,:),...
        barra_carrega_mom(i,:),TL2G(:,:,i),L,extrems);
    Retorn_esforcos_barres(i,:)=(KL-kriggeom_loc_ev(:,:,i))*...
        vector_mov_local-Nepfor(i,:)'; % Esforços en extrems de cada barra
                                       % en eixos locals.
    Ax_prim(i,j+1)=Retorn_esforcos_barres(i,1);

    % A coninuació evaula la matriu geomètrica local i l'acobla a la matriu
    % geomètrica global. Aquesta matriu geomètrica (K2) serà la que el
    % mètode utilitzarà per convinar-la amb l'anterior (K1).
    kriggeom_cte_ev(:,:,i)=(Ax_prim(i,j+1)/(30*L))*kriggeom_cte_tot(:,:,i);
    kriggeom_loc_ev(:,:,i)=(Ax_prim(i,j+1)/(30*L))*kriggeom_cte(:,:,i);
    K_Geom_Glob_cte_prim=K_Geom_Glob_cte_prim+...
        ensamblatge_matriu_rigidesa(kriggeom_cte_ev(:,:,i),ns,sistema,...
        n1,n2);

end

    % Diferencial dels moviments al final del pas. Aquest moviments
    % serviran per obtenir la matriu de rigidesa tangen a l'inici del pas
    % següent..
    d_vector_moviments(1:nn,1)=(K_Glob(1:nn,1:nn)-(K_Geom_Glob_cte...
        (1:nn,1:nn)*alpha1+K_Geom_Glob_cte_prim(1:nn,1:nn)*alpha2))...
        \dp(1:nn,1);
    sum_dif_mov=sum_dif_mov+d_vector_moviments;

    while det(K_Glob(1:nn,1:nn)-K_Geom_Glob_cte(1:nn,1:nn))<0 % Si el
                                      % determinant de la matriu geomètrica
                                      % global es superior al determinant
                                      % de la matriu de rigidesa global
        K_Geom_Glob_cte=[]; % Dimensionament de la matriu goemètrica global
                            % buida
        errmet=1; % Valor pertinent per la variable d'error del mètode
        return % Retorna l'execució a la funció calculinestabilitat.m
    end

end

    % Un cop s'han realitzat tots els increments trobarem la matriu
    % geomètrica global a partir dels desplaçaments finals dels nussos.

    % Calcula les reaccions als recolzaments de l'estructura.
    Vector_reaccions=K_Glob*sum_dif_mov-Vector_esforcos;

    % Retorn de reaccions i moviments en l'ordre original dels nusos i
    % eixos globals.
    [Retorn_reaccions,Retorn_moviments]=ordre_inicial_nusos(sistema,...
        Vector_reaccions,sum_dif_mov,n,Ncond);
    Retorn_esforcos_barres=zeros(N,12); % Càlcul dels esforços en extrems
                                        % de barra.
    K_Geom_Glob_cte=zeros(ns);
for i=1:N
    n1=Lnods(i,1); % node inferior
    n2=Lnods(i,2); % node superior
    extrems=Lartyp(i,1)*2+Lartyp(i,2); % rígid-rígid =0
                                       % rígid-articulat =1
                                       % articulat-rígid =2
                                       % articulat-articulat =3
    vector_barra=Ncoord(n2,1:3)-Ncoord(n1,1:3); % Vector de la barra.
    L=norm(vector_barra); % Longitut de la barra.
    T=[TL2G(:,:,i),zeros(3);zeros(3),TL2G(:,:,i)]; % Matriu canvi de base
                                                   % completa 6x6

    % Càlcul de la submatriu de rigidesa en eixos locals.
    KL=calcul_submatrius_rigidesa(Lmecp(i,:),extrems,L);
    vector_mov_local=[T'*Retorn_moviments(n1*6-5:n1*6);T'*...
        Retorn_moviments(n2*6-5:n2*6)]; % Vector moviments en eixos locals.
    Nepfor(i,:)=calcul_carregues_barra_local(barra_carrega_rep(i,:),...
        barra_carrega_rep_proj(i,:),barra_carrega_punt(i,:),...
        barra_carrega_mom(i,:),TL2G(:,:,i),L,extrems);
    Retorn_esforcos_barres(i,:)=(KL-kriggeom_loc_ev(:,:,i))*...
        vector_mov_local-Nepfor(i,:)'; % Esforços en extrems de cada barra
                                       % en eixos locals.
    Ax_fin(i,1)=Retorn_esforcos_barres(i,1);

    % A coninuació evaula la matriu geomètrica local i l'acobla a la matriu
    % geomètrica global. Aquesta matriu geomètrica serà la que la funció
    % retornarà.
    kriggeom_cte_ev(:,:,i)=(Ax_fin(i,1)/(30*L))*kriggeom_cte_tot(:,:,i);
    kriggeom_loc_ev(:,:,i)=(Ax_fin(i,1)/(30*L))*kriggeom_cte(:,:,i);
    K_Geom_Glob_cte=K_Geom_Glob_cte+ensamblatge_matriu_rigidesa...
        (kriggeom_cte_ev(:,:,i),ns,sistema,n1,n2);

end

end

FUNCIÓ PER REORDENAR ELS MODES DE VINCLAMENT

% Funció que ordena els vectors del mode de vinclament perquè aquest
% tinguin l'ordre dels nussos

function v_reord=reordenacio_vec(v,sistema,n,Ncond)

v_reord=zeros(n*6,1);

for i=1:n
    for j=1:6
        k=sistema(6*i+j-6);
        if k
            v_reord(6*i+j-6)=v(k);
        end
    end
    if Ncond(i,1)==1 % Si és rodet inclinat torna desplaçaments del sistema
                     % d'eixos locals del rodet al global
        Tr=calcul_canvi_base(Ncond(i,2:4)); % Calcula matriu canvi de base
                                            % del rodet.
        T=[Tr,zeros(3);zeros(3),Tr];
        v_reord(6*i-5:6*i)=T*v_reord(6*i-5:6*i);
    end
end
end