%
%  thick_plate_mass_stiff.m  ver 1.0  June 6, 2012
%
function[mass_local,stiff_local,area]=thick_plate_mass_stiff(nodex,nodey,node1,node2,node3,node4,mu,E,h,ijk)
%
mass_local =zeros(12,12);
stiff_local=zeros(12,12);
%
x1=nodex(node1(ijk));
x2=nodex(node2(ijk));
x3=nodex(node3(ijk));
x4=nodex(node4(ijk));
%
y1=nodey(node1(ijk));
y2=nodey(node2(ijk));
y3=nodey(node3(ijk));
y4=nodey(node4(ijk));
%
e3=(1-mu)/2;
%
%
khat=5/6;
%
%
a2= ((x3+2*x2-3*x1)*y4+(x2-x4)*y3+(-2*x4-x3+3*x1)*y2+(3*x4-3*x2)*y1)/36;
%
ab= ((x3+x2-2*x1)*y4+(-x4+2*x2-x1)*y3+(-x4-2*x3+3*x1)*y2+(2*x4+x3-3*x2)*y1)/72;
%
ac= ((x3-x1)*y4+(x2-x4)*y3+(x1-x3)*y2+(x4-x2)*y1)/72;
%
ad= ((2*x3+x2-3*x1)*y4+(-2*x4+x2+x1)*y3+(-x4-x3+2*x1)*y2+(3*x4-x3-2*x2)*y1)/72;
%
b2= ((x3-x1)*y4+(-x4+3*x2-2*x1)*y3+(3*x1-3*x3)*y2+(x4+2*x3-3*x2)*y1)/36;
%
bc= ((2*x3-x2-x1)*y4+(-2*x4+3*x2-x1)*y3+(x4-3*x3+2*x1)*y2+(x4+x3-2*x2)*y1)/72;
%
bd=((x3-x1)*y4+(x2-x4)*y3+(x1-x3)*y2+(x4-x2)*y1)/72;
%
c2=((3*x3-2*x2-x1)*y4+(3*x2-3*x4)*y3+(2*x4-3*x3+x1)*y2+(x4-x2)*y1)/36;
%
cd=((3*x3-x2-2*x1)*y4+(-3*x4+2*x2+x1)*y3+(x4-2*x3+x1)*y2+(2*x4-x3-x2)*y1)/72;
%
d2=((3*x3-3*x1)*y4+(-3*x4+x2+2*x1)*y3+(x1-x3)*y2+(3*x4-2*x3-x2)*y1)/36;
%
q=h^2/12;
%
a2q=a2*q; 
%
abq=ab*q; 
%
acq=ac*q; 
%
adq=ad*q; 
%
b2q=b2*q; 
%
bcq=bc*q;
%
bdq=bd*q;
%
c2q=c2*q;
%
cdq=cd*q;
%
d2q=d2*q;
%
mass=[a2,  0,  0, ab,  0,  0, ac,  0,  0, ad,  0,  0;...
       0,a2q,  0,  0,abq,  0,  0,acq,  0,  0,adq,  0;...
       0,  0,a2q,  0,  0,abq,  0,  0,acq,  0,  0,adq;...
      ab,  0,  0, b2,  0,  0, bc,  0,  0, bd,  0,  0;...
       0,abq,  0,  0,b2q,  0,  0,bcq,  0,  0,bdq,  0;...
       0,  0,abq,  0,  0,b2q,  0,  0,bcq,  0,  0,bdq;...
      ac,  0,  0, bc,  0,  0, c2,  0,  0, cd,  0,  0;...
       0,acq,  0,  0,bcq,  0,  0,c2q,  0,  0,cdq,  0;...
       0,  0,acq,  0,  0,bcq,  0,  0,c2q,  0,  0,cdq;...
      ad,  0,  0, bd,  0,  0, cd,  0,  0, d2,  0,  0;...
       0,adq,  0,  0,bdq,  0,  0,cdq,  0,  0,d2q,  0;...
       0,  0,adq,  0,  0,bdq,  0,  0,cdq,  0,  0,d2q];
%
%%% M=mass*rho*h;
%%% M
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Jinv11=-(2*y4+2*y3-2*y2-2*y1)/((x3-x1)*y4+(x2-x4)*y3+(x1-x3)*y2+(x4-x2)*y1);
Jinv12=-(2*y4-2*y3-2*y2+2*y1)/((x3-x1)*y4+(x2-x4)*y3+(x1-x3)*y2+(x4-x2)*y1);
Jinv21= (2*x4+2*x3-2*x2-2*x1)/((x3-x1)*y4+(x2-x4)*y3+(x1-x3)*y2+(x4-x2)*y1);
Jinv22= (2*x4-2*x3-2*x2+2*x1)/((x3-x1)*y4+(x2-x4)*y3+(x1-x3)*y2+(x4-x2)*y1);
%
disp 'ref 3'
%
detJ=((x4-x3-x2+x1)*(-y4-y3+y2+y1))/16-((-x4-x3+x2+x1)*(y4-y3-y2+y1))/16;
%
disp 'ref 4'
%        
a1=@(r,s)(( Jinv11*(1+s)+Jinv12*(1+r))/4);
a2=@(r,s)((-Jinv11*(1+s)+Jinv12*(1-r))/4);  
a3=@(r,s)((-Jinv11*(1-s)-Jinv12*(1-r))/4);
a4=@(r,s)(( Jinv11*(1-s)-Jinv12*(1+r))/4);
%
b1=@(r,s)(( Jinv21*(1+s)+Jinv22*(1+r))/4);
b2=@(r,s)((-Jinv21*(1+s)+Jinv22*(1-r))/4);
b3=@(r,s)((-Jinv21*(1-s)-Jinv22*(1-r))/4);
b4=@(r,s)(( Jinv21*(1-s)-Jinv22*(1+r))/4);
%
core2_2=@(r,s) ((  a1(r,s).^2*e3       +b1(r,s).^2          ) );
core2_3=@(r,s) (( -a1(r,s).*b1(r,s)*mu -a1(r,s).*b1(r,s)*e3 ) );
core2_5=@(r,s) ((  a1(r,s).*a2(r,s)*e3 +b1(r,s).*b2(r,s)    ) );
core2_6=@(r,s) (( -a2(r,s).*b1(r,s)*mu -a1(r,s).*b2(r,s)*e3 ) );
core2_8=@(r,s) ((  a1(r,s).*a3(r,s)*e3 +b1(r,s).*b3(r,s)    ) );
core2_9=@(r,s) (( -a3(r,s).*b1(r,s)*mu -a1(r,s).*b3(r,s)*e3 ) );
core2_11=@(r,s)((  a1(r,s).*a4(r,s)*e3 +b1(r,s).*b4(r,s)    ) );
core2_12=@(r,s)(( -a4(r,s).*b1(r,s)*mu -a1(r,s).*b4(r,s)*e3 ) );
%
core3_3=@(r,s) ((  b1(r,s).^2*e3       +a1(r,s).^2          ) );
core3_5=@(r,s) (( -a1(r,s).*b2(r,s)*mu -a2(r,s).*b1(r,s)*e3 ) );
core3_6=@(r,s) ((  b1(r,s).*b2(r,s)*e3 +a1(r,s).*a2(r,s)    ) );
core3_8=@(r,s) (( -a1(r,s).*b3(r,s)*mu -a3(r,s).*b1(r,s)*e3 ) );
core3_9=@(r,s) ((  b1(r,s).*b3(r,s)*e3 +a1(r,s).*a3(r,s)    ) );
core3_11=@(r,s)(( -a1(r,s).*b4(r,s)*mu -a4(r,s).*b1(r,s)*e3 ) );
core3_12=@(r,s)((  b1(r,s).*b4(r,s)*e3 +a1(r,s).*a4(r,s)    ) );
%
core5_5=@(r,s) ((  a2(r,s).^2*e3       +b2(r,s).^2          ) );
core5_6=@(r,s) (( -a2(r,s).*b2(r,s)*mu -a2(r,s).*b2(r,s)*e3 ) );
core5_8=@(r,s) ((  a2(r,s).*a3(r,s)*e3 +b2(r,s).*b3(r,s)    ) );
core5_9=@(r,s) (( -a3(r,s).*b2(r,s)*mu -a2(r,s).*b3(r,s)*e3 ) );
core5_11=@(r,s)((  a2(r,s).*a4(r,s)*e3 +b2(r,s).*b4(r,s)    ) );
core5_12=@(r,s)(( -a4(r,s).*b2(r,s)*mu -a2(r,s).*b4(r,s)*e3 ) );
%
core6_6=@(r,s) ((  b2(r,s).^2*e3       +a2(r,s).^2           ) );
core6_8=@(r,s) (( -a2(r,s).*b3(r,s)*mu -a3(r,s).*b2(r,s)*e3  ) );
core6_9=@(r,s) ((  b2(r,s).*b3(r,s)*e3 +a2(r,s).*a3(r,s)     ) );
core6_11=@(r,s)(( -a2(r,s).*b4(r,s)*mu -a4(r,s).*b2(r,s)*e3  ) );
core6_12=@(r,s)((  b2(r,s).*b4(r,s)*e3 +a2(r,s).*a4(r,s)     ) );
%
core8_8=@(r,s) ((  a3(r,s).^2*e3       +b3(r,s).^2           ) );
core8_9=@(r,s) (( -a3(r,s).*b3(r,s)*mu -a3(r,s).*b3(r,s)*e3  ) );
core8_11=@(r,s)((  a3(r,s).*a4(r,s)*e3 +b3(r,s).*b4(r,s)     ) );
core8_12=@(r,s)(( -a4(r,s).*b3(r,s)*mu -a3(r,s).*b4(r,s)*e3  ) );
%
core9_9=@(r,s) ((  b3(r,s).^2*e3       +a3(r,s).^2           ) );
core9_11=@(r,s)(( -a3(r,s).*b4(r,s)*mu -a4(r,s).*b3(r,s)*e3  ) );
core9_12=@(r,s)((  b3(r,s).*b4(r,s)*e3 +a3(r,s).*a4(r,s)     ) );
%
core11_11=@(r,s)(( a4(r,s).^2*e3       +b4(r,s).^2           ) );
core11_12=@(r,s)((-a4(r,s).*b4(r,s)*mu -a4(r,s).*b4(r,s)*e3  ) );
%
core12_12=@(r,s)(( b4(r,s).^2*e3       +a4(r,s).^2           ) );
%%
%
disp 'ref 5'
%
k2_2 = dblquad(core2_2,-1,1,-1,1);
k2_3 = dblquad(core2_3,-1,1,-1,1);
k2_5 = dblquad(core2_5,-1,1,-1,1);
k2_6 = dblquad(core2_6,-1,1,-1,1);
k2_8 = dblquad(core2_8,-1,1,-1,1);
k2_9 = dblquad(core2_9,-1,1,-1,1);
k2_11 =dblquad(core2_11,-1,1,-1,1);
k2_12 =dblquad(core2_12,-1,1,-1,1);
%
k3_3 = dblquad(core3_3,-1,1,-1,1);
k3_5 = dblquad(core3_5,-1,1,-1,1);
k3_6 = dblquad(core3_6,-1,1,-1,1);
k3_8 = dblquad(core3_8,-1,1,-1,1);
k3_9 = dblquad(core3_9,-1,1,-1,1);
k3_11 =dblquad(core3_11,-1,1,-1,1);
k3_12 =dblquad(core3_12,-1,1,-1,1);
%
k5_5 = dblquad(core5_5,-1,1,-1,1);
k5_6 = dblquad(core5_6,-1,1,-1,1);
k5_8 = dblquad(core5_8,-1,1,-1,1);
k5_9 = dblquad(core5_9,-1,1,-1,1);
k5_11 =dblquad(core5_11,-1,1,-1,1);
k5_12 =dblquad(core5_12,-1,1,-1,1);
%
k6_6 = dblquad(core6_6,-1,1,-1,1);
k6_8 = dblquad(core6_8,-1,1,-1,1);
k6_9 = dblquad(core6_9,-1,1,-1,1);
k6_11 =dblquad(core6_11,-1,1,-1,1);
k6_12 =dblquad(core6_12,-1,1,-1,1);
%
k8_8 = dblquad(core8_8,-1,1,-1,1);
k8_9 = dblquad(core8_9,-1,1,-1,1);
k8_11 =dblquad(core8_11,-1,1,-1,1);
k8_12 =dblquad(core8_12,-1,1,-1,1);
%
k9_9 = dblquad(core9_9,-1,1,-1,1);
k9_11 =dblquad(core9_11,-1,1,-1,1);
k9_12 =dblquad(core9_12,-1,1,-1,1);
%
k11_11 = dblquad(core11_11,-1,1,-1,1);
k11_12 = dblquad(core11_12,-1,1,-1,1);
%
k12_12 = dblquad(core12_12,-1,1,-1,1);
%
%%%
%
K1=zeros(12,12);
%
K1(2,2)=k2_2;
K1(2,3)=k2_3;
K1(2,5)=k2_5;
K1(2,6)=k2_6;
K1(2,8)=k2_8;
K1(2,9)=k2_9;
K1(2,11)=k2_11;
K1(2,12)=k2_12;
%
K1(3,3)=k3_3;
K1(3,5)=k3_5;
K1(3,6)=k3_6;
K1(3,8)=k3_8;
K1(3,9)=k3_9;
K1(3,11)=k3_11;
K1(3,12)=k3_12;
%
K1(5,5)=k5_5;
K1(5,6)=k5_6;
K1(5,8)=k5_8;
K1(5,9)=k5_9;
K1(5,11)=k5_11;
K1(5,12)=k5_12;
%
K1(6,6)=k6_6;
K1(6,8)=k6_8;
K1(6,9)=k6_9;
K1(6,11)=k6_11;
K1(6,12)=k6_12;
%
K1(8,8)=k8_8;
K1(8,9)=k8_9;
K1(8,11)=k8_11;
K1(8,12)=k8_12;
%
K1(9,9)=k9_9;
K1(9,11)=k9_11;
K1(9,12)=k9_12;
%
K1(11,11)=k11_11;
K1(11,12)=k11_12;
%
K1(12,12)=k12_12;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
v11=@(r,s)(( Jinv21*(1+s)+Jinv22*(1+r))/4);
v12=@(r,s)(( -(1+r).*(1+s))/4);
%
v14=@(r,s)(( -Jinv21*(1+s)+Jinv22*(1-r))/4);
v15=@(r,s)(( -(1-r).*(1+s))/4);
%
v17=@(r,s)(( -Jinv21*(1-s)-Jinv22*(1-r))/4);
v18=@(r,s)(( -(1-r).*(1-s))/4);
%
v110=@(r,s)(( Jinv21*(1-s)-Jinv22*(1+r))/4);
v111=@(r,s)(( -(1+r).*(1-s))/4);
%%
%%
v21=@(r,s)(( Jinv11*(1+s)+Jinv12*(1+r))/4);
%
v23=@(r,s)(( (1+r).*(1+s))/4);
v24=@(r,s)((  -Jinv11*(1+s)+Jinv12*(1-r))/4);
%
v26=@(r,s)(( (1-r).*(1+s))/4);
v27=@(r,s)(( -Jinv11*(1-s)-Jinv12*(1-r))/4);
%
v29=@(r,s)(( (1-r).*(1-s))/4);
v210=@(r,s)(( Jinv11*(1-s)-Jinv12*(1+r))/4);
%
v212=@(r,s)(( (1+r).*(1-s))/4);
%
%%
%
kore1_1=@(r,s) ((  v21(r,s).^2      +v11(r,s).^2            ) );
kore1_2=@(r,s) ((  v11(r,s).*v12(r,s)                       ) );
kore1_3=@(r,s) ((  v21(r,s).*v23(r,s)                       ) );
kore1_4=@(r,s) ((  v21(r,s).*v24(r,s) + v11(r,s).*v14(r,s)  ) );
kore1_5=@(r,s) ((  v11(r,s).*v15(r,s)                       ) );
kore1_6=@(r,s) ((  v21(r,s).*v26(r,s)                       ) );
kore1_7=@(r,s) ((  v21(r,s).*v27(r,s) + v11(r,s).*v17(r,s)  ) );
kore1_8=@(r,s) ((  v11(r,s).*v18(r,s)                       ) );
kore1_9=@(r,s) ((  v21(r,s).*v29(r,s)                       ) );
kore1_10=@(r,s)((  v21(r,s).*v210(r,s) + v11(r,s).*v110(r,s)) );
kore1_11=@(r,s)((  v11(r,s).*v111(r,s)                      ) );
kore1_12=@(r,s)((  v21(r,s).*v212(r,s)                      ) );
%
kore2_2=@(r,s) ((  v12(r,s).^2                              ) );

kore2_4=@(r,s) ((  v12(r,s).*v14(r,s)                       ) );
kore2_5=@(r,s) ((  v12(r,s).*v15(r,s)                       ) );

kore2_7=@(r,s) ((  v12(r,s).*v17(r,s)                       ) );
kore2_8=@(r,s) ((  v12(r,s).*v18(r,s)                       ) );
kore2_10=@(r,s)((  v110(r,s).*v12(r,s)                      ) );
kore2_11=@(r,s)((  v111(r,s).*v12(r,s)                      ) );

%
kore3_3=@(r,s) ((  v23(r,s).^2                              ) );
kore3_4=@(r,s) ((  v23(r,s).*v24(r,s)                       ) );

kore3_6=@(r,s) ((  v23(r,s).*v26(r,s)                       ) );
kore3_7=@(r,s) ((  v23(r,s).*v27(r,s)                       ) );

kore3_9=@(r,s) ((  v23(r,s).*v29(r,s)                       ) );
kore3_10=@(r,s)((  v210(r,s).*v23(r,s)                      ) );

kore3_12=@(r,s)((  v212(r,s).*v23(r,s)                      ) );
%
kore4_4=@(r,s) ((  v24(r,s).^2      +v14(r,s).^2            ) );
kore4_5=@(r,s) ((  v14(r,s).*v15(r,s)                       ) );
kore4_6=@(r,s) ((  v24(r,s).*v26(r,s)                       ) );
kore4_7=@(r,s) ((  v24(r,s).*v27(r,s) + v14(r,s).*v17(r,s)  ) );
kore4_8=@(r,s) ((  v14(r,s).*v18(r,s)                       ) );
kore4_9=@(r,s) ((  v24(r,s).*v29(r,s)                       ) );
kore4_10=@(r,s)((  v210(r,s).*v24(r,s) + v110(r,s).*v14(r,s)) );
kore4_11=@(r,s)((  v111(r,s).*v14(r,s)                      ) );
kore4_12=@(r,s)((  v212(r,s).*v24(r,s)                      ) );
%
kore5_5=@(r,s) ((  v15(r,s).^2                              ) );

kore5_7=@(r,s) ((  v15(r,s).*v17(r,s)                       ) );
kore5_8=@(r,s) ((  v15(r,s).*v18(r,s)                       ) );

kore5_10=@(r,s)((  v110(r,s).*v15(r,s)                      ) );
kore5_11=@(r,s)((  v111(r,s).*v15(r,s)                      ) );

%
kore6_6=@(r,s) ((  v26(r,s).^2                              ) );
kore6_7=@(r,s) ((  v26(r,s).*v27(r,s)                       ) );

kore6_9=@(r,s) ((  v26(r,s).*v29(r,s)                       ) );
kore6_10=@(r,s)((  v210(r,s).*v26(r,s)                      ) );

kore6_12=@(r,s)((  v212(r,s).*v26(r,s)                      ) );
%
kore7_7=@(r,s) ((  v27(r,s).^2      +v17(r,s).^2            ) );
kore7_8=@(r,s) ((  v17(r,s).*v18(r,s)                       ) );
kore7_9=@(r,s) ((  v27(r,s).*v29(r,s)                       ) );
kore7_10=@(r,s)((  v210(r,s).*v27(r,s) + v110(r,s).*v17(r,s)) );
kore7_11=@(r,s)((  v111(r,s).*v17(r,s)                      ) );
kore7_12=@(r,s)((  v212(r,s).*v27(r,s)                      ) );
%
kore8_8=@(r,s) ((  v18(r,s).^2                              ) );

kore8_10=@(r,s)((  v110(r,s).*v18(r,s)                      ) );
kore8_11=@(r,s)((  v111(r,s).*v18(r,s)                      ) );

%
kore9_9=@(r,s) ((  v29(r,s).^2                              ) );
kore9_10=@(r,s)((  v210(r,s).*v29(r,s)                      ) );

kore9_12=@(r,s)((  v212(r,s).*v29(r,s)                      ) );
%
kore10_10=@(r,s)((  v210(r,s).^2      +v110(r,s).^2         ) );
kore10_11=@(r,s)((  v110(r,s).*v111(r,s)                    ) );
kore10_12=@(r,s)((  v210(r,s).*v212(r,s)                    ) );
%
kore11_11=@(r,s)((  v111(r,s).^2                            ) );

%
kore12_12=@(r,s)((  v212(r,s).^2                            ) );
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
K2(1,1)=dblquad(kore1_1,-1,1,-1,1);
K2(1,2)=dblquad(kore1_2,-1,1,-1,1);
K2(1,3)=dblquad(kore1_3,-1,1,-1,1);
K2(1,4)=dblquad(kore1_4,-1,1,-1,1);
K2(1,5)=dblquad(kore1_5,-1,1,-1,1);
K2(1,6)=dblquad(kore1_6,-1,1,-1,1);
K2(1,7)=dblquad(kore1_7,-1,1,-1,1);
K2(1,8)=dblquad(kore1_8,-1,1,-1,1);
K2(1,9)=dblquad(kore1_9,-1,1,-1,1);
K2(1,10)=dblquad(kore1_10,-1,1,-1,1);
K2(1,11)=dblquad(kore1_11,-1,1,-1,1);
K2(1,12)=dblquad(kore1_12,-1,1,-1,1);
%
K2(2,2)=dblquad(kore2_2,-1,1,-1,1);

K2(2,4)=dblquad(kore2_4,-1,1,-1,1);
K2(2,5)=dblquad(kore2_5,-1,1,-1,1);

K2(2,7)=dblquad(kore2_7,-1,1,-1,1);
K2(2,8)=dblquad(kore2_8,-1,1,-1,1);

K2(2,10)=dblquad(kore2_10,-1,1,-1,1);
K2(2,11)=dblquad(kore2_11,-1,1,-1,1);

%
K2(3,3)=dblquad(kore3_3,-1,1,-1,1);
K2(3,4)=dblquad(kore3_4,-1,1,-1,1);

K2(3,6)=dblquad(kore3_6,-1,1,-1,1);
K2(3,7)=dblquad(kore3_7,-1,1,-1,1);

K2(3,9)=dblquad(kore3_9,-1,1,-1,1);
K2(3,10)=dblquad(kore3_10,-1,1,-1,1);

K2(3,12)=dblquad(kore3_12,-1,1,-1,1);
%
K2(4,4)=dblquad(kore4_4,-1,1,-1,1);
K2(4,5)=dblquad(kore4_5,-1,1,-1,1);
K2(4,6)=dblquad(kore4_6,-1,1,-1,1);
K2(4,7)=dblquad(kore4_7,-1,1,-1,1);
K2(4,8)=dblquad(kore4_8,-1,1,-1,1);
K2(4,9)=dblquad(kore4_9,-1,1,-1,1);
K2(4,10)=dblquad(kore4_10,-1,1,-1,1);
K2(4,11)=dblquad(kore4_11,-1,1,-1,1);
K2(4,12)=dblquad(kore4_12,-1,1,-1,1);
%
K2(5,5)=dblquad(kore5_5,-1,1,-1,1);

K2(5,7)=dblquad(kore5_7,-1,1,-1,1);
K2(5,8)=dblquad(kore5_8,-1,1,-1,1);

K2(5,10)=dblquad(kore5_10,-1,1,-1,1);
K2(5,11)=dblquad(kore5_11,-1,1,-1,1);

%
K2(6,6)=dblquad(kore6_6,-1,1,-1,1);
K2(6,7)=dblquad(kore6_7,-1,1,-1,1);

K2(6,9)=dblquad(kore6_9,-1,1,-1,1);
K2(6,10)=dblquad(kore6_10,-1,1,-1,1);

K2(6,12)=dblquad(kore6_12,-1,1,-1,1);
%
K2(7,7)=dblquad(kore7_7,-1,1,-1,1);
K2(7,8)=dblquad(kore7_8,-1,1,-1,1);
K2(7,9)=dblquad(kore7_9,-1,1,-1,1);
K2(7,10)=dblquad(kore7_10,-1,1,-1,1);
K2(7,11)=dblquad(kore7_11,-1,1,-1,1);
K2(7,12)=dblquad(kore7_12,-1,1,-1,1);
%
K2(8,8)=dblquad(kore8_8,-1,1,-1,1);

K2(8,10)=dblquad(kore8_10,-1,1,-1,1);
K2(8,11)=dblquad(kore8_11,-1,1,-1,1);

%
K2(9,9)=dblquad(kore9_9,-1,1,-1,1);
K2(9,10)=dblquad(kore9_10,-1,1,-1,1);

K2(9,12)=dblquad(kore9_12,-1,1,-1,1);
%
K2(10,10)=dblquad(kore10_10,-1,1,-1,1);
K2(10,11)=dblquad(kore10_11,-1,1,-1,1);
K2(10,12)=dblquad(kore10_12,-1,1,-1,1);
%
K2(11,11)=dblquad(kore11_11,-1,1,-1,1);

%
K2(12,12)=dblquad(kore12_12,-1,1,-1,1);
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 
K1=K1*detJ;
K2=K2*detJ;
%
K1=K1*(E*h^3)/(12*(1-mu^2));
K2=K2*(E*h*khat)/(2*(1+mu));
%
%  Make stiffness matrices symmetric
%
K=K1+K2;
%
num=12;
%
for i=1:num
    for j=i+1:num
        K(j,i)=K(i,j);        
    end
end
%
area=4*detJ;
%
mass_local=mass;
stiff_local=K;