kill(all);

## Auxiliary definitions

I define two auxiliary functions for expressing the moments, where linear
on a branch

f( clarabella ) := 1 - clarabella ;

g( orazio ) := orazio ;

## Bending and torsional moments

### AB beam segment

Moment components at nodes

M1Ap :  fullratsimp( 0*P + 0*YA + 0*XA     - 1*PsiA + 0*F );
M2Ap :  fullratsimp( 0*H + 0*ZA + 1*ThetaA + 0*PhiA + 0*I );
MtAp :  fullratsimp( 0*H + 0*ZA + 0*ThetaA + 1*PhiA + 0*I );

M1Bm :  fullratsimp( 0*P + 0*YA + h*XA     - 1*PsiA + 0*F );
M2Bm :  fullratsimp( 0*H + h*ZA + 1*ThetaA + 0*PhiA + 0*I );
MtBm :  fullratsimp( 0*H + 0*ZA + 0*ThetaA + 1*PhiA + 0*I );

Moment components along the beam segment

M1 : fullratsimp( M1Ap*f(s/h) + M1Bm*g(s/h));
M2 : fullratsimp( M2Ap*f(s/h) + M2Bm*g(s/h));
Mt : fullratsimp( MtAp*f(s/h) + MtBm*g(s/h));

Strain energy calculation

U : integrate(
    M1^2/2/E/J + M2^2/2/E/J +Mt^2/2/G/Kt,
    s,0,h
)$
U : expand( U );

### BC beam segment

Moment components at nodes

M1Bp : M1Bm;
M2Bp : fullratsimp( 0*H + h*ZA + 1*ThetaA + 0*PhiA + 0*I);
MtBp : fullratsimp( 0*H + 0*ZA + 0*ThetaA + 1*PhiA + 0*I);

M1Cm : fullratsimp( 0*P + 0*YA + a*XA -PsiA - (a-h)*F);
M2Cm : fullratsimp( 0*H + a*ZA + 1*ThetaA + 0*PhiA + 0*I);
MtCm : fullratsimp( 0*H + 0*ZA + 0*ThetaA + 1*PhiA + 0*I);

Moment components along the beam segment

M1 : fullratsimp( M1Bp*f(s/(a-h)) + M1Cm*g(s/(a-h)));
M2 : fullratsimp( M2Bp*f(s/(a-h)) + M2Cm*g(s/(a-h)));
Mt : fullratsimp( MtBp*f(s/(a-h)) + MtCm*g(s/(a-h)));

Strain energy calculation AND ACCUMULATION

U : 
U + integrate(
    M1^2/2/E/J + M2^2/2/E/J +Mt^2/2/G/Kt,
    s,0,a-h
)$
U : expand( U );

### CD beam segment

Moment components at nodes

M1Cp : M1Cm ;
M2Cp : fullratsimp( 0*H + 0*ZA + 0*ThetaA - 1*PhiA + 0*I);
MtCp : fullratsimp( 0*H + a*ZA + 1*ThetaA + 0*PhiA + 0*I);

M1Dm : fullratsimp( 0*P - b*YA + a*XA -PsiA - (a-h)*F);
M2Dm : fullratsimp(-b*H + b*ZA + 0*ThetaA - 1*PhiA + 0*I);
MtDm : fullratsimp( 0*H + a*ZA + 1*ThetaA + 0*PhiA + 0*I);

Moment components along the beam segment

M1 : fullratsimp( M1Cp*f(s/b) + M1Dm*g(s/b));
M2 : fullratsimp( M2Cp*f(s/b) + M2Dm*g(s/b));
Mt : fullratsimp( MtCp*f(s/b) + MtDm*g(s/b));

Strain energy calculation AND ACCUMULATION

U : 
U + integrate(
    M1^2/2/E/J + M2^2/2/E/J +Mt^2/2/G/Kt,
    s,0,b
)$
U : expand( U );

### DE beam segment

Moment components at nodes

M1Dp : M1Dm;
M2Dp : fullratsimp( 0*H - a*ZA - 1*ThetaA + 0*PhiA + 0*I);
MtDp : fullratsimp(-b*H + b*ZA + 0*ThetaA - 1*PhiA + 0*I);

M1Em : fullratsimp( a*P - b*YA + 0*XA -PsiA + h*F);
M2Em : fullratsimp(-a*H + 0*ZA - 1*ThetaA + 0*PhiA - a*I);
MtEm : fullratsimp(-b*H + b*ZA + 0*ThetaA - 1*PhiA + 0*I);

Moment components along the beam segment

M1 : fullratsimp( M1Dp*f(s/a) + M1Em*g(s/a));
M2 : fullratsimp( M2Dp*f(s/a) + M2Em*g(s/a));
Mt : fullratsimp( MtDp*f(s/a) + MtEm*g(s/a));

Strain energy calculation AND ACCUMULATION

U : 
U + integrate(
    M1^2/2/E/J + M2^2/2/E/J +Mt^2/2/G/Kt,
    s,0,a
)$
U : expand( U );

## some further checks on the internal action components at nodes

double check all the nodal values!

Fig 2 graphs

l:a-h$

 [ M1A, M1B, M1C, M1D, M1E] :
 [  0 ,  0 ,  0 ,  0 ,  a ]*P
+[  0 ,  0 ,  0 , -b , -b ]*YA
+[  0 ,  h ,  a ,  a ,  0 ]*XA
+[ -1 , -1 , -1 , -1 , -1 ]*PsiA
+[  0 ,  0 , -l , -l , +h ]*F;

comparison with the previously defined M1

fullratsimp(
    [M1A,M1B,M1C,M1D]-[M1Ap,M1Bp,M1Cp,M1Dp]
);

fullratsimp(
    [M1B,M1C,M1D,M1E]-[M1Bm,M1Cm,M1Dm,M1Em]
);

a null difference, WOW!!! 

Figure 3 graphs, bending moment component 2

fullratsimp(
  - [M2Ap,M2Bm,M2Bp,M2Cm,M2Cp,M2Dm,M2Dp,M2Em]
  + [   0,   0,   0,   0,   0,  -b,   0,  -a]*H
  + [   0,  +h,  +h,  +a,   0,  +b,  -a,   0]*ZA
  + [  +1,  +1,  +1,  +1,   0,   0,  -1,  -1]*ThetaA
  + [   0,   0,   0,   0,  -1,  -1,   0,   0]*PhiA
  + [   0,   0,   0,   0,   0,   0,   0,  -a]*I
);

torsional moment

fullratsimp(
  - [MtAp,MtBm,MtBp,MtCm,MtCp,MtDm,MtDp,MtEm]
  + [   0,   0,   0,   0,   0,   0,  -b,  -b]*H
  + [   0,   0,   0,   0,  +a,  +a,  +b,  +b]*ZA
  + [   0,   0,   0,   0,  +1,  +1,   0,   0]*ThetaA
  + [  +1,  +1,  +1,  +1,   0,   0,  -1,  -1]*PhiA
  + [   0,   0,   0,   0,   0,   0,   0,   0]*I
);

## Solution of the statically indederminated problem

Reaction forces at A calculation

eqns : [
    diff( U , XA ) = 0,
    diff( U , YA ) = 0,
    diff( U , ZA ) = 0,
    diff( U , ThetaA ) = 0,
    diff( U , PhiA ) = 0,
    diff( U , PsiA )=0
] ;

a way to assign the solution to the unknowns

sols : linsolve( eqns , [ XA , YA, ZA , ThetaA, PhiA , PsiA ] ) ;

[ XA , YA, ZA , ThetaA, PhiA , PsiA ] : 
ev(
    [ XA , YA, ZA , ThetaA, PhiA , PsiA ] 
    , sols
)$

another way to assign the solution to the unknown, relying on the 
"globalsolve" control variable that, if set to "true", switches the linsolve
function behavior from "not assigning" to "assigning".

linsolve( eqns , [ XA , YA, ZA , ThetaA, PhiA , PsiA ] ) , globalsolve=true;

## Evaluation of the deflections at the monitored d.o.f.s

I redefine the strain energy as a function of F,P,H,I alone

U : fullratsimp( ev( U ) ) $

Deflection calculation at the B point, in the negative x direction

utB

diff( U , F ) $
% , F=0, I=0$
fullratsimp(%)$
utB: % ;

Deflection calculation at the C point, in the negative x direction

utC

diff( U , P )$
% , F=0, I=0$
fullratsimp(%)$
utC: %;

check for the limit for point B approaching point C

limit( defB , h, a);

ok, it behaves as expected.

out-of-plane deflections at C and D, named
wtC and wtD

diff( U , H )$
% , F=0, I=0$
fullratsimp( % )$
wtC : %;

diff( U , I )$
% , F=0, I=0$
fullratsimp( % )$
wtD : %;

## Back to the evaluation of the constraint reactions

evaluation of the reaction forces and moments at A

ev([ XA , YA, ZA , ThetaA, PhiA , PsiA ] , [ F=0,I=0 ])$
fullratsimp(%)$
[ XA , YA, ZA , ThetaA, PhiA , PsiA ] : %;

evaluation of the reaction forces and moments at E

eqtx : -P-F+XA+XE=0$
eqty :      YA+YE=0$
eqtz : -H-I+ZA+ZE=0$

eqrEx : +a*H + a*I         + ThetaA + ThetaE = 0$
eqrEy : +b*H - b*ZA        + PhiA   + PhiE   = 0$
eqrEz : -a*P - h*F  + YA*b + PsiA   + PsiE   = 0$

linsolve( 
    ev(
        [ eqtx,eqty,eqtz,eqrEx,eqrEy,eqrEz],
        F=0,I=0
    ) , 
    [XE,YE,ZE,ThetaE,PhiE,PsiE] 
), globalsolve=true;

## A numerical test case, and a comparison with FEA analysis

reference sizing for the parametrically defined dimensions

dim:[
    E=210000,
    nu=0.3,
    a=750,
    h=600,
    b=450,
    t=2,
    d=38,
    P=1000,
    H=850
]$
dim:append(
    dim,
    ev(
        [J=%pi*d^3*t/8,Kt=%pi*d^3*t/4,G=E/2/(1+nu)],
        dim    
    )
);

monitored deflections, 

[utB,utC,wtC,wtD], dim$
%,numer;

comparison data from FEA, and the relative deviation

(%-[2.07931,2.471,7.97458,5.23294])/max(abs(%));

reaction forces at A

[XA,YA,ZA],dim$
%,numer;

comparison data from FEA, and the relative deviation

(%-[500,757.575,616.608])/max(abs(%));

reaction moments at A

[ThetaA,PhiA,PsiA],dim$
%,numer;

comparison data from FEA, and the relative deviation

(%-[-410763,-52513.3,204546])/max(abs(%));

reaction forces at E

[XE,YE,ZE],dim$
%,numer;

comparison data from FEA, and the relative deviation

(%-[500,-757.575,233.392])/max(abs(%));

reaction moments at E

[ThetaE,PhiE,PsiE],dim$
%,numer;

comparison data from FEA, and the relative deviation

(%-[-226737,-52513.3,204546])/max(abs(%));

The FEA doublecheck pointed out a forgotten "b*H" contribution to the "eqrEy"
equilibrium equation. I added it and now it is OK.

More attention is required in casting those equations, since a doublecheck
is not always feasible!!