kill(all);

List of symbols
delta_B  =   d

rollbar intrusion has to be monitored at the point B, as long as the reaction forces at the rollbar supports A and E.
=> d=?, XA=?, YA=?, CA=?

s is a curvilinear abscissa that spans the beam axis from point to point, 
e.g. from A to B, from B to C, etc.

rendo "l" quota dipendendente da "a" e "h"

l : a-h;

linear shape functions in the [0,1] interval (auxiliary functions for interpolation);

N0(0)=1, N0(1)=0
N1(0)=0, N1(1)=1

two different constructions for defining a functions, NOT STRICTLY equivalent in general, equivalent in the present case

the expression at the define second argument is first evaluated, then assigned to N0. 
Hence, evaluation occurs at definition time.
"pippo", a dummy variable, is treated as a placeholder for the first (and only) 
argument passed to the function at call time. 

define(
    N0( pippo ),
    1-pippo
);

define(
    N1( xi ),
    xi
);

controllo si comportino come desiderato

[N0(0),N0(1/2),N0(1)];
[N1(0),N1(1/2),N1(1)];

si assume positivo il momento flettente che tende le fibre interne al portale/rollbar

Bending moment definitions

MfAB :
  ( bending moment contributions at each segment first  extreme A) * N0(xi)
+ ( bending moment contributions at each segment second extreme B) * N1(xi)

bending moment contributions inserted in the following order
. force P alone
. reaction force YA alone
. reaction force XA alone
. reaction moment CA alone
. fictitious force F alone

           |------- first extreme -------|        |------- second extreme ------|
           |                             |        |                             |
           |  P .   YA .   XA .  CA .   F|        |  P .   YA .   XA .  CA .   F|

MfAB   :   (0*P + 0*YA + 0*XA -1*CA  +0*F)*N0(xi)+(0*P  +0*YA  +h*XA -1*CA  +0*F)*N1(xi)$
MfBC   :                     ev(MfAB,xi=1)*N0(xi)+(0*P  +0*YA  +a*XA -1*CA  -l*F)*N1(xi)$
MfCD   :                     ev(MfBC,xi=1)*N0(xi)+(0*P  -b*YA  +a*XA -1*CA  -l*F)*N1(xi)$
MfDE   :                     ev(MfCD,xi=1)*N0(xi)+(a*P  -b*YA  +0*XA -1*CA  +h*F)*N1(xi)$

qui sopra ho sfruttato la PECULIARE continuita' della componente simmetrica del momento flettente rispetto al piano della struttura ("entro piano") del momento flettente

dovrei integrare in ds = [diff(s,xi)] * dxi = [lunghezza_tratto_di_trave]*dxi 

UAB : integrate( MfAB^2/2/E/J * h, xi, 0,1    )$
UBC : integrate( MfBC^2/2/E/J * l, xi, 0,1    )$
UCD : integrate( MfCD^2/2/E/J * b, xi, 0,1    )$
UDE : integrate( MfDE^2/2/E/J * a, xi, 0,1    )$

calcolo l'energia potenziale elastica totale

U : expand(UAB+UBC+UCD+UDE);

utilizzo il secondo teorema di Castigliano per calcolare rotazioni e frecce

dP : diff(U,P)$

dF : diff(U,F)$

uA : diff(U, XA)$

vA : diff(U, YA)$

rA : diff(U, CA)$

valuto le reazioni vincolari imponendo nulli gli spostamenti e le rotazioni all'incastro

linsolve( [ uA=0, vA=0, rA=0 ] , [ XA , YA, CA ]  ) , globalsolve=true;

ho ricavato le reazioni vincolari in funzione di P ed F, e' arrivato il momento di imporre la natura fittizia di F

[dP,dF,XA,YA,CA] : 
factor(
    ev( [dP,dF,XA,YA,CA] , F=0 ,eval, eval,eval)
);

mantengo "a" come unica dimensione caratteristica, ed esprimo le dimensioni
h e b in funzione di rapporti adimensionali
p = h/a
r = b/a

[dP,dF,XA,YA,CA]$
ratsubst( r , b/a , %  )$
ratsubst( p , h/a , %  )$
factor(%);
[dP,dF,XA,YA,CA] : %;

definisco ai fini di normalizzazione dei risultati tipo freccia 
- una struttura di riferimento (il singolo montante direttamente caricato)
- una freccia significativa presa su quella struttura (la freccia in apice)

MfCArif : P * s$
UCArif : integrate( MfCArif^2/2/E/J , s , 0, a)$
dPrif : diff(UCArif, P);

definisco un rapporto di rigidezza relativa tra portale e struttura di riferimento
come rapporto inverso tra le frecce al punto C pari carico

krel : dPrif/dP;

wxplot2d( [ krel ] , [ r , 0.01 , 100 ] , logx );  