kill(all);

assume(a>0,E>0);

aint(f):=
     integrate(integrate( f(x,y) ,x,0*a,2*a),y,0*a,1*a)
    +integrate(integrate( f(x,y) ,x,0*a,1*a),y,1*a,3*a)$

area: aint( lambda([x,y],1) );

SxG: aint( lambda([x,y],x-xG) );

SyG: aint( lambda([x,y],y-yG) );

linsolve([SxG=0,SyG=0],[xG,yG]),globalsolve=true;

Jxx : aint( lambda([x,y],(y-yG)*(y-yG) ) );

Jxy : aint( lambda([x,y],(y-yG)*(x-xG) ) );

Jyy : aint( lambda([x,y],(x-xG)*(x-xG) ) );

Jmat : matrix([E*Jyy,E*Jxy],[E*Jxy,E*Jxx])/(E*Jxx*E*Jyy-(E*Jxy)^2);

[[vals,vmult],vecs]:eigenvectors(Jmat);

linsolve(E*J11=1/vals[1],J11);

linsolve(E*J22=1/vals[2],J22);

vecs[1][1]$
[%[1],%[2]]/sqrt(%[1]^2+%[2]^2)$
[v1x,v1y]:%;

vecs[2][1]$
[%[1],%[2]]/sqrt(%[1]^2+%[2]^2)$
[v2x,v2y]:%;

alpha_rad:atan2(v1y,v1x),numer;
alpha_grad: alpha_rad/%pi*180,numer;

nodes

ni:[1,3,4,9,10]$
nx:[0,2,2,1, 0]*a-xG$
ny:[0,0,1,3, 3]*a-yG$
N:length(ni)$

[Mx,My]:[1,0]*a*(E*a^2)$
[cx,cy]:[Mx*E*Jyy+My*E*Jxy,Mx*E*Jxy+My*E*Jxx]/(E*Jxx*E*Jyy-(E*Jxy)^2);
ne:-cy*nx+cx*ny;
[lmax(ne),lmin(ne)]$
%,numer;

[Mx,My]:[0,1]*a*(E*a^2)$
[cx,cy]:[Mx*E*Jyy+My*E*Jxy,Mx*E*Jxy+My*E*Jxx]/(E*Jxx*E*Jyy-(E*Jxy)^2);
ne:-cy*nx+cx*ny;
[lmax(ne),lmin(ne)]$
%,numer;

[Mx,My]:[v1x,v1y]*a*(E*a^2)$
[cx,cy]:[Mx*E*Jyy+My*E*Jxy,Mx*E*Jxy+My*E*Jxx]/(E*Jxx*E*Jyy-(E*Jxy)^2);
ne:-cy*nx+cx*ny;
[lmax(ne),lmin(ne)]$
%,numer;

[Mx,My]:[v2x,v2y]*a*(E*a^2)$
[cx,cy]:[Mx*E*Jyy+My*E*Jxy,Mx*E*Jxy+My*E*Jxx]/(E*Jxx*E*Jyy-(E*Jxy)^2);
ne:-cy*nx+cx*ny;
[lmax(ne),lmin(ne)]$
%,numer;