kill(all);

/*
relative representation of b and h

h is kept dimensional
b is defined through the upsilon = b/h ratio
t is defined through the epsilon = t/h ratio
*/

b : h * upsilon;

t : h * epsilon;

/*
Dimensions are declared positive, in order to simplify abs(x) with x
*/

assume(upsilon>0,epsilon>0,h>0);

/*
for the sake of simplicity, wall thickness is assumed constant.
*/

/*
two coordinate systems are employed:

OXY is centered at the web midpoint
Gxy is centered at the (initially unknown) section center of gravity
*/

/*
profile nodal point coordinates according to the OXY reference system
*/

XA : b $ YA: h/2$
XB : 0 $ YB: h/2$
XC : 0 $ YC:-h/2$
XD : b $ YD:-h/2$

/*
profile segment lengths
*/

l_AB : sqrt((XA-XB)^2 + (YA-YB)^2)$
l_BC : sqrt((XB-XC)^2 + (YB-YC)^2)$
l_CD : sqrt((XC-XD)^2 + (YC-YD)^2)$

/*
(x,y) coordinates for points along the profile segments are parametrically defined
with respect to a normalized arlength coordinate xi that spans from 0 to 1 
while moving along the segment from the first to the second endpoint.
The Gxy centroidal coordinate system is here employed, with unknown XG,YG. 
*/

x_AB(xi) := (XB-XG)*xi + (XA-XG)*(1-xi)$
y_AB(xi) := (YB-YG)*xi + (YA-YG)*(1-xi)$
x_BC(xi) := (XC-XG)*xi + (XB-XG)*(1-xi)$
y_BC(xi) := (YC-YG)*xi + (YB-YG)*(1-xi)$
x_CD(xi) := (XD-XG)*xi + (XC-XG)*(1-xi)$
y_CD(xi) := (YD-YG)*xi + (YC-YG)*(1-xi)$

/*
Area calculation by integration along the section
*/

area :  integrate(t * l_AB ,xi,0,1)
      + integrate(t * l_BC ,xi,0,1)
      + integrate(t * l_CD ,xi,0,1);

/*
First order (or "statical") moments of area, with respect to the Center of Gravity
*/

Mx  :   integrate(y_AB(xi) * t * l_AB ,xi,0,1)
      + integrate(y_BC(xi) * t * l_BC ,xi,0,1)
      + integrate(y_CD(xi) * t * l_CD ,xi,0,1);

My  :   integrate(x_AB(xi) * t * l_AB ,xi,0,1)
      + integrate(x_BC(xi) * t * l_BC ,xi,0,1)
      + integrate(x_CD(xi) * t * l_CD ,xi,0,1);

/*
the first order moment of area is known to be null
if the reference system is centroidal; such condition is employed
to evaluate the center of gravity coordinates with respect to OXY
*/

linsolve([Mx=0,My=0],[XG,YG]), globalsolve=true;

/*
second order moments of area, with respect to Gxy
*/

Jxx : fullratsimp(
        integrate(y_AB(xi)^2 * t * l_AB ,xi,0,1)
      + integrate(y_BC(xi)^2 * t * l_BC ,xi,0,1)
      + integrate(y_CD(xi)^2 * t * l_CD ,xi,0,1)
    );

Jyy : fullratsimp(
        integrate(x_AB(xi)^2 * t * l_AB ,xi,0,1)
      + integrate(x_BC(xi)^2 * t * l_BC ,xi,0,1)
      + integrate(x_CD(xi)^2 * t * l_CD ,xi,0,1)
    );

Jxy :   integrate(x_AB(xi)*y_AB(xi) * t * l_AB ,xi,0,1)
      + integrate(x_BC(xi)*y_BC(xi) * t * l_BC ,xi,0,1)
      + integrate(x_CD(xi)*y_CD(xi) * t * l_CD ,xi,0,1);

/*
due to the symmetric nature of the cross section, the product (not "mixed") moment
of area is zero.
*/

/*
we now define the alpha,beta coefficients of the equation

d sigma_z
---------  = alpha(x,y,...) Sy - beta(x,y,...) *Sx
    dz

Material is assumed homogeneous and isotropic along the section.
*/

define( alpha(x,y) , fullratsimp ( (-Jxy * x + Jyy*y)/(Jxx*Jyy-Jxy^2) ));

define( beta(x,y) , fullratsimp ( (-Jxx * x + Jxy*y)/(Jxx*Jyy-Jxy^2) ));

/*
for each point along the section, the axial stress component variation 
with increasing z is defined.
Such variation causes a disequilibrium that induces the compensating shear actions. 
*/

define( dsigmaz_dz(x,y) , fullratsimp( alpha(x,y) * Sy - beta(x,y) * Sx ) );

/*
shear stress at the A end, which is a free surface in the case of an open section.
In the case of a closed thin wall section, such value is generally nonzero, and it 
is to be treated as an unknown parameter.
*/

tau_A : 0;

/*
AB segment : The tau_iz shear stress component is evaluated as a function of the
normalized arc length coordinate xi, based on the axial equilibrium of the portion
of section spanning from A (s=0) to the current point s=xi*l_AB along the segment
*/

define(
    tau_AB(xi),
    fullratsimp(
         ( 
              tau_A *t 
            + integrate(
                dsigmaz_dz(x_AB(dummy),y_AB(dummy)) * t * l_AB,
                dummy,0,xi   
              )   
     )/t
    )
);

/*
BC segment: the shear stress component is evaluated as for segment AB; 
at the B segment end, the shear stress transmitted from the previous segment is applied. 

The same normalized arc length coordinate is considered, whereas the absolute arclength 
coordinate (if needed) may be defined as

s = l_AB + xi*l_BC, ds = l_BC * dxi

*/

define(
    tau_BC(xi),
    fullratsimp(
         ( 
              tau_AB(1) *t 
            + integrate(
                dsigmaz_dz(x_BC(dummy),y_BC(dummy)) * t * l_BC,
                dummy,0,xi   
              )   
  )/t
    )
);

/*
BC segment:  the shear stress component is evaluated as for segments AB and CD.

here,

s = l_AB + l_BC + xi*l_CD, ds = l_CD * dxi
*/

define(
    tau_CD(xi),
    fullratsimp(
         ( 
              tau_BC(1) *t 
            + integrate(
             dsigmaz_dz(x_CD(dummy),y_CD(dummy)) * t * l_CD,
                dummy,0,xi   
              )   
          )/t
    )
);

/*
check : is tau_iz equal to zero at the D end??
*/

tau_D : fullratsimp(tau_CD(1));

/*
wow. ok.
*/

/*
further check : let's plot the tau_zs variation along the cross section segments.
sample dimensions are required for extracting numerical results to be plotted.
*/

dim : [h=80,upsilon=40/80,epsilon=6/80];

/*
Shear resultants Sx,Sy are applied for plotting purposes s.t. the nominal shear stress, 
avaluated as the load divided by the area, is unitary.
*/

wxplot2d(
    ev(
        [
            tau_AB(xi),
            tau_BC(xi),
            tau_CD(xi)
        ],dim, Sy=area*1, Sx=0, infeval
    ),
    [xi,0,1],
    [legend,"AB","BC","CD"]
    );

wxplot2d(
    ev(
        [tau_AB(xi),tau_BC(xi),tau_CD(xi)],
        dim,  Sy=0,  Sx=area*1 ,infeval
    ),
    [xi,0,1],
    [legend,"AB","BC","CD"]
    );

/*
I evaluate the peak value for the actual shear stress vs. nominal stress
ratio along the AB segment, and along the section. 
*/

tmp : ev(
            tau_AB(xi),
            dim,  Sy=0,  Sx=area*1 ,infeval
        )$

linsolve(
    diff(
        tmp,
        xi
    )=0,
    xi
)$

ev(tmp , % , numer);

/*
in order to evaluate resultants, a unit vector is to be defined for each
segment that orients the tau_iz actions along the plane.

vx = - diff ( x, s) = -diff (x, xi) / l;
vy = - diff ( y, s) = -diff (y, xi) / l;

those vectors are a function of xi (and should be defined accordingly) 
in the case of curved segments.

the negative sign is due to the opposite "inward" and "growing s" directions; 
tau_zi are oriented inward, whereas the

[ diff ( x, s) , diff ( y, s) ]

vector is oriented towards growing s

*/

vx_AB : - diff(x_AB(xi),xi)/l_AB;
vy_AB : - diff(y_AB(xi),xi)/l_AB;

vx_BC : - diff(x_BC(xi),xi)/l_BC;
vy_BC : - diff(y_BC(xi),xi)/l_BC;

vx_CD : - diff(x_CD(xi),xi)/l_CD;
vy_CD : - diff(y_CD(xi),xi)/l_CD;

/*
x component of the resultant force of the tau_iz stress distribution. 
*/

Fx : integrate(
          vx_AB * tau_AB(xi) * t * l_AB
        + vx_BC * tau_BC(xi) * t * l_BC
        + vx_CD * tau_CD(xi) * t * l_CD,
        xi,0,1
    );

/*
y component of the resultant force of the tau_iz stress distribution. 
*/

Fy : integrate(
          vy_AB * tau_AB(xi) * t * l_AB
        + vy_BC * tau_BC(xi) * t * l_BC
        + vy_CD * tau_CD(xi) * t * l_CD,
        xi,0,1
    )$
Fy : fullratsimp(Fy);

/*
axial z component of the resultant moment of the tau_iz stress 
distribution.
*/

MGz : integrate(
          (vx_AB * (-y_AB(xi)) + vy_AB * (+x_AB(xi))) * tau_AB(xi) * t * l_AB
        + (vx_BC * (-y_BC(xi)) + vy_BC * (+x_BC(xi))) * tau_BC(xi) * t * l_BC
        + (vx_CD * (-y_CD(xi)) + vy_CD * (+x_CD(xi))) * tau_CD(xi) * t * l_CD,
        xi,0,1
    )$
MGz : fullratsimp(MGz);

/*
The evaluation follows of the shear center coordinates (e,f), based on 
an equal resultant moment condition with respect to the CoG.
*/

eq : MGz = Sx * (-f) + Sy * (+e);

eq_res : rhs(eq) - lhs(eq);

/*
such equation holds for null shear components, 

we now impose that such null residual
condition is also constant for arbitrary varying Sx, Sy
*/

ev( eq_res, Sx=0, Sy=0);

linsolve( 
    [
        diff(eq_res,Sx)=0, 
        diff(eq_res,Sy)=0
    ]
,[e,f] ), globalsolve=true;

/*
check based on comparison with known results from literature, e.g.
https://theconstructor.org/structural-engg/analysis/shear-centre-with-examples/3677/
*/

fullratsimp(
      ( -(XG + b^2*h^2*t/4/Jxx) ) 
    - ( e ) 
);

/*
ok
*/

/*
the shear contribution to the elastic strain energy per unit beam length is evaluated 
in the following 
*/

dU_dz :
    integrate(
        tau_AB(xi)^2/2/G * t* l_AB   
      + tau_BC(xi)^2/2/G * t* l_BC   
      + tau_CD(xi)^2/2/G * t* l_CD,
      xi,0,1  
    );

/*
a nominal counterpart is defined, 
that contains the three corrective chi factors
*/

dU_dz_nom : (chi_x * Sx^2 + chi_y * Sy^2 + chi_xy *Sx*Sy)/2/G/area;

eq : dU_dz = dU_dz_nom;

eq_res : lhs(eq) - rhs(eq);

/*
since it is not allowed to differentiate any expression with respect to 
non-atomic terms as Sx^2,Sy^2,Sx*Sy, a workaround is applied in the following
consisting in a change of variable.
*/

eq_res : ratsubst( Sxq  , Sx^2  , eq_res  )$
eq_res : ratsubst( Syq  , Sy^2  , eq_res  )$
eq_res : ratsubst( SxSy , Sx*Sy , eq_res  )$

eq_res;

/*
check: no residual Sx, Sy terms within eq_res
*/

diff( eq_res , Sx);
diff( eq_res , Sy);

/*
ok
*/

/*
final value for the three corrective shear strain energy coefficients
*/

linsolve(
    [
        diff(eq_res, Sxq) =0,
        diff(eq_res, Syq) =0,
        diff(eq_res, SxSy)=0
    ],
    [ chi_x, chi_y, chi_xy]
), globalsolve=true;

ev([ chi_x, chi_y, chi_xy] , dim ,numer);