kill(all);
make_random_state (true)$ 
/* feed the pseudo-random number generator with a time-based seed*/

placeholder for data to be provided based on the FE analysis
(randomly defined)

xxxFILLMExxx():= random(2.0)-1.0;

silence the "rat: replaced ... " warnings

ratprint: false$

profile area, extracted with "UTILS->MORE->FACE AREA",
and associated to the base metal elements alone

area:499.96826;

in the following, (X,Y) are the in-cross-section coordinates with respect to the model
origin,
(x,y) the in-cross-section coordinates with respect to the center of axial elasticity

## Pure unit z elongation loadcase.

ez=1, kx=0, ky=0.
generic RBE2 master node position.

results from model: hex20profile_ffcd24_ctwsymm_handleatO.mfd

from the output file        
9203   0.0000       0.0000      3.50078E+07 -7.12641E+05 -6.67034E+07   0.0000    

Reaction force Z (EA axial modulus)

EA:3.50078E+07;

Reaction moments X and Y

uezMx :-7.12641E+05;

uezMy :-6.67034E+07;

current coordinates for the center of axial elasticity
(the handle is assumed currently located at (XH,YH)=(0,0) )

[XG,YG] : [-uezMy,+uezMx]/EA;

Before proceeding further, the master node may be shifted accordingly.

procedure to feed the move menu automatically
(the current position of the handle is assumed (0,0,whatever) )

str:
"*move_reset
*set_move_translations ~a ~a 0
"$
ostream:openw("set_move_to_G.proc")$
printf(ostream,str,ev(XG,numer),ev(YG,numer))$
close(ostream)$

after shifting, the reaction forces/moments become
        9203   0.0000       0.0000      3.50078E+07 -0.45097      -38.253       0.0000
the reaction moment comps Mx and My are not vanished, but
the decreased of ~ 6 orders of magnitude

## Pure unit kx and ky curvature loadcases

pure, unit kx curvature (ez=0,kx=1,ky=0).
Centroidal RBE2 master node position.

results from model: 
    2.20673E-11   -8.39407E-12    1.00511E-07    2.75708E+10    5.12261E+07   
0.00000E+00

ukxMx:2.75708E+10;

ukxMy:5.12261E+07;

pure, unit ky curvature (ez=0,kx=0,ky=1).
Centroidal RBE2 master node position

results from model:
    1.40149E-11   -1.77582E-11    6.67967E-09    5.12261E+07    1.12605E+10   
0.00000E+00

ukyMx: 5.12261E+07;

ukyMy:  1.12605E+10;

derived moduli

EJxx:ukxMx$
EJyy:ukyMy$
EJxy:-(ukxMy+ukyMx)/2$

flexural stiffness matrix

EJ: matrix(
        [ EJxx,-EJxy],
        [-EJxy, EJyy]
    );

flexural compliance matrix

invEJ:invert(EJ);

eigendecomposition of the EJ matrix

[[vals,mults],vecs] : eigenvectors(EJ)$

eigenvalues are returned in unknown order, 
let's identify the higher and the lower.

if vals[1] >= vals[2] then (
i1:1,i2:2    
) else (
i1:2,i2:1            
)$

EJ11: vals[i1],numer;

EJ22: vals[i2],numer;

principal directions

vecs[i1][1],numer$ /* eigenvalues are assumed to be distinct*/
dir1 : % / sqrt( %[1]^2 + %[2]^2);

vecs[i2][1],numer$
dir2 : % / sqrt( %[1]^2 + %[2]^2);

Gz rotation of the PROFILE for having the global X,Y aligned with with 1,2,
respectively

rotGz : -atan2(dir1[2],dir1[1]); /*radians*/

rotGz*180/%pi,numer; /*degrees*/

## Pure 
qz=Ez*(a*y+b*x)
global load

Before proceeding further, the following set of b.c. must be defined
- skewsymm. constraints are applied at Z=0
- a free-warping RBE2 link is attached to the profile at Z=l, whose centroidal master
node is simply supported in both the x,y, in-plane directions. 
- I also locked the z rotation to a) let the line of action of the shear forces
deviating from X=XG,Y=YG, Z=l, and b) impose a no-twist, pure shear condition.
- a w=0 positioning constraint inhibits the residual tz rigid body motion.

we scale the qz distribution s.t. the Qx, Qy shear values equate the product of 
* the area (metal only), and of 
* a unit tau (1MPa) shear stress [alternatively, unit eq. vMises stress, and hence
tau=1/sqrt(3)]

uQ:1*area;

unit nominal stress Qx (uQx) case

uQxqz:Ez*matrix([y,-x]).invEJ.matrix([Qy],[-Qx]),Qx=uQ,Qy=0,x=X-XG,y=Y-YG;

unit nominal stress Qy (uQy) case

uQyqz:Ez*matrix([y,-x]).invEJ.matrix([Qy],[-Qx]),Qx=0,Qy=uQ,x=X-XG,y=Y-YG;

create a procedure file to feed mentat

procedure file template, for a single table

str:
"*new_md_table 2 1
*table_name ~a
*set_md_table_type 1
x0_coordinate
*set_md_table_type 2
y0_coordinate
*set_md_table_method_formula
*md_table_formula ~a
"$

procedure file creation

WARNING: possibly pre-existing same-named tables **must** be manually 
         deleted before executing the procedure in Mentat

ostream:openw("set_qz_over_Ez_table.proc")$
printf(ostream,str,"uQxqz_div_Ez",ev(uQxqz,Ez=1,X=v1,Y=v2,numer))$
printf(ostream,str,"uQyqz_div_Ez",ev(uQyqz,Ez=1,X=v1,Y=v2,numer))$
close(ostream)$

moment result collection for the shear center coordinates retrieval
energy result collection for the chi_x chi_y chi_xy definition

unit Qy loadcase

uQyMz:xxxFILLMExxx();

uQyStrainEnergy:xxxFILLMExxx();

unit Qx loadcase

uQxMz:xxxFILLMExxx();

uQxStrainEnergy:xxxFILLMExxx();

both unit Qx Qy loadcase,
energy result collection for the chi_x chi_y chi_xy definition

uQbStrainEnergy:xxxFILLMExxx();

/home/ebertocchi/shear_center_coordinates_for_maxima_FE_evaluation.png

with reference to the image above, Mz should equate the transport moment moving from S
to C, thus leading to Mz=Qx*(-dy)+Qy*(+dx) for any arbitrary Qx,Qy loadings.

linsolve(
    [
         uQyMz= 0*(-dy)+uQ*(+dx),
         uQxMz=uQ*(-dy)+ 0*(+dx)
    ],
    [dx,dy]
),numer,globalsolve=true;

we may now move the RBE2 handle from the current position to the shear center

procedure to feed the move menu automatically

str:
"*move_reset
*set_move_translations ~a ~a 0
"$
ostream:openw("set_move_to_C.proc")$
printf(ostream,str,ev(dx,numer),ev(dy,numer))$
close(ostream)$

Shear contribution to the elastic strain energy (lineic) density 

reference shear modulus, for a reasonable normalization

G:70000/2/(1+0.3)$

length of the profile modeled portion

l:40$

general quadratic for m

u(Qx,Qy) := 1/2 * 
matrix( [Qx,Qy] )
.
matrix( 
        [chix, chixy/2],
        [chixy/2, chiy]
)
.
matrix( 
    [Qx],
    [Qy] 
) / (G*area);

expand(u(Qx,Qy)),area=_A,G=_G;

I substituite the FE reported values for the strain energy in the cases
* unit stress Qx, zero Qy
* unit stress Qy, zero Qx
* unit stress both Qx and Qy (coupled case)

linsolve(
    [
        u(uQ, 0)*l = uQxStrainEnergy,
        u( 0,uQ)*l = uQyStrainEnergy,
        u(uQ,uQ)*l = uQbStrainEnergy
    ],[chix,chixy,chiy]
),globalsolve=true$
[chix,chixy,chiy],numer;

[[vals,mults],vecs]:
eigenvectors(
    matrix( 
        [chix, chixy/2],
        [chixy/2, chiy]
    )
)$

if vals[1] >= vals[2] then (
i1:1,i2:2    
) else (
i1:2,i2:1            
)$

principal chi factors for shear

chi1sh: vals[i1],numer;

chi2sh: vals[i2],numer;

principal directions for shear

vecs[i1][1],numer;
dir1sh : % / sqrt( %[1]^2 + %[2]^2);

vecs[i2][1],numer;
dir2sh : % / sqrt( %[1]^2 + %[2]^2);

## Pure torsion

torsional stiffness calculation: unit twist loadcase

utwMz: xxxFILLMExxx();

GKt : utwMz;