3D elasticity complex with mixed elements

[jmv2009]

3D elasticity complex with mixed elements with weak symmetry enforcement

(only "primal" method, based on only displacement dofs, is demonstrated in Freefem documentation)

Describe the solution you'd like
Would like to see Nedelec type II (linear face elements) for the elasticity complex in 3D.

Or better, actually all of these:
https://web.archive.org/web/20190325195633/http://femtable.org/femtable.pdf

Below is lowest order elasticity complex with mixed elements with weak symmetry, enforced by Lagrange multipliers. Literature ref:
https://www.doi.org/10.1007/0-387-38034-5_3

Below is 2D example:
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

load "Element_Mixte";

int nn=5;
mesh Th= square(4*nn,nn,[4*x,y]);
fespace Vh(Th, [BDM1,BDM1,P0,P0,P0]);

real young=2;
real nu=0.3;

real mu=young/2/(1+nu);
real lambda=young*nu/(1+nu)/(1-2*nu);

real mu1= 0.5/mu;
real lambda1 = -mu1*lambda/(2*lambda+2*mu);
real rho = 1;
real gravity = -1*rho;
 
Vh [usxx,usyx,usxy,usyy,ux,uy,up],[vsxx,vsyx,vsxy,vsyy,vx,vy,vp];
solve fld([usxx,usyx,usxy,usyy,ux,uy,up],[vsxx,vsyx,vsxy,vsyy,vx,vy,vp])
  = int2d(Th)((mu1+lambda1)*usxx*vsxx+mu1*usxy*vsxy+mu1*usyx*vsyx+(mu1+lambda1)*usyy*vsyy+lambda1*usxx*vsyy+lambda1*usyy*vsxx+ux*(dx(vsxx)+dy(vsyx))+uy*(dx(vsxy)+dy(vsyy))+(dx(usxx)+dy(usyx))*vx+(dx(usxy)+dy(usyy))*vy+(usxy-usyx)*vp+up*(vsxy-vsyx)  
  )
  +  int2d(Th) ( -gravity*vy)
  +on(1,2,3,usxx=0,usyx=0,usxy=0,usyy=0)
 //+on(4,ux=0,uy=0)
  ;
 plot(ux, wait=1);
 mesh th = movemesh(Th, [x+0.001*ux, y+0.001*uy]);
plot(th, wait=true);