%
%  triangular_plate_local.m  ver 1.0  August 21, 2012
%
function[kl,ml,area]=triangular_plate_local(D,x1,x2,x3,y1,y2,y3,term_integ)
%
area=x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2);
area=abs(area/2);
TA=2*area;
%
b1=y2-y3;
b2=y3-y1;
b3=y1-y2;
%
c1=x3-x2;
c2=x1-x3;
c3=x2-x1;
%
Len1=sqrt( (x3-x2)^2 + (y3-y2)^2 );
Len2=sqrt( (x1-x3)^2 + (y1-y3)^2 );
Len3=sqrt( (x2-x1)^2 + (y2-y1)^2 );
%
mu1=(Len3^2-Len2^2)/Len1^2;
mu2=(Len1^2-Len3^2)/Len2^2;
mu3=(Len2^2-Len1^2)/Len3^2;
%
%%%%%%%%%%%%%%%%
%
%% out1=sprintf('\n x1=%8.4g  x2=%8.4g  x3=%8.4g ',x1,x2,x3);
%% disp(out1);
%% out1=sprintf(' y1=%8.4g  y2=%8.4g  y3=%8.4g \n',y1,y2,y3);
%% disp(out1);
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Local Mass Matrix
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
term_integ=term_integ*TA;
%
%
for i=1:5
    for j=1:5
        for k=1:5
            ii=i-1;
            jj=j-1;
            kk=k-1;     
            data=term_integ(i,j,k);
            index=sprintf('term_%d_%d_%d',ii,jj,kk);
            tt.(index)=data;
        end
    end
end    
%
m=zeros(9,9);
%
%  1,1
m(1,1)=tt.term_2_0_0;
%
% 1,2
m(1,2)=tt.term_1_1_0;
% 
% 1,3
m(1,3)=tt.term_1_0_1;
%
% 1,4
m(1,4)=tt.term_2_1_0;
% 
% 1,5
m(1,5)=tt.term_1_1_1;
%
% 1,6
m(1,6)=tt.term_2_0_1;
% 
% 1,7
m(1,7)= ((1-3*mu3)*tt.term_3_1_1)/2 ...
       +tt.term_3_1_0;
%
% 1,8
m(1,8)=+((1-mu1)*tt.term_2_2_1)/2 ...
       +tt.term_1_2_1 ...
       +((1+3*mu1)*(tt.term_3_1_1-tt.term_2_1_2))/2;
% 
% 1,9
%
m(1,9)= +((1-mu2)*tt.term_2_1_2)/2 ...
        +tt.term_2_0_2 ...
        +((1+3*mu2)*(tt.term_2_2_1-tt.term_3_1_1))/2;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 2,2
m(2,2)=tt.term_0_2_0;
% 
% 2,3
m(2,3)=tt.term_0_1_1;
%
% 2,4
m(2,4)=tt.term_1_2_0;
% 
% 2,5
m(2,5)=tt.term_0_2_1;
%
% 2,6
m(2,6)=tt.term_1_1_1;
% 
% 2,7
m(2,7)= ((1+3*mu3)*(tt.term_1_2_2-tt.term_1_3_1))/2 ...
        +((1-mu3)*tt.term_2_2_1)/2 ...
        +tt.term_2_2_0; 
%
% 2,8
m(2,8)= +((1-mu1)*tt.term_1_3_1)/2 ...
        +tt.term_0_3_1;
% 
% 2,9
m(2,9)= ((1-mu2)*tt.term_1_2_2)/2 ...
        +tt.term_1_1_2 ...
        +((1+3*mu2)*(tt.term_1_3_1-tt.term_2_2_1))/2;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 3,3
m(3,3)=tt.term_0_0_2;
%
% 3,4
m(3,4)=tt.term_1_1_1;
% 
% 3,5
m(3,5)=tt.term_0_1_2;
%
% 3,6
m(3,6)=tt.term_1_0_2;
% 
% 3,7
m(3,7)= ((1+3*mu3)*(tt.term_1_1_3-tt.term_1_2_2))/2 ...
       +((1-mu3)*tt.term_2_1_2)/2 ...
       +tt.term_2_1_1;
%
% 3,8
m(3,8)= +((1-mu1)*tt.term_1_2_2)/2 ...
        +tt.term_0_2_2 ...
        +((1+3*mu1)*(tt.term_2_1_2-tt.term_1_1_3))/2;
% 
% 3,9
%
m(3,9)= +((1-mu2)*tt.term_1_1_3)/2 ...  
        +tt.term_1_0_3;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 4,4
m(4,4)=tt.term_2_2_0;
% 
% 4,5
m(4,5)=tt.term_1_2_1;
%
% 4,6
m(4,6)=tt.term_2_1_1;
% 
% 4,7
m(4,7)=  ((1+3*mu3)*(tt.term_2_2_2-tt.term_2_3_1))/2 ...
        +((1-mu3)*tt.term_3_2_1)/2 ...
        +tt.term_3_2_0;
%
% 4,8
m(4,8)= +((1-mu1)*tt.term_2_3_1)/2 ...
        +tt.term_1_3_1 ...
        +((1+3*mu1)*(tt.term_3_2_1-tt.term_2_2_2))/2;
% 
% 4,9
%
m(4,9)= +((1-mu2)*tt.term_2_2_2)/2 ...
        +tt.term_2_1_2;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 5,5
m(5,5)=tt.term_0_2_2;
%
% 5,6
m(5,6)=tt.term_1_1_2;
% 
% 5,7
m(5,7)= +((1-mu3)*tt.term_2_2_2)/2 ...
        +tt.term_2_2_1;
%
% 5,8
m(5,8)= +((1-mu1)*tt.term_1_3_2)/2 ...
        +tt.term_0_3_2 ...
        +((1+3*mu1)*(tt.term_2_2_2-tt.term_1_2_3))/2;
% 
% 5,9
m(5,9)= +((1-mu2)*tt.term_1_2_3)/2 ...       
        +tt.term_1_1_3 ...
        +((1+3*mu2)*(tt.term_1_3_2-tt.term_2_2_2))/2;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 6,6
m(6,6)=tt.term_2_0_2;
% 
% 6,7
m(6,7)= ((1+3*mu3)*(tt.term_2_1_3-tt.term_2_2_2))/2 ...
       +((1-mu3)*tt.term_3_1_2)/2 ...
       +tt.term_3_1_1;
%
% 6,8
m(6,8)= +((1-mu1)*tt.term_2_2_2)/2 ...
        +tt.term_1_2_2;
% 
% 6,9
%
m(6,9)=  +((1-mu2)*tt.term_2_1_3)/2 ...  
         +tt.term_2_0_3 ...  
         +((1+3*mu2)*(tt.term_2_2_2-tt.term_3_1_2))/2;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 7,7
m(7,7)= +(2.25*mu3^2+1.5*mu3+0.25)*(tt.term_2_2_4+tt.term_2_4_2) ...
        -(4.5*mu3^2+3*mu3+0.5)*(tt.term_2_3_3) ...
        +(0.25*mu3^2-0.5*mu3+0.25)*(tt.term_4_2_2) ...   
        +(1+3*mu3)*(tt.term_3_2_2-tt.term_3_3_1) ... 
        +(1-mu3)*tt.term_4_2_1 ...
        +tt.term_4_2_0;
%
% 7,8
m(7,8)= -(2.25*mu1*mu3+0.75*mu3+0.75*mu1+0.25)*tt.term_2_2_4 ...
        +(1.5*mu1*mu3+1.5*mu3+0.5*mu1+0.5 )*tt.term_2_3_3 ...
        +(1.5*mu3+0.5)*(tt.term_1_3_3-tt.term_1_4_2) ...
        +(3.0*mu1*mu3+1.0*mu3)*tt.term_3_2_3 ...
        +(0.75*mu1*mu3-0.75*mu3+0.25*mu1-0.25)*tt.term_2_4_2 ...   
        -(2.0*mu1*mu3+1.0*mu3+1.0*mu1)*tt.term_3_3_2 ...
        +(-0.5*mu3+0.5)*tt.term_2_3_2 ...
        +(-0.75*mu1*mu3-0.25*mu3+0.75*mu1+0.25)*tt.term_4_2_2 ...
        +(-0.5*mu1+0.5 )*tt.term_3_3_1 ...
        +tt.term_2_3_1 ...
        +(1.5*mu1+0.5)*(tt.term_4_2_1-tt.term_3_2_2);
% 
% 7,9
m(7,9)= +(-0.75*mu2*mu3+0.75*mu3-0.25*mu2+0.25)*tt.term_2_2_4 ...
        +(1.5*mu3+0.5)*(tt.term_2_1_4-tt.term_2_2_3)...
        +(3.0*mu2*mu3+1.0*mu2)*tt.term_2_3_3 ...
        -(2.0*mu2*mu3+1.0*mu3+1.0*mu2)*tt.term_3_2_3 ...
        +(-0.5*mu3+0.5)*tt.term_3_1_3 ...
        -(2.25*mu2*mu3+0.75*mu3+0.75*mu2+0.25)*tt.term_2_4_2 ...
        +(1.5*mu2*mu3+0.5*mu3+1.5*mu2+0.5)*tt.term_3_3_2 ...
        +(0.75*mu2*mu3+0.25*mu3-0.75*mu2-0.25)*tt.term_4_2_2 ...
        +(-0.5*mu2+0.5)*tt.term_3_2_2 ...
        +tt.term_3_1_2 ...
        +(1.5*mu2+0.5)*(tt.term_3_3_1-tt.term_4_2_1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 8,8
m(8,8)=  +(2.25*mu1^2+1.5*mu1+0.25)*(tt.term_2_2_4+tt.term_4_2_2) ...
         -(1+3*mu1)*tt.term_1_3_3 ...
         -(4.5*mu1^2+3*mu1+0.5)*tt.term_3_2_3 ...
         +(0.25*mu1^2-0.5*mu1+0.25)*tt.term_2_4_2 ...
         +(1-mu1)*tt.term_1_4_2 ...
         +tt.term_0_4_2 ...
         +(1+3*mu1)*tt.term_2_3_2;
% 
% 8,9  
m(8,9)=  0.25*(3*mu1*mu2+mu2-3*mu1-1)*tt.term_2_2_4 ...
        -0.5*(1+3*mu1)*tt.term_2_1_4 ...
        -(2*mu1*mu2+mu2+mu1)*tt.term_2_3_3 ...
        +0.5*((1-mu2))*tt.term_1_3_3 ...
        +0.5*(3*mu1*mu2+mu2+3*mu1+1)*tt.term_3_2_3 ...
        +0.5*(1-mu1)*tt.term_2_2_3 ... 
        +tt.term_1_2_3 ...
        +0.5*((1+3*mu1)*tt.term_3_1_3) ...
        +0.25*(-3*mu1*mu2+3*mu2-mu1+1)*tt.term_2_4_2 ...
        +0.5*(1+3*mu2)*tt.term_1_4_2 ...
        +(3*mu1*mu2+mu1)*tt.term_3_3_2 ...
        -0.5*(1+3*mu2)*tt.term_2_3_2 ...
        -0.25*(9*mu1*mu2+3*mu2+3*mu1+1)*tt.term_4_2_2;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 9,9  
m(9,9)= (0.25*mu2^2-0.5*mu2+0.25)*(tt.term_2_2_4) ...  
       +(1-mu2)*tt.term_2_1_4 ...
       +tt.term_2_0_4 ...
       +(1+3*mu2)*(tt.term_2_2_3-tt.term_3_1_3)...
       +(2.25*mu2^2+1.5*mu2+0.25)*(tt.term_2_4_2+tt.term_4_2_2) ...
       -(4.5*mu2^2+3*mu2+0.5)*(tt.term_3_3_2);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Local Stiffness Matrix
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%
%
TA=b1*c2-b2*c1;

%%%
%
d1d1_1=@(L1,L2,L3)0;
d1d1_2=@(L1,L2,L3)0;
d1d1_3=@(L1,L2,L3)0;
d1d1_4=@(L1,L2,L3)0;
d1d1_5=@(L1,L2,L3)0;
d1d1_6=@(L1,L2,L3)0;
d1d1_7=@(L1,L2,L3)(1-mu3)*L2*L3+2*L2;
d1d1_8=@(L1,L2,L3)(3*mu1+1)*L2*L3;
d1d1_9=@(L1,L2,L3)(-3*mu2-1)*L2*L3;
%
%               
d1d2_1=@(L1,L2,L3)0;
d1d2_2=@(L1,L2,L3)0;
d1d2_3=@(L1,L2,L3)0;
d1d2_4=@(L1,L2,L3)1;
d1d2_5=@(L1,L2,L3)0;
d1d2_6=@(L1,L2,L3)0;
d1d2_7=@(L1,L2,L3)( 0.5*L3*((3*mu3+1)*L3-(3*mu3+1)*L2+(1-mu3)*L1)...
                   +0.5*(-3*mu3-1)*L2*L3+0.5*(1-mu3)*L1*L3+2*L1);
d1d2_8=@(L1,L2,L3)( 0.5*L3*((-3*mu1-1)*L3+(1-mu1)*L2+(3*mu1+1)*L1)...
                   +0.5*(1-mu1)*L2*L3+0.5*(3*mu1+1)*L1*L3);
d1d2_9=@(L1,L2,L3)( 0.5*L3*((1-mu2)*L3+(3*mu2+1)*L2+(-3*mu2-1)*L1)...
                   +0.5*(3*mu2+1)*L2*L3+0.5*(-3*mu2-1)*L1*L3);
%
% 
d1d3_1=@(L1,L2,L3)0;
d1d3_2=@(L1,L2,L3)0;  
d1d3_3=@(L1,L2,L3)0;  
d1d3_4=@(L1,L2,L3)0;  
d1d3_5=@(L1,L2,L3)0;  
d1d3_6=@(L1,L2,L3)1;  
d1d3_7=@(L1,L2,L3)( 0.5*L2*((3*mu3+1)*L3-(3*mu3+1)*L2+(1-mu3)*L1)...
                   +0.5*(3*mu3+1)*L2*L3+0.5*(1-mu3)*L1*L2);  
d1d3_8=@(L1,L2,L3)( 0.5*L2*((-3*mu1-1)*L3+(1-mu1)*L2+(3*mu1+1)*L1)...
                   +0.5*(-3*mu1-1)*L2*L3+0.5*(3*mu1+1)*L1*L2); 
d1d3_9=@(L1,L2,L3)( 0.5*L2*((1-mu2)*L3+(3*mu2+1)*L2+(-3*mu2-1)*L1)...
                   +0.5*(1-mu2)*L2*L3+2*L3+0.5*(-3*mu2-1)*L1*L2);             
%
d2d2_1=@(L1,L2,L3)0;
d2d2_2=@(L1,L2,L3)0;
d2d2_3=@(L1,L2,L3)0;
d2d2_4=@(L1,L2,L3)0;
d2d2_5=@(L1,L2,L3)0;
d2d2_6=@(L1,L2,L3)0;
d2d2_7=@(L1,L2,L3)(-3*mu3-1)*L1*L3;
d2d2_8=@(L1,L2,L3)(1-mu1)*L1*L3+2*L3;
d2d2_9=@(L1,L2,L3)(3*mu2+1)*L1*L3;
%
d2d3_1=@(L1,L2,L3)0; 
d2d3_2=@(L1,L2,L3)0;  
d2d3_3=@(L1,L2,L3)0;  
d2d3_4=@(L1,L2,L3)0;  
d2d3_5=@(L1,L2,L3)1;  
d2d3_6=@(L1,L2,L3)0;  
d2d3_7=@(L1,L2,L3)(0.5*L1*((3*mu3+1)*L3-(3*mu3+1)*L2+(1-mu3)*L1)...
                  +0.5*(3*mu3+1)*L1*L3+0.5*(-3*mu3-1)*L1*L2);  
d2d3_8=@(L1,L2,L3)(0.5*L1*((-3*mu1-1)*L3+(1-mu1)*L2+(3*mu1+1)*L1)...
                  +0.5*(-3*mu1-1)*L1*L3+0.5*(1-mu1)*L1*L2+2*L2);  
d2d3_9=@(L1,L2,L3)(0.5*L1*((1-mu2)*L3+(3*mu2+1)*L2+(-3*mu2-1)*L1)...
                  +0.5*(1-mu2)*L1*L3+0.5*(3*mu2+1)*L1*L2);              
% 
d3d3_1=@(L1,L2,L3)0;
d3d3_2=@(L1,L2,L3)0;
d3d3_3=@(L1,L2,L3)0;
d3d3_4=@(L1,L2,L3)0;
d3d3_5=@(L1,L2,L3)0;
d3d3_6=@(L1,L2,L3)0;
d3d3_7=@(L1,L2,L3)(3*mu3+1)*L1*L2;
d3d3_8=@(L1,L2,L3)(-3*mu1-1)*L1*L2;
d3d3_9=@(L1,L2,L3)(1-mu2)*L1*L2+2*L1;
%
%
FA2=4*area^2;
%
DXDX_1=@(L1,L2,L3)(     b1^2*d1d1_1(L1,L2,L3)... 
                   + 2*b1*b2*d1d2_1(L1,L2,L3)...
                   + 2*b1*b3*d1d3_1(L1,L2,L3)...
                   +    b2^2*d2d2_1(L1,L2,L3)... 
                   + 2*b2*b3*d2d3_1(L1,L2,L3)...
                   +    b3^2*d3d3_1(L1,L2,L3))/FA2;
%
DXDX_2=@(L1,L2,L3)(     b1^2*d1d1_2(L1,L2,L3)... 
                   + 2*b1*b2*d1d2_2(L1,L2,L3)...
                   + 2*b1*b3*d1d3_2(L1,L2,L3)...
                   +    b2^2*d2d2_2(L1,L2,L3)... 
                   + 2*b2*b3*d2d3_2(L1,L2,L3)...
                   +    b3^2*d3d3_2(L1,L2,L3))/FA2;
%
DXDX_3=@(L1,L2,L3)(     b1^2*d1d1_3(L1,L2,L3)... 
                   + 2*b1*b2*d1d2_3(L1,L2,L3)...
                   + 2*b1*b3*d1d3_3(L1,L2,L3)...
                   +    b2^2*d2d2_3(L1,L2,L3)... 
                   + 2*b2*b3*d2d3_3(L1,L2,L3)...
                   +    b3^2*d3d3_3(L1,L2,L3))/FA2;
%
DXDX_4=@(L1,L2,L3)(     b1^2*d1d1_4(L1,L2,L3)... 
                   + 2*b1*b2*d1d2_4(L1,L2,L3)...
                   + 2*b1*b3*d1d3_4(L1,L2,L3)...
                   +    b2^2*d2d2_4(L1,L2,L3)... 
                   + 2*b2*b3*d2d3_4(L1,L2,L3)...
                   +    b3^2*d3d3_4(L1,L2,L3))/FA2;
%
DXDX_5=@(L1,L2,L3)(     b1^2*d1d1_5(L1,L2,L3)... 
                   + 2*b1*b2*d1d2_5(L1,L2,L3)...
                   + 2*b1*b3*d1d3_5(L1,L2,L3)...
                   +    b2^2*d2d2_5(L1,L2,L3)... 
                   + 2*b2*b3*d2d3_5(L1,L2,L3)...
                   +    b3^2*d3d3_5(L1,L2,L3))/FA2;
%
DXDX_6=@(L1,L2,L3)(     b1^2*d1d1_6(L1,L2,L3)... 
                   + 2*b1*b2*d1d2_6(L1,L2,L3)...
                   + 2*b1*b3*d1d3_6(L1,L2,L3)...
                   +    b2^2*d2d2_6(L1,L2,L3)... 
                   + 2*b2*b3*d2d3_6(L1,L2,L3)...
                   +    b3^2*d3d3_6(L1,L2,L3))/FA2;
%
DXDX_7=@(L1,L2,L3)(     b1^2*d1d1_7(L1,L2,L3)... 
                   + 2*b1*b2*d1d2_7(L1,L2,L3)...
                   + 2*b1*b3*d1d3_7(L1,L2,L3)...
                   +    b2^2*d2d2_7(L1,L2,L3)... 
                   + 2*b2*b3*d2d3_7(L1,L2,L3)...
                   +    b3^2*d3d3_7(L1,L2,L3))/FA2;
%
DXDX_8=@(L1,L2,L3)(     b1^2*d1d1_8(L1,L2,L3)... 
                   + 2*b1*b2*d1d2_8(L1,L2,L3)...
                   + 2*b1*b3*d1d3_8(L1,L2,L3)...
                   +    b2^2*d2d2_8(L1,L2,L3)... 
                   + 2*b2*b3*d2d3_8(L1,L2,L3)...
                   +    b3^2*d3d3_8(L1,L2,L3))/FA2;
%
DXDX_9=@(L1,L2,L3)(     b1^2*d1d1_9(L1,L2,L3)... 
                   + 2*b1*b2*d1d2_9(L1,L2,L3)...
                   + 2*b1*b3*d1d3_9(L1,L2,L3)...
                   +    b2^2*d2d2_9(L1,L2,L3)... 
                   + 2*b2*b3*d2d3_9(L1,L2,L3)...
                   +    b3^2*d3d3_9(L1,L2,L3))/FA2;
%%%
%
DYDY_1=@(L1,L2,L3)(     c1^2*d1d1_1(L1,L2,L3)... 
                   + 2*c1*c2*d1d2_1(L1,L2,L3)...
                   + 2*c1*c3*d1d3_1(L1,L2,L3)...
                   +    c2^2*d2d2_1(L1,L2,L3)... 
                   + 2*c2*c3*d2d3_1(L1,L2,L3)...
                   +    c3^2*d3d3_1(L1,L2,L3))/FA2;
%
DYDY_2=@(L1,L2,L3)(     c1^2*d1d1_2(L1,L2,L3)... 
                   + 2*c1*c2*d1d2_2(L1,L2,L3)...
                   + 2*c1*c3*d1d3_2(L1,L2,L3)...
                   +    c2^2*d2d2_2(L1,L2,L3)... 
                   + 2*c2*c3*d2d3_2(L1,L2,L3)...
                   +    c3^2*d3d3_2(L1,L2,L3))/FA2;
%
DYDY_3=@(L1,L2,L3)(     c1^2*d1d1_3(L1,L2,L3)... 
                   + 2*c1*c2*d1d2_3(L1,L2,L3)...
                   + 2*c1*c3*d1d3_3(L1,L2,L3)...
                   +    c2^2*d2d2_3(L1,L2,L3)... 
                   + 2*c2*c3*d2d3_3(L1,L2,L3)...
                   +    c3^2*d3d3_3(L1,L2,L3))/FA2;
%
DYDY_4=@(L1,L2,L3)(     c1^2*d1d1_4(L1,L2,L3)... 
                   + 2*c1*c2*d1d2_4(L1,L2,L3)...
                   + 2*c1*c3*d1d3_4(L1,L2,L3)...
                   +    c2^2*d2d2_4(L1,L2,L3)... 
                   + 2*c2*c3*d2d3_4(L1,L2,L3)...
                   +    c3^2*d3d3_4(L1,L2,L3))/FA2;
%
DYDY_5=@(L1,L2,L3)(     c1^2*d1d1_5(L1,L2,L3)... 
                   + 2*c1*c2*d1d2_5(L1,L2,L3)...
                   + 2*c1*c3*d1d3_5(L1,L2,L3)...
                   +    c2^2*d2d2_5(L1,L2,L3)... 
                   + 2*c2*c3*d2d3_5(L1,L2,L3)...
                   +    c3^2*d3d3_5(L1,L2,L3))/FA2;
%
DYDY_6=@(L1,L2,L3)(     c1^2*d1d1_6(L1,L2,L3)... 
                   + 2*c1*c2*d1d2_6(L1,L2,L3)...
                   + 2*c1*c3*d1d3_6(L1,L2,L3)...
                   +    c2^2*d2d2_6(L1,L2,L3)... 
                   + 2*c2*c3*d2d3_6(L1,L2,L3)...
                   +    c3^2*d3d3_6(L1,L2,L3))/FA2;
%
DYDY_7=@(L1,L2,L3)(     c1^2*d1d1_7(L1,L2,L3)... 
                   + 2*c1*c2*d1d2_7(L1,L2,L3)...
                   + 2*c1*c3*d1d3_7(L1,L2,L3)...
                   +    c2^2*d2d2_7(L1,L2,L3)... 
                   + 2*c2*c3*d2d3_7(L1,L2,L3)...
                   +    c3^2*d3d3_7(L1,L2,L3))/FA2;
%
DYDY_8=@(L1,L2,L3)(     c1^2*d1d1_8(L1,L2,L3)... 
                   + 2*c1*c2*d1d2_8(L1,L2,L3)...
                   + 2*c1*c3*d1d3_8(L1,L2,L3)...
                   +    c2^2*d2d2_8(L1,L2,L3)... 
                   + 2*c2*c3*d2d3_8(L1,L2,L3)...
                   +    c3^2*d3d3_8(L1,L2,L3))/FA2;
%
DYDY_9=@(L1,L2,L3)(     c1^2*d1d1_9(L1,L2,L3)... 
                   + 2*c1*c2*d1d2_9(L1,L2,L3)...
                   + 2*c1*c3*d1d3_9(L1,L2,L3)...
                   +    c2^2*d2d2_9(L1,L2,L3)... 
                   + 2*c2*c3*d2d3_9(L1,L2,L3)...
                   +    c3^2*d3d3_9(L1,L2,L3))/FA2;
%%
DXDY_1=@(L1,L2,L3)(          b1*c1*d1d1_1(L1,L2,L3)... 
                   + (b1*c2+b2*c1)*d1d2_1(L1,L2,L3)...
                   + (b1*c3+b3*c1)*d1d3_1(L1,L2,L3)...
                   +         b2*c2*d2d2_1(L1,L2,L3)... 
                   + (b2*c3+b3*c2)*d2d3_1(L1,L2,L3)...
                   +         b3*c3*d3d3_1(L1,L2,L3))/FA2;
%
DXDY_2=@(L1,L2,L3)(          b1*c1*d1d1_2(L1,L2,L3)... 
                   + (b1*c2+b2*c1)*d1d2_2(L1,L2,L3)...
                   + (b1*c3+b3*c1)*d1d3_2(L1,L2,L3)...
                   +         b2*c2*d2d2_2(L1,L2,L3)... 
                   + (b2*c3+b3*c2)*d2d3_2(L1,L2,L3)...
                   +         b3*c3*d3d3_2(L1,L2,L3))/FA2;
%
DXDY_3=@(L1,L2,L3)(          b1*c1*d1d1_3(L1,L2,L3)... 
                   + (b1*c2+b2*c1)*d1d2_3(L1,L2,L3)...
                   + (b1*c3+b3*c1)*d1d3_3(L1,L2,L3)...
                   +         b2*c2*d2d2_3(L1,L2,L3)... 
                   + (b2*c3+b3*c2)*d2d3_3(L1,L2,L3)...
                   +         b3*c3*d3d3_3(L1,L2,L3))/FA2;
%
DXDY_4=@(L1,L2,L3)(          b1*c1*d1d1_4(L1,L2,L3)... 
                   + (b1*c2+b2*c1)*d1d2_4(L1,L2,L3)...
                   + (b1*c3+b3*c1)*d1d3_4(L1,L2,L3)...
                   +         b2*c2*d2d2_4(L1,L2,L3)... 
                   + (b2*c3+b3*c2)*d2d3_4(L1,L2,L3)...
                   +         b3*c3*d3d3_4(L1,L2,L3))/FA2;
%
DXDY_5=@(L1,L2,L3)(          b1*c1*d1d1_5(L1,L2,L3)... 
                   + (b1*c2+b2*c1)*d1d2_5(L1,L2,L3)...
                   + (b1*c3+b3*c1)*d1d3_5(L1,L2,L3)...
                   +         b2*c2*d2d2_5(L1,L2,L3)... 
                   + (b2*c3+b3*c2)*d2d3_5(L1,L2,L3)...
                   +         b3*c3*d3d3_5(L1,L2,L3))/FA2;
%
DXDY_6=@(L1,L2,L3)(          b1*c1*d1d1_6(L1,L2,L3)... 
                   + (b1*c2+b2*c1)*d1d2_6(L1,L2,L3)...
                   + (b1*c3+b3*c1)*d1d3_6(L1,L2,L3)...
                   +         b2*c2*d2d2_6(L1,L2,L3)... 
                   + (b2*c3+b3*c2)*d2d3_6(L1,L2,L3)...
                   +         b3*c3*d3d3_6(L1,L2,L3))/FA2;
%
DXDY_7=@(L1,L2,L3)(          b1*c1*d1d1_7(L1,L2,L3)... 
                   + (b1*c2+b2*c1)*d1d2_7(L1,L2,L3)...
                   + (b1*c3+b3*c1)*d1d3_7(L1,L2,L3)...
                   +         b2*c2*d2d2_7(L1,L2,L3)... 
                   + (b2*c3+b3*c2)*d2d3_7(L1,L2,L3)...
                   +         b3*c3*d3d3_7(L1,L2,L3))/FA2;
%
DXDY_8=@(L1,L2,L3)(          b1*c1*d1d1_8(L1,L2,L3)... 
                   + (b1*c2+b2*c1)*d1d2_8(L1,L2,L3)...
                   + (b1*c3+b3*c1)*d1d3_8(L1,L2,L3)...
                   +         b2*c2*d2d2_8(L1,L2,L3)... 
                   + (b2*c3+b3*c2)*d2d3_8(L1,L2,L3)...
                   +         b3*c3*d3d3_8(L1,L2,L3))/FA2;
%
DXDY_9=@(L1,L2,L3)(          b1*c1*d1d1_9(L1,L2,L3)... 
                   + (b1*c2+b2*c1)*d1d2_9(L1,L2,L3)...
                   + (b1*c3+b3*c1)*d1d3_9(L1,L2,L3)...
                   +         b2*c2*d2d2_9(L1,L2,L3)... 
                   + (b2*c3+b3*c2)*d2d3_9(L1,L2,L3)...
                   +         b3*c3*d3d3_9(L1,L2,L3))/FA2;
%
%%%
W=1/3;
%
LN=[ 4 1 1; 1 4 1; 1 1 4]/6;
%
asum=zeros(9,9);
%
for i=1:3
%
  L1=LN(i,1);
  L2=LN(i,2);
  L3=LN(i,3);
%
 DXM=[ DXDX_1(L1,L2,L3)...  
       DXDX_2(L1,L2,L3)...
       DXDX_3(L1,L2,L3)...
       DXDX_4(L1,L2,L3)...
       DXDX_5(L1,L2,L3)...       
       DXDX_6(L1,L2,L3)...
       DXDX_7(L1,L2,L3)...
       DXDX_8(L1,L2,L3)...
       DXDX_9(L1,L2,L3) ];
%
 DYM=[ DYDY_1(L1,L2,L3)...  
       DYDY_2(L1,L2,L3)...
       DYDY_3(L1,L2,L3)...
       DYDY_4(L1,L2,L3)...
       DYDY_5(L1,L2,L3)...       
       DYDY_6(L1,L2,L3)...
       DYDY_7(L1,L2,L3)...
       DYDY_8(L1,L2,L3)...
       DYDY_9(L1,L2,L3) ];
%
DXYM=[ DXDY_1(L1,L2,L3)...  
       DXDY_2(L1,L2,L3)...
       DXDY_3(L1,L2,L3)...
       DXDY_4(L1,L2,L3)...
       DXDY_5(L1,L2,L3)...       
       DXDY_6(L1,L2,L3)...
       DXDY_7(L1,L2,L3)...
       DXDY_8(L1,L2,L3)...
       DXDY_9(L1,L2,L3) ];
%
    V=[ DXM ; DYM; 2*DXYM ];
%
    asum=asum+V'*D*V;
%
end
%
k=asum*W*area;
%
%
[cbar]=triangular_specht_cbar(mu1,mu2,mu3,b1,b2,b3,c1,c2,c3,TA);
%
ml=cbar'*m*cbar;
kl=cbar'*k*cbar;
%
for i=1:9
    for j=i+1:9
        ml(j,i)=ml(i,j);
        kl(j,i)=kl(i,j);        
    end
end