function M=matrix_Mass_Cortex
global Invariant ini_Val Material

    

M=zeros(ini_Val.ncontrol);
connect=Invariant.connect;
aux = 2./Invariant.gaussp.*Invariant.gweight;

for iel=1:size(connect,1)
    ig=(1:ini_Val.gaussno);
    igglob=(iel-1)*ini_Val.gaussno+ig;
    M_loc=Invariant.Nshape(igglob,:)'*diag(aux(igglob))*Invariant.Nshape(igglob,:);
    M(connect(iel,:),connect(iel,:)) = ...
        M(connect(iel,:),connect(iel,:)) + M_loc;
end


M = M * Material.eta_c/Material.R0;




