function [cmm]=Q4_2Dcmm(X, Y,t, rho)
%clear all;                                   


syms Xi Eta real;



XiC=[-1 1 1 -1];
EtaC=[-1 -1 1 1];

Ni=[ ( (XiC(2)-Xi)*(EtaC(4)-Eta) )/( (XiC(2)-XiC(1))*(EtaC(4)-EtaC(1))),    ((XiC(1)-Xi)*(EtaC(3)-Eta))/( (XiC(1)-XiC(2))*(EtaC(3)-EtaC(2))),      ((XiC(4)-Xi)*(EtaC(2)-Eta))/( (XiC(4)-XiC(3))*(EtaC(2)-EtaC(3))),      ((XiC(3)-Xi)*(EtaC(1)-Eta))/( (XiC(3)-XiC(4))*(EtaC(1)-EtaC(4))) ];

x=Ni*X';
y=Ni*Y';

J=[ diff(x,Xi), diff(y,Xi); diff(x,Eta), diff(y,Eta)];

N=[Ni(1)    0    Ni(2)    0  Ni(3) 0 Ni(4) 0;
       0         Ni(1)  0    Ni(2)  0   Ni(3)  0  Ni(4)];

NTNrhoJ1=subs(N'*N*det(J)*rho*t,{Xi, Eta},{ -1/sqrt(3), -1/sqrt(3)});
NTNrhoJ2=subs(N'*N*det(J)*rho*t,{Xi, Eta},{ 1/sqrt(3), -1/sqrt(3)});
NTNrhoJ3=subs(N'*N*det(J)*rho*t,{Xi, Eta},{ 1/sqrt(3), 1/sqrt(3)});
NTNrhoJ4=subs(N'*N*det(J)*rho*t,{Xi, Eta},{ -1/sqrt(3), 1/sqrt(3)});
cmm=NTNrhoJ1+NTNrhoJ2+NTNrhoJ3+NTNrhoJ4;





   