Definition of the Hermitian polynomilas

(Subscript[n, 1])[ξ_]:=(1/4) (2-3 ξ+(ξ^3));
(Subscript[n, 2])[ξ_]:= (L/8) (1-ξ-(ξ^2)+(ξ^3));
(Subscript[n, 3])[ξ_]:=(1/4) (2+3 ξ-(ξ^3));
(Subscript[n, 4])[ξ_]:=-(L/8) (1+ξ-(ξ^2)-(ξ^3));

Vector of interpolation function

Nbend=Table[{(Subscript[n, i])[ξ]},{i,4}];
Bbend=(4/(L^2)) D[Nbend,{ξ,2}];
Bbendt= Transpose[Bbend];

Construction of bending stiffness matrix

Kebend=Iz Ex L/2 Integrate[Bbend.Bbendt,{ξ,-1,1}];

Definition of interpolation functions for tension-compression

(Subscript[n, 5])[ξ_]:=(1/2) (1-ξ);
(Subscript[n, 6])[ξ_]:=(1/2) (1+ξ);

Vector of interpolation functions, tension-compression

Ntens=Table[{(Subscript[n, i])[ξ]},{i,5,6}];
Ntenst=Transpose[Ntens];
Btens=(2/L) D[Ntens,ξ];
Btenst=Transpose[Btens];

Construction od extensional stiffness matrix

Ketens=A Ex L/2 Integrate[Btens.Btenst,{ξ,-1,1}];

Construction of the stiffness matrix, general beam element (bending+tension/compression)

Ke=Table[0,{6},{6}];
Ke[[{1,4},{1,4}]]=Ketens;
Ke[[{2,3,5,6},{2,3,5,6}]]=Kebend;

Construction of the element stiffness matrices in the local (element) coordinate system (CS)

Do[(Subscript[Keloc, i])=Ke/.L->(Subscript[L, i]),{i,2}]

Definition of the transformation matrices

Do[{
(Subscript[T, i])={{(Subscript[c, i]),-(Subscript[s, i]),0,0,0,0},{(Subscript[s, i]),(Subscript[c, i]),0,0,0,0},{0,0,1,0,0,0},{0,0,0,(Subscript[c, i]),-(Subscript[s, i]),0},{0,0,0,(Subscript[s, i]),(Subscript[c, i]),0},{0,0,0,0,0,1}};
(Subscript[Tt, i])=Transpose[(Subscript[T, i])]
},{i,2}]

Check the orthogonality of the transformation matrices

Do[Print[MatrixForm[(Subscript[T, i]).(Subscript[Tt, i])]],{i,2}]

Transformation of the stiffness matrices to the global CS

Do[(Subscript[KE, i])=(Subscript[T, i]).(Subscript[Keloc, i]).(Subscript[Tt, i]),{i,2}]

Construction of the structural level matrices by expanding the element matrices

K=Table[0,{9},{9}];
K[[1;;6,1;;6]]=K[[1;;6,1;;6]]+(Subscript[KE, 1]);
K[[4;;9,4;;9]]=K[[4;;9,4;;9]]+(Subscript[KE, 2]);

Definition of the local line load function and the load vector for element 1

P[ξ_]:=-((((Subscript[p, 1])+(Subscript[p, 2])))/2)+((((Subscript[p, 1])-(Subscript[p, 2])))/2) ξ;
F1ploc=(L/2) Integrate[{0,P[ξ] (Subscript[n, 1])[ξ], P[ξ] (Subscript[n, 2])[ξ],0,P[ξ] (Subscript[n, 3])[ξ],P[ξ] (Subscript[n, 4])[ξ]},{ξ,-1,1}]/.L->(Subscript[L, 1]);

Transformation of the load vector to the global CS

F1p=(Subscript[T, 1]).F1ploc;

Construction of the structural force and nodal displacement vector

Ay=.;M1=.;Bx=.;By=.;u1=.;u2=.;v1=.;v2=.;φ1=.;φ2=.;
F={0,Ay,M1,0,-Fc,0,Bx,By,0}+Flatten[{F1p,{0,0,0}}];
U={u1,0,0,u2,v2,φ2,0,0,φ3};

Definition of the geometrical and material properties (units are in SI)

a=1;b=1.2;cpar=1.1;(Subscript[L, 1])=(Sqrt[(a^2)+(cpar^2)]);(Subscript[L, 2])=b;Iz=21 (10^(-8));Ex=200 (10^9);A=150 (10^(-6));(Subscript[p, 1])=600;(Subscript[p, 2])=800;Fc=10000;

Definition of the nodal coordinates, plot of the structure

(Subscript[x, 1])=0;(Subscript[y, 1])=0;(Subscript[x, 2])=a;(Subscript[y, 2])=cpar;(Subscript[x, 3])=a+b;(Subscript[y, 3])=cpar;
Nodes=Table[{(Subscript[x, i]),(Subscript[y, i])},{i,3}];
Graphics[{{Line[Nodes]}}]

Definition of the transformation matrix components

Do[(Subscript[c, i])=(((Subscript[x, i+1])-(Subscript[x, i]))/(Subscript[L, i]));(Subscript[s, i])=(((Subscript[y, i+1])-(Subscript[y, i]))/(Subscript[L, i])),{i,2}]

Construction and solution of the structural equation system, nodal displacements

adof={1,4,5,6,9};
Kcon=K[[adof,adof]];
Fcon=F[[adof]];
Ucon=LinearSolve[Kcon,Fcon]

Calculation and assignment of the reaction

U[[adof]]=Ucon;
fixdof=Complement[Table[i,{i,9}],adof];
sol=Solve[(K.U)[[fixdof]]==F[[fixdof]],{Ay,M1,Bx,By}]

Calculationof the element displacement vectors in the global and local CSs

(Subscript[ue, 1])=U[[1;;6]];
(Subscript[ue, 2])=U[[4;;9]];
Do[(Subscript[ueloc, i])=(Subscript[Tt, i]).(Subscript[ue, i]),{i,2}]

Calculation of the internal force vectors (= beam diagram data), local CS

(Subscript[Feintloc, 1])=(Subscript[Keloc, 1]).(Subscript[ueloc, 1])-F1ploc
(Subscript[Feintloc, 2])=(Subscript[Keloc, 2]).(Subscript[ueloc, 2])