function K=matrix_Stiffness_Cortex
global Invariant ini_Val Material


K=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;
    K_loc=Invariant.dNsh(igglob,:)'*diag(aux(igglob))*Invariant.dNsh(igglob,:);
    K(connect(iel,:),connect(iel,:)) = ...
        K(connect(iel,:),connect(iel,:)) + K_loc;
end

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