%% 质量矩阵组装 Mass matrix assembly
% {M_{\sigma } }_{ij}^e ={\int ~}_{\Omega ~e} N_i^e \cdot \sigma^e {\cdot N}_j^e dV
% The assembly does not take i\omega \mu_0 into account
% Jin, 2014.  8.3.4 equation 8.68

function M = get_M_v2(abcd,vol,elems2edges,edgeLength,elemEdgeLength,sign,sig)

b=abcd(:,5:8);
c=abcd(:,9:12);
d =abcd(:,13:16);
neles    = size(elems2edges,1);
nedges   = size(edgeLength,1);
M_ij     = zeros(36*neles,1);
ii       = zeros(36*neles,1);
jj       = zeros(36*neles,1);

f = get_f(b,c,d,sig,neles) ;

index   = 0;
for ei= 1:6
    for ej = 1:6
        
        ii(index+1:index+neles) = elems2edges(:,ei);
        jj(index+1:index+neles) = elems2edges(:,ej);
        length_i = elemEdgeLength(:,ei);   %第i边的边长 The length of ith edge
        length_j = elemEdgeLength(:,ej);   %第j边边长 The length of jth edge
        Sign_ij = sign(:,ei).*sign(:,ej);
        
        if ei==1 && ej==1
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,2,2)-2*f(:,1,2)+f(:,1,1))./(360*vol);
            
        elseif (ei==1 && ej==2) || (ei==2 && ej==1) 
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(2*f(:,2,3)-f(:,2,1)-f(:,1,3)+f(:,1,1))./(720*vol);
            
        elseif (ei==1 && ej==3) || (ei==3 && ej==1)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(2*f(:,2,4)-f(:,2,1)-f(:,1,4)+f(:,1,1))./(720*vol);
            
        elseif (ei==1 && ej==4) || (ei==4 && ej==1)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,2,3)-f(:,2,2)-2*f(:,1,3)+f(:,1,2))./(720*vol);
            
        elseif (ei==1 && ej==5) || (ei==5 && ej==1)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,2,2)-f(:,2,4)-f(:,1,2)+2*f(:,1,4))./(720*vol);
            
        elseif (ei==1 && ej==6) || (ei==6 && ej==1)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,2,4)-f(:,2,3)-f(:,1,4)+f(:,1,3))./(720*vol);
            
        elseif (ei==2 && ej==2)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,3,3)-f(:,1,3)+f(:,1,1))./(360*vol);
            
        elseif (ei==2 && ej==3) || (ei==3 && ej==2)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(2*f(:,3,4)-f(:,1,3)-f(:,1,4)+f(:,1,1))./(720*vol);
            
        elseif (ei==2 && ej==4) || (ei==4 && ej==2)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,3,3)-f(:,2,3)-f(:,1,3)+2*f(:,1,2))./(720*vol);
            
        elseif (ei==2 && ej==5) || (ei==5 && ej==2)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,2,3)-f(:,3,4)-f(:,1,2)+f(:,1,4))./(720*vol);
            
        elseif (ei==2 && ej==6) || (ei==6 && ej==2)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,1,3)-f(:,3,3)-2*f(:,1,4)+f(:,3,4))./(720*vol);
            
        elseif (ei==3 && ej==3)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,4,4)-f(:,1,4)+f(:,1,1))./(360*vol);
            
        elseif (ei==3 && ej==4) || (ei==4 && ej==3)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,3,4)-f(:,2,4)-f(:,1,3)+f(:,1,2))./(720*vol);
            
        elseif (ei==3 && ej==5) || (ei==5 && ej==3)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,2,4)-f(:,4,4)-2*f(:,1,2)+f(:,1,4))./(720*vol);
            
        elseif (ei==3 && ej==6) || (ei==6 && ej==3)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,4,4)-f(:,3,4)-f(:,1,4)+2*f(:,1,3))./(720*vol);
            
        elseif (ei==4 && ej==4)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,3,3)-f(:,2,3)+f(:,2,2))./(360*vol);
            
        elseif (ei==4 && ej==5) || (ei==5 && ej==4)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,2,3)-2*f(:,3,4)-f(:,2,2)+f(:,2,4))./(720*vol);
            
        elseif (ei==4 && ej==6) || (ei==6 && ej==4)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,3,4)-f(:,3,3)-2*f(:,2,4)+f(:,2,3))./(720*vol);
            
        elseif (ei==5 && ej==5)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,2,2)-f(:,2,4)+f(:,4,4))./(360*vol);
            
        elseif (ei==5 && ej==6) || (ei==6 && ej==5)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,2,4)-2*f(:,2,3)-f(:,4,4)+f(:,3,4))./(720*vol);
            
        elseif (ei==6 && ej==6)
            M_ij(index+1:index+neles ) = (length_i.*length_j).*(f(:,4,4)-f(:,3,4)+f(:,3,3))./(360*vol);
        end

        M_ij(index+1:index+neles) = Sign_ij.*M_ij(index+1:index+neles);
        index = index + neles;
    end
end

M = sparse(ii,jj,M_ij,nedges,nedges);

end

%
% Jin,2014, page 301: fij=bi*bj + ci*cj + di*dj
% Cao Xiaoyue, 2019. PhD Thesis, Jilin University, Page: 38-40
% \begin{array}{c}
% f_{i,j} =b_i b_j \sigma_{xx} +b_i c_j \sigma_{xy} +b_i d_j \sigma_{xz} +\\
% {\ \ \ \ \ \ \ \ \ c_i b}_j \sigma_{yx} +c_i c_j \sigma_{yy} +c_i d_j \sigma_{yz} +\\
% {\ \ \ \ d}_i b_j \sigma_{zx} +d_i c_j \sigma_{zy} +b_i d_j \sigma_{xz}
% \end{array}     ;    \sigma =\left\lbrack \begin{array}{ccc}
% \sigma_{\textrm{xx}}  & \sigma_{\textrm{xy}}  & \sigma_{\textrm{xz}} \\
% \sigma_{\textrm{yx}}  & \sigma_{\textrm{yy}}  & \sigma_{\textrm{yz}} \\
% \sigma_{\textrm{zx}}  & \sigma_{\textrm{zy}}  & \sigma_{\textrm{zz}}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{ccc}
% \sigma_1  & \sigma_4  & \sigma_5 \\
% \sigma_7  & \sigma_2  & \sigma_6 \\
% \sigma_8  & \sigma_9  & \sigma_3
% \end{array}\right\rbrack =\left\lbrack {\textrm{sig}}_1 ,\ldotp \ldotp \ldotp \ldotp {\textrm{sig}}_9 \right\rbrack

function f = get_f(b,c,d,sig,neles)

ani = size(sig,2);
f=zeros(neles,4,4);

if ani==  1  % Isotropic
    for i =1:4
        for j =1:4
            f(:,i,j) =  (b(:,i).*b(:,j) +c(:,i).*c(:,j) +d(:,i).*d(:,j)) .*sig;
        end
    end
elseif  ani== 3 % axis anisotropy, Diagonal elements of the tensor
    for i =1:4
        for j =1:4
            f(:,i,j) = b(:,i).*b(:,j).*sig(:,1) + c(:,i).*c(:,j).*sig(:,2) + d(:,i).*d(:,j).*sig(:,3) ;
        end
    end
elseif  ani== 6 % Arbitrary anisotropy, symmetric tensor
    for i =1:4
        for j =1:4
            f(:,i,j) = ...
                b(:,i).*b(:,j).*sig(:,1) + b(:,i).*c(:,j).*sig(:,4) + b(:,i).*d(:,j).*sig(:,5)+...
                c(:,i).*b(:,j).*sig(:,4) + c(:,i).*c(:,j).*sig(:,2) + c(:,i).*d(:,j).*sig(:,6)+...
                d(:,i).*b(:,j).*sig(:,5) + d(:,i).*c(:,j).*sig(:,6) + d(:,i).*d(:,j).*sig(:,3)  ;
        end
    end
elseif  ani==9 % tensor
    for i =1:4
        for j =1:4
            f(:,i,j) = ...
                b(:,i).*b(:,j).*sig(:,1) + b(:,i).*c(:,j).*sig(:,4) + b(:,i).*d(:,j).*sig(:,5)+...
                c(:,i).*b(:,j).*sig(:,7) + c(:,i).*c(:,j).*sig(:,2) + c(:,i).*d(:,j).*sig(:,6)+...
                d(:,i).*b(:,j).*sig(:,8) + d(:,i).*c(:,j).*sig(:,9) + d(:,i).*d(:,j).*sig(:,3)  ;
        end
    end
    
end
end
