function F=vector_rhs_psi(psi_1,phi)
global Invariant ini_Val Material dt

aux0 = phi/Material.phi_bar;
aux0(aux0>3) = 3;
aux = psi_1/dt + Material.k1 * subplus(aux0 - Material.enrichment_treshold);
aux = aux .* Invariant.gweight.*Invariant.gaussp*Material.R0;


F=zeros(ini_Val.ncontrol,1);
connect=Invariant.connect;


for iel=1:size(connect,1)
    ig=(1:ini_Val.gaussno);
    igglob=(iel-1)*ini_Val.gaussno+ig;
    F_loc= (aux(igglob))*Invariant.Nshape(igglob,:);
    F(connect(iel,:)) = ...
        F(connect(iel,:)) + F_loc(:);
end
