kill(all);

The reaction forces and moments at the D support are obtained
as a function of F and of the reaction forces and moments 
at A.
THIS PASSAGE IS, FOR THE CASE UNDER SCRUTINY, NOT STRICTLY
REQUIRED FOR THE DEFINITION OF THE STRESS RESULTANT DIAGRAMS

eqrDx: ThetaA+ThetaD+F*b=0;

eqrDz: PsiA + PsiD + F*a-VA*a=0;

eqty: VA+VD-F=0;

linsolve([eqrDx,eqrDz,eqty],[ThetaD,PsiD,VD]), globalsolve=true;

bending moment and torsional moment contributions along the 
various frame members due to F, VA, PsiA, ThetaA;
The other reaction components at a (i.e. UA, WA and PhiA) are
zero in case of skew-symmetric loads with respect to the
frame plane.

Mf_AB : ( 0)*F + (    -s)*VA + ( 0)*PsiA + (-1)*ThetaA$
Mt_AB : ( 0)*F + (     0)*VA + (-1)*PsiA + ( 0)*ThetaA$
Mf_BC : (+s)*F + (    -s)*VA + ( 1)*PsiA + ( 0)*ThetaA$
Mt_BC : ( 0)*F + (    -b)*VA + ( 0)*PsiA + (-1)*ThetaA$
Mf_CD : (+s)*F + (+(b-s))*VA + ( 0)*PsiA + (+1)*ThetaA$
Mt_CD : (+a)*F + (    -a)*VA + (+1)*PsiA + ( 0)*ThetaA$

A single bending moment component is considered (the other is
null due to skew-symmetry).
The case of two nonzero bending moment components is to be
considered in general.
Axial and shear loads are not considered since their internal
energy contributions are assumed negligible.

internal strain energy along the structure

U : 
  integrate( Mf_AB^2/2/E/J + Mt_AB^2/2/G/Kt , s, 0, b)
+ integrate( Mf_BC^2/2/E/J + Mt_BC^2/2/G/Kt , s, 0, a)
+ integrate( Mf_CD^2/2/E/J + Mt_CD^2/2/G/Kt , s, 0, b);

I just obtained the internal energy
 U(F, VA, ThetaA, PsiA )
as a function of the applied load and of the parametric
reactions.

I evaluate the deflections and the rotations at the A
support by employing the second Castigliano theorem;
due to the constrained nature of the associated dofs,
those deflections have to be set as zero.

eq_comp :[
    diff( U , VA     ) = 0,
    diff( U , ThetaA ) = 0,
    diff( U , PsiA   ) = 0 
];

a linear system of equations is obtained, whose solution
will define the reactions

linsolve( eq_comp , [ VA, ThetaA, PsiA ] ) , globalsolve=true;

In order to evaluate the load aligned deflection at B, the strain
energy is reformulated as a function of the external force alone

U(F, VA(F), ThetaA(F), PsiA(F)) = U(F)

U : ev( U ) $

in the absence of the previous instruction, diff(U,F) would be 
erroneously evaluated under the wrong assumptions

diff(VA,F)=0
diff(ThetaA,F)=0
diff(PsiA,F)=0

since the dependence of VA, ThetaA, PsiA on F has not been
explicited in any way.

The second Castigliano Theorem may be applied in evaluating 
the deflection

delta_F : fullratsimp( diff ( U , F ) );

reference dimensions

dim:	[
            a=800,
            b=1000, 
            E=210000,
            G=210000/2/(1+3/10),
            J=(40^4-36^4)*%pi/64,
            Kt=(40^4-36^4)*%pi/32
        ];

ev( [ VA, ThetaA, PsiA, delta_F ] , dim);

numerical evaluation

ev(%,numer);