kill(all);

## Auxiliary definitions

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

f( normalizedabscissa ) := 1 - normalizedabscissa ;

g( normalizedabscissa ) := normalizedabscissa ;

## Bending and torsional moments

### AB beam segment

Moment components at nodes

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

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

Moment components along the beam segment

M1 : ratsimp( M1Ap*f(s/h) + M1Bm*g(s/h));
M2 : ratsimp( M2Ap*f(s/h) + M2Bm*g(s/h));
Mt : ratsimp( 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, /*ds*/
    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 );

please doublecheck that the previous expression is NOT a function of any reaction force which may be expressed as a function of 
- the external loads
- the parametrically defined reaction forces
based on equilibrium alone, i.e. it must NOT be a function of the reactions of the principal structure residual constraints.

## 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
);

## Application of the Castigliano theorem to derive the attentioned displacements

    uA:diff(U,    XA,1)$
    vA:diff(U,    YA,1)$
    wA:diff(U,    ZA,1)$
thetaA:diff(U,ThetaA,1)$
  phiA:diff(U,  PhiA,1)$
  psiA:diff(U,  PsiA,1)$

utB:diff(U,F)$
utC:diff(U,P)$
wtC:diff(U,H)$
wtD:diff(U,I)$

I enforce the kinematic contraints at A

direct assignment:
linsolve( 
    [ uA=0 , vA=0 , wA=0, thetaA=0, phiA=0, psiA=0 ] , 
    [ XA   , YA   , ZA  , ThetaA  , PhiA  , PsiA   ] 
) , globalsolve=true;

def_RA : linsolve( 
    [ uA=0 , vA=0 , wA=0, thetaA=0, phiA=0, psiA=0 ] , 
    [ XA   , YA   , ZA  , ThetaA  , PhiA  , PsiA   ] 
);

define the expected displacement results

utB$
% , def_RA $
% , F=0, I=0 $ /* here I declare the fictitious nature of F and I*/
factor(%) $
utB : % ;

utC$
% , def_RA $
% , F=0, I=0 $
factor(%) $
utC : % ;

let's try to normalize those results
reference displacement
                      | PLACEHOLDER
\\|                   V
\\|====================
\\|
                  <-s-|
  |-------- a --------|

define(
    refdisp(PLACEHOLDER), 
    diff(
        integrate( (PLACEHOLDER*s)^2/2/E/J , s, 0, a ),
        PLACEHOLDER
    )
);

normalized utB and utC

nutB : utB/refipdisp;