function M=matrix_Mass(Aux,C,C1,C2)
global Invariant ini_Val Material dt

if isfield(ini_Val,'iopen')
    iopen = ini_Val.iopen;
else
    iopen = 0;
end
 
    

M=zeros(ini_Val.ncontrol);
connect=Invariant.connect;
aux = Aux.au .* C(:,1) .* 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

% if (iopen == 3)
%     M(:,end) = 0;
%     M(end,:) = 0;
% end

M = M * Material.eta_a0;




