kill(all);

make floating point <-> exact rational arithmetics warning silent.

ratprint:false;

define a general linear in x, y (not z)
displacement (u,v only) / rotation field
the drillin rot. is constant
the w displacement component is actually quadratic

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;  

hypotheses related to given locations, "O"rigin

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

the displacement coefficients are defined based on the desired
strain components, and rigid body motions

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 in the Irons' patch test
which considers 
i) uniform membrane stretching
ii) uniform curvature

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

defgzx:diff(w,x)-(gzx-φ); /* =0 from the plate theory paragraph */
defgzy:diff(w,y)-(gzy+θ); /*    residual of such equation */

cast a lot of equations in order to obtain 
all the residual unknown coefficients
from curvature values and oop shear components

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;

All the coefficients that define the displacement and the rotation field
are now expressed as a function of 
- membrane strains
- curvatures
- positioning constraints

now I can assign arbitrary values to each of the aforementioned quantities

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$

refer to a specific geometry for the patch, from MacNeal1985proposed

reference thickness and offset

h:0.001$ o:0$

reference rectangular domain, and its area

lx:0.24$ly:0.12$

patch perimeter corners

atA:[x=0 , y=0   ]$
atB:[x=lx, y=0   ]$
atC:[x=lx, y=ly  ]$
atD:[x=0 , y=ly  ]$
atE:[x=lx/2 , y=0  ]$ /*AB midpoint*/
atF:[x=lx   , y=ly/2  ]$ /*BC midpoint*/
atG:[x=lx/2 , y=ly  ]$ /* ... */
atH:[x=0    , y=ly/2  ]$

area:lx*ly$

[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;

[disps,rots],atE;

[disps,rots],atF;

[disps,rots],atG;

[disps,rots],atH;

Now I evaluate the exact value of the strain energy according
the the aforementione strain/curvature field and patch size
plus strain values at top, mid, bottom

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;

material properties

E:1000000$ nu:1/4$
G:E/2/(1+nu)$

plane stress isotropic hooke law stiffness matrix

D:E/(1-nu^2)*matrix([1,nu,0],[nu,1,0],[0,0,0])+G*matrix([0,0,0],[0,0,0],[0,0,1])$

laminate constitutive matrix minors

a:integrate(D*1  ,z,o-h/2,o+h/2)$
b:integrate(D*z  ,z,o-h/2,o+h/2)$ bt:transpose(b)$
c:integrate(D*z^2,z,o-h/2,o+h/2)$

matrix form for "e" and "k"

e:matrix([ex],[ey],[gxy])$ et:transpose(e)$
k:matrix([kx],[ky],[kxy])$ kt:transpose(k)$

total strain energy

U:1/2*(et.(a.e+b.k)+kt.(bt.e+c.k))*area;