%
x1=0;
y1=0;
%
x2=1;
y2=0;
%
x3=1;
y3=1;
%
x4=0;
y4=1;
%
xc=(x1+x2+x3+x4)/4;
yc=(y1+y2+y3+y4)/4;
%
L12=sqrt((x1-x2)^2+(y1-y2)^2);
L23=sqrt((x2-x3)^2+(y2-y3)^2);
L13=sqrt((x1-x3)^2+(y1-y3)^2);
s=(L12+L23+L13)/2;
area1=sqrt( s*(s-L12)*(s-L23)*(s-L13));
%
L14=sqrt((x1-x4)^2+(y1-y4)^2);
L34=sqrt((x3-x4)^2+(y3-y4)^2);
s=(L14+L34+L13)/2;
area2=sqrt( s*(s-L14)*(s-L34)*(s-L13));
%
area=area1+area2;
%
num=100;
nump=101;
%
dr=(1-(-1))/num;
ds=(1-(-1))/num;
%
mass=zeros(12,12);
stiff=zeros(12,12);
%
%  Stiffness
%
for(i=1:nump)
    r=-1+(i-1)*dr;
    for(j=1:nump)
        s=-1+(j-1)*ds; 
%
        x=(1+r)*(1+s)*x1+(1-r)*(1+s)*x2+(1-r)*(1-s)*x3+(1+r)*(1-s)*x4;
        x=x/4;
        y=(1+r)*(1+s)*y1+(1-r)*(1+s)*y2+(1-r)*(1-s)*y3+(1+r)*(1-s)*y4;
        y=y/4;
%        
        J(1,1)=(1+s)*x1 -(1+s)*x2 -(1-s)*x3 +(1-s)*x4;
        J(1,2)=(1+s)*y1 -(1+s)*y2 -(1-s)*y3 +(1-s)*y4;
        J(2,1)=(1+r)*x1 +(1-r)*x2 -(1-r)*x3 -(1+r)*x4;
        J(2,2)=(1+r)*y1 +(1-r)*y2 -(1-r)*y3 -(1+y)*x4;
        J=J/4;
        Jh=inv(J);
    end
end    
%
%  Mass
%
for(i=1:nump)
    r=-1+(i-1)*dr;
    for(j=1:nump)
       s=-1+(j-1)*ds; 
%
       x=(1+r)*(1+s)*x1+(1-r)*(1+s)*x2+(1-r)*(1-s)*x3+(1+r)*(1-s)*x4;
       x=x/4;
       y=(1+r)*(1+s)*y1+(1-r)*(1+s)*y2+(1-r)*(1-s)*y3+(1+r)*(1-s)*y4;
       y=y/4;
%
       CM(1,1)=1;
       CM(1,2)=y-yc;
       CM(1,3)=x-xc;
       CM(2,2)=(y-yc)^2;
       CM(2,3)=(x-xc)*(y-yc);
       CM(3,3)=(x-xc)^2;
%
       CM(2,1)=CM(1,2);
       CM(3,1)=CM(1,3);
       CM(3,2)=CM(2,3);
%
       PT=zeros(12,3);
%
       PT(1,1)=(1+r)*(1+s);
       PT(2,2)=PT(1,1);
       PT(3,3)=PT(1,1);
%
       PT(4,1)=(1-r)*(1+s);
       PT(5,2)=PT(4,1);
       PT(6,3)=PT(4,1);
%
       PT(7,1)=(1-r)*(1-s);
       PT(8,2)=PT(7,1);
       PT(9,3)=PT(7,1);
%
       PT(10,1)=(1+r)*(1-s);
       PT(11,2)=PT(10,1);
       PT(12,3)=PT(10,1);
       PT=PT/4;
%
       P=PT';
%
       AAA=(PT*(CM*P));
%
       if(i>1 && i <nump)
           AAA=AAA*dr;
       else
           AAA=AAA*dr/2;
       end
       if(j>1 && j <nump)
           AAA=AAA*ds; 
       else
           AAA=AAA*ds/2;       
       end  
       mass=mass+AAA;    
%
    end
end
mass=area*mass/4;
mass
area