(*definiciones*)

(*matriz de paso: globales -> locales*)
T[x_,y_]:={{x,y,0,0,0,0},{-y,x,0,0,0,0},{0,0,1,0,0,0},
		{0,0,0,x,y,0},{0,0,0,-y,x,0},{0,0,0,0,0,1}};
MatrixForm[T[(Subscript[λ, X]),(Subscript[λ, Y])]]

(*matriz de rigidez elemental*)
k[A_,Iz_,L_,vE_]:={{A vE/L,0,0,-A vE/L,0,0},
		{0,12 vE Iz/L^3,6 vE Iz/L^2,0,-12 vE Iz/L^3,6 vE Iz/L^2},
		{0,6 vE Iz/L^2,4 vE Iz/L,0,-6 vE Iz/L^2,2 vE Iz/L},
		{-A vE/L,0,0,A vE/L,0,0},
		{0,-12 vE Iz/L^3,-6 vE Iz/L^2,0,12 vE Iz/L^3,-6 vE Iz/L^2},
		{0,6 vE Iz/L^2,2 vE Iz/L,0,-6 vE Iz/L^2,4 vE Iz/L}};
MatrixForm[k[A,Iz,L,vE]]

(* HEB 300 *)
A1=149.1 10^-4;
I1=25166 10^-8;

(* IPE 450 *)
A2=98.8 10^-4;
I2=33740 10^-8;

(* discretizacion *)

nd=5;(*num de nudos*)
nb=4;(*num de barras*)
(*datos*)
A={A1,A2,A2,A1};(*area*)
Iz={I1,I2,I2,I1}; (*momento de inercia*)
h=.; (*canto de la seccion*)
vE=2.1 10^11;(*modulo de Young*)
α=.
L=.;(*longitud*)

(*coord de los nudos*)
xy=Table[0,{nd},{2}];
Part[xy,2,2]=7;
Part[xy,3,1]=12.5;
Part[xy,3,2]=9.5;
Part[xy,4,1]=25;
Part[xy,4,2]=7;
Part[xy,5,1]=25;
MatrixForm[xy]
azul={PlotStyle->RGBColor[0, 0, 1]};
ops={azul,PlotJoined->True,AspectRatio->Automatic};
ListPlot[xy,ops]

(*definicion de elementos/barras*)
barras=Table[0,{nb},{2}];
Part[barras,1,1]=1;
Part[barras,1,2]=2;
Part[barras,2,1]=2;
Part[barras,2,2]=3;
Part[barras,3,1]=3;
Part[barras,3,2]=4;
Part[barras,4,1]=5;
Part[barras,4,2]=4;
MatrixForm[barras]

(*matriz de rigidez elemental (globales)*)
Ke=Table[0.0,{nb},{1},{6},{6}];


For[b=1,b<=nb,b++,{
		Print["Elemento = ",b],
	i=barras[[b]][[1]],
	j=barras[[b]][[2]],
	L=√((xy[[j]][[1]]-xy[[i]][[1]])^2+(xy[[j]][[2]]-xy[[i]][[2]])^2),
	(Subscript[λ, X])=(xy[[j]][[1]]-xy[[i]][[1]])/L,
	Print["	\!\(\*SubscriptBox[\(λ\), \(\(X\)\(\\ \)\)]\)= ",N[(Subscript[λ, X]),3]],
	(Subscript[λ, Y])=(xy[[j]][[2]]-xy[[i]][[2]])/L,
	Print["	\!\(\*SubscriptBox[\(λ\), \(\(Y\)\(\\ \)\)]\)= ",N[(Subscript[λ, Y]),3]],
	Part[Ke,b,1]= Transpose[T[(Subscript[λ, X]),(Subscript[λ, Y])]].k[A[[b]],Iz[[b]],L,vE].T[(Subscript[λ, X]),(Subscript[λ, Y])],
	Print[MatrixForm[Ke[[b]][[1]]]]
	}
]

(*matriz de rigidez de la estructura*)
s=3; (*gdl/nudo, tipo de problema*)
ngdl=s nd;(*grados de libertad*)
K=Table[0.0,{ngdl},{ngdl}];

For[b=1,b<=nb,b++,{
	(*Print["Elemento = ",b];*)
	i=barras[[b]][[1]];
	j=barras[[b]][[2]];
	(*submatriz 11*)
	For[f=1,f<=s,f++,(*filas*)
				For[c=1,c<=s,c++,(*columnas*)
						Part[K,s(i-1)+f,s(i-1)+c]+=Part[Ke,b,1,f,c]
					]
		]

		(*submatriz 12*)
		For[f=1,f<=s,f++,
				For[c=1,c<=s,c++,Part[K,s(i-1)+f,s(j-1)+c]+=Part[Ke,b,1,f,s+c]
					]
		]
		
		(*submatriz 21*)
		For[f=1,f<=s,f++,
				For[c=1,c<=s,c++,Part[K,s(j-1)+f,s(i-1)+c]+=Part[Ke,b,1,s+f,c]
					]
		]
	
		(*submatriz 22*)
	For[f=1,f<=s,f++,
				For[c=1,c<=s,c++,Part[K,s(j-1)+f,s(j-1)+c]+=Part[Ke,b,1,s+f,s+c]
					]
		]
	
		}
]
	MatrixForm[N[K,7]]

(*condiciones de contorno*)
u=Array[Λ,{ngdl}]; (*desplazamientos y giros*)

nodo=1;
dir=1;(*direccion*)
u[[s nodo-1(s-dir)]]=0.0;
dir=2;(*direccion*)
u[[s nodo-1(s-dir)]]=0.0;
dir=3;(*direccion*)
u[[s nodo-1(s-dir)]]=0.0;

nodo=5;
dir=1;(*direccion*)
u[[s nodo-1(s-dir)]]=0.0;
dir=2;(*direccion*)
u[[s nodo-1(s-dir)]]=0.0;
dir=3;(*direccion*)
u[[s nodo-1(s-dir)]]=0.0;

F=Array[W,{ngdl}]; (*fuerzas y momentos en los nudos*)

P0=10000;
nodo=2; (*nodo*)
dir=1;(*direccion*)
F[[s nodo-1(s-dir)]]=P0;
dir=2;(*direccion*)
F[[s nodo-1(s-dir)]]=-P0;
dir=3;(*direccion*)
F[[s nodo-1(s-dir)]]=0.0;

nodo=3; (*nodo*)
dir=1;(*direccion*)
F[[s nodo-1(s-dir)]]=0.0;
dir=2;(*direccion*)
F[[s nodo-1(s-dir)]]=-P0;
dir=3;(*direccion*)
F[[s nodo-1(s-dir)]]=0.0;

nodo=4; (*nodo*)
dir=1;(*direccion*)
F[[s nodo-1(s-dir)]]=P0;
dir=2;(*direccion*)
F[[s nodo-1(s-dir)]]=-P0;
dir=3;(*direccion*)
F[[s nodo-1(s-dir)]]=0.0;

MatrixForm[u] (*desplazamientos*)

MatrixForm[F](*fuerzas*)

sol=Solve[F==K.u][[1]];
u=u/.sol
F=F/.sol

(*esfuerzos: axil,cortante y flector*)
esf=Table[0.0,{nb},{1},{3}];


For[b=1,b<=nb,b++,{
		Print["Elemento = ",b];
	i=barras[[b]][[1]];
	j=barras[[b]][[2]];
	L=√((xy[[j]][[1]]-xy[[i]][[1]])^2+(xy[[j]][[2]]-xy[[i]][[2]])^2);
	(Subscript[λ, X])=(xy[[j]][[1]]-xy[[i]][[1]])/L;
	(Subscript[λ, Y])=(xy[[j]][[2]]-xy[[i]][[2]])/L;
	gi=Table[s (i-1)+l,{l,s}];
	gj=Table[s (j-1)+l,{l,s}];
	gdl=Join[gi,gj];
	Δ=u[[gdl]];
 Part[esf,b,1]=k[A[[b]],Iz[[b]],L,vE].T[(Subscript[λ, X]),(Subscript[λ, Y])].Δ;
 Part[esf,b,1]=Simplify[ Part[esf,b,1]];
	Print[MatrixForm[N[Part[esf,b,1],7]]]
		}
]

(*Pandeo Global: funciones de estabilidad*)

(*matriz de rigidez elemental en funcion del Axil*)
k[A_,Iz_,L_,vE_,x_]:=Module[{Nx=x,a,Δ,C1,C2,C3,C4},
          a=Sqrt[Abs[Nx]/(vE Iz)];

Δ=2(1-Cos[a L])-a L Sin[a L];
C1[Nx_]:=Piecewise[{{a^2 L^2 (1-Cos[a L])/Δ,Nx< 0},{a^2 L^2/(a L/Tanh[a L/2]-2),Nx>0},{6,Nx==0}}];
C2[Nx_]:=Piecewise[{{(a L Sin[a L]-a^2 L^2 Cos[a L])/Δ,Nx< 0},{a L(a L Cosh[a L]-Sinh[a L])/(2(1-Cosh[a L])+a L Sinh[a L]),Nx>0},{4,Nx==0}}];
C3[Nx_]:=Piecewise[{{(a^2 L^2-a L Sin[a L])/Δ,Nx< 0},{a L(Sinh[a L]-a L)/(2(1-Cosh[a L])+a L Sinh[a L]),Nx>0},{2,Nx==0}}];
C4[Nx_]:=Piecewise[{{a^3 L^3 Sin[a L]/Δ,Nx< 0},{a^3 L^3/(a L-2 Tanh[a L/2]),Nx>0},{12,Nx==0}}];

{{A vE/L,0,0,-A vE/L,0,0},
          		{0,C4[Nx] vE Iz/L^3,C1[Nx] vE Iz/L^2,0,-C4[Nx] vE Iz/L^3,C1[Nx] vE Iz/L^2},
          		{0,C1[Nx] vE Iz/L^2,C2[Nx] vE Iz/L,0,-C1[Nx] vE Iz/L^2,C3[Nx] vE Iz/L},
          		{-A vE/L,0,0,A vE/L,0,0},
          		{0,-C4[Nx] vE Iz/L^3,-C1[Nx] vE Iz/L^2,0,C4[Nx] vE Iz/L^3,-C1[Nx] vE Iz/L^2},
          		{0,C1[Nx] vE Iz/L^2,C3[Nx] vE Iz/L,0,-C1[Nx] vE Iz/L^2,C2[Nx] vE Iz/L}} 

      ];

(* axiles *)
Nx=Table[0,nb];
For[b=1,b<=nb,b++,
Nx[[b]]=esf[[b]][[1]][[4]];
];
Nx

(*matriz de rigidez elemental (globales)*)
Ke=Table[0.0,{nb},{1},{6},{6}];


For[b=1,b<=nb,b++,{
		Print["Elemento = ",b],
	i=barras[[b]][[1]],
	j=barras[[b]][[2]],
	L=√((xy[[j]][[1]]-xy[[i]][[1]])^2+(xy[[j]][[2]]-xy[[i]][[2]])^2),
	(Subscript[λ, X])=(xy[[j]][[1]]-xy[[i]][[1]])/L,
	Print["	\!\(\*SubscriptBox[\(λ\), \(\(X\)\(\\ \)\)]\)= ",N[(Subscript[λ, X]),3]],
	(Subscript[λ, Y])=(xy[[j]][[2]]-xy[[i]][[2]])/L,
	Print["	\!\(\*SubscriptBox[\(λ\), \(\(Y\)\(\\ \)\)]\)= ",N[(Subscript[λ, Y]),3]],
	Part[Ke,b,1]= Transpose[T[(Subscript[λ, X]),(Subscript[λ, Y])]].Simplify[k[A[[b]],Iz[[b]],L,vE,λ Nx[[b]]],λ>0].T[(Subscript[λ, X]),(Subscript[λ, Y])],
	Print[MatrixForm[Ke[[b]][[1]]]]
	}
];

(*matriz de rigidez de la estructura*)
s=3; (*gdl/nudo, tipo de problema*)
ngdl=s nd;(*grados de libertad*)
K=Table[0.0,{ngdl},{ngdl}];

For[b=1,b<=nb,b++,{
	(*Print["Elemento = ",b];*)
	i=barras[[b]][[1]];
	j=barras[[b]][[2]];
	(*submatriz 11*)
	For[f=1,f<=s,f++,(*filas*)
				For[c=1,c<=s,c++,(*columnas*)
						Part[K,s(i-1)+f,s(i-1)+c]+=Part[Ke,b,1,f,c]
					]
		]

		(*submatriz 12*)
		For[f=1,f<=s,f++,
				For[c=1,c<=s,c++,Part[K,s(i-1)+f,s(j-1)+c]+=Part[Ke,b,1,f,s+c]
					]
		]
		
		(*submatriz 21*)
		For[f=1,f<=s,f++,
				For[c=1,c<=s,c++,Part[K,s(j-1)+f,s(i-1)+c]+=Part[Ke,b,1,s+f,c]
					]
		]
	
		(*submatriz 22*)
	For[f=1,f<=s,f++,
				For[c=1,c<=s,c++,Part[K,s(j-1)+f,s(j-1)+c]+=Part[Ke,b,1,s+f,s+c]
					]
		]
	
		}
]
	MatrixForm[N[K,7]]

(*resolucion*)
orden=0;
For[i=1,i<=ngdl,i++,If[u[[i]]==0.0,orden++]]

cono=Table[0,{orden}]; (*gdl conocidos*)
For[i=1;a=1,i<=ngdl,i++,If[u[[i]]==0.0,{Part[cono,a]=i,a++}]]

kk=Table[i,{i,1,ngdl}];
incog=Complement[kk, cono];(*gdl incognita*)
cargas=F[[incog]];
MatrixForm[cargas];
M=K[[incog,incog]];
MatrixForm[N[M,3]]

kk=Det[M];

λcr=λ/.FindRoot[kk,{λ,30}]//Quiet

(* beta - barra 4 *)

b=4;
i=barras[[b]][[1]];
j=barras[[b]][[2]];
L=√((xy[[j]][[1]]-xy[[i]][[1]])^2+(xy[[j]][[2]]-xy[[i]][[2]])^2);
Ncr=-λcr Nx[[b]];
β2=β/.Solve[(Pi/(β L))^2 vE Iz[[b]]==Ncr,β][[2]]//Quiet