kill(all);

u:Aux*x+Auy*y+bu*1;
v:Avx*x+Avy*y+bv*1;
w:Awxx*x*x+Awxy*x*y+Awyy*y*y+Awx*x+Awy*y+bw*1;
θ:Aθx*x+Aθy*y+bθ*1;
φ:Aφx*x+Aφy*y+bφ*1;
ψ:bψ*1;

atO:[x=0,y=0];

linsolve(
    [ 
        ex = diff(u,x),
        ey = diff(v,y),
        gxy= diff(u,y)+diff(v,x),
        ψO=-diff(u,y)+diff(v,x), /*rotation term around Oz*/
        uO=ev(u,atO),
        vO=ev(v,atO)        
    ],
    [
        Aux,Avy,Auy,Avx,bu,bv
    ]
),globalsolve=true;

the out-of-plane shear components induce nonuniform bending in the absence of distributed actions, and hence they must be zero

[gzx,gzy]  :[0,0]$

defgzx:diff(w,x)-(gzx-φ);
defgzy:diff(w,y)-(gzy+θ);

sols:linsolve(
    [ 
        kx = +diff(φ,x),
        ky = -diff(θ,y),
        kxy= +diff(φ,y)-diff(θ,x),
        ev(defgzx,atO)=0, /* the defgzx equation must be verified for all (x,y) */
        diff(defgzx,x)=0,
        diff(defgzx,y)=0,        
        ev(defgzy,atO)=0, /* the defgzy equation equation must be verified for all (x,y) */
        diff(defgzy,x)=0,
        diff(defgzy,y)=0,        
        wO=ev(w,atO),
        θO=ev(θ,atO),
        φO=ev(φ,atO),        
        ψO=ev(ψ,atO)       
    ],
    [
        Aθx,Aθy,bθ,
        Aφx,Aφy,bφ,
        Awxx,Awxy,Awyy,Awx,Awy,bw,
        bψ
        
    ]
),globalsolve=true;

reference rigid body motion, sampled at the origin

[uO,vO,wO]:[0,0,0]$
[θO,φO,ψO]:[0,0,0]$

reference uniform deformation values

[ex,ey,gxy]:[0,0,0]$
[kx,ky,kxy]:[-1/6,1,4/6]*1e-3$

reference thickness and offset

h:0.001$ o:0$

reference mesh

atA:[x=0   ,y=0   ]$
atB:[x=0.24,y=0   ]$
atC:[x=0.24,y=0.12]$
atD:[x=0   ,y=0.12]$

[u,v,w,θ,φ,ψ]:ev([u,v,w,θ,φ,ψ],infeval)$

displacement and rotation fields

disps:[u,v,w];

rots:[θ,φ,ψ];

values at the ABCD nodes

[disps,rots],atA;

[disps,rots],atB;

[disps,rots],atC;

[disps,rots],atD;

in-plane strain values

ztop:o+h/2$
zmid:o$
zbtm:o-h/2$

top, mid and bottom values

[ex,ey,gxy]+[kx,ky,kxy]*ztop;

[ex,ey,gxy]+[kx,ky,kxy]*zmid;

[ex,ey,gxy]+[kx,ky,kxy]*zbtm;