"Parameters:";
Par={
a->0.6,
b->0.4,
(Subscript[A, 1])->100 (10^(-6)),
(Subscript[A, 2])->120 (10^(-6)),
(Subscript[A, 3])->90 (10^(-6)),
(Subscript[Ex, 1])->210 (10^9),
(Subscript[Ex, 2])->180 (10^9),
(Subscript[Ex, 3])->195 (10^9),
F1->100 (10^3)
};
"Topology of elements:";
(Subscript[Tp, 1])={1,2};
(Subscript[Tp, 2])={2,3};
(Subscript[Tp, 3])={1,3};
"Nodal coordinates:";
(Subscript[Nod, 1])={0,0};
(Subscript[Nod, 2])={a,0};
(Subscript[Nod, 3])={a,b};
"Constraints on the nodes:";
Con={0,1,1,1,0,0};
"Load vector:";
F={0,Ay,Bx,By,F1,0};
"Number of free DoFs:";
Nu=0;
Do[If[Con[[i]]<1,Nu=Nu+1],{i,Dimensions[Con][[1]]}]
"Storing the free DoFs:";
DoF=ConstantArray[0,Nu];
ix=0;
Do[If[Con[[j]]<1,ix=ix+1;DoF[[ix]]=j],{j,Dimensions[Con][[1]]}]
"Element lengths, and deflection:";
Do[{
i1=(Subscript[Tp, i])[[1]];
i2=(Subscript[Tp, i])[[2]];
(Subscript[L, i])=(Sqrt[((((Subscript[Nod, i1])[[1]]-(Subscript[Nod, i2])[[1]]))^2)+((((Subscript[Nod, i1])[[2]]-(Subscript[Nod, i2])[[2]]))^2)]);
(Subscript[s, i])=((Subscript[Nod, i2])[[2]]-(Subscript[Nod, i1])[[2]])/(Subscript[L, i]);
(Subscript[c, i])=((Subscript[Nod, i2])[[1]]-(Subscript[Nod, i1])[[1]])/(Subscript[L, i]);
},
{i,3}]
"Empty global stiffness matrix:";
Km=Table[0,{i,6},{j,6}];
"Filling up of the global stiffness matrix:";
Do[
{(Subscript[Ke, i])=(Subscript[A, i]) (Subscript[Ex, i])/(Subscript[L, i]) {{1,-1},{-1,1}},{i,3};
(Subscript[T, i])={{(Subscript[c, i]),0},{(Subscript[s, i]),0},{0,(Subscript[c, i])},{0,(Subscript[s, i])}};
(Subscript[Keg, i])=(Subscript[T, i]).(Subscript[Ke, i]).Transpose[(Subscript[T, i])];
Ind=Table[2 (Subscript[Tp, i])[[IntegerPart[j/2-1/2]+1]]-1+If[OddQ[j]==True,0,1],{j,4}];
Do[Km[[Ind[[j]],Ind[[k]]]]=Km[[Ind[[j]],Ind[[k]]]]+(Subscript[Keg, i])[[j,k]],{j,1,4},{k,1,4}]
},
{i,3}
];
"Solution of the condensed equation";
Fc=F[[DoF]];
Kc=Km[[DoF,DoF]];
Uc=Inverse[Kc].Fc;
"Global displacement vector";
U=ConstantArray[0,6];
ix=0;
Do[If[Con[[i]]==0,ix=ix+1;U[[i]]=Uc[[ix]],0],{i,Dimensions[Con][[1]]}]
"Displacements:"
FullSimplify[U]
U/.Par
"Forces:"
F=FullSimplify[Km.U]
F/.Par