
function R = Residual_mod_Ener_Aug(Aux,C,C1,C2,phi,d_phi,phi_1,dd_phi)
global Invariant ini_Val Material Constraints dt

if isfield(Material,'Lmd')
    Lmd = Material.Lmd;
else
    Lmd = 0;
end

% For the term added to the energy
if isfield(Material,'alpha')
    alpha = Material.alpha;
else
    alpha = 0;
end

if isfield(Material,'beta')
    beta = Material.beta;
else
    beta = 0;
end

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

if isfield(Constraints,'FixChem')
    F_Chem = Constraints.FixChem;
else
    F_Chem = 0;
end

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


Chi = Material.Chi;
RT_a0 = Material.RT_a0;
C0 = Constraints.C0;

% To Compute R1, R2, R3 and R4
tmp1 = (C(:,1).* Aux.au);
tmp2 = (C(:,1)./ Aux.au);
tmp3 = ((C1(:,1).*Aux.au - C(:,1).*Aux.dau)./(Aux.au).^2);

tmp4 = beta./(phi)./(phi) + alpha./(Material.phi_max-phi).^2;
% tmp4 = alpha*beta^2*sech(beta*(Material.phi_max-phi)).^2.*tanh(beta*(Material.phi_max-phi));
% tmp6 = -gamma*delta^2.*sech(delta*(phi+0.05)).^2.*tanh(delta*(phi+0.05));

aux1 = tmp1.* (Material.eta_a0*phi) .* Invariant.gweight;

aux2 = - tmp1.* (Material.eta_a0*phi_1) .* Invariant.gweight;

aux3_a = tmp3.*(Material.k0*C0.*phi.*Aux.Hu).*Invariant.gweight*dt;
aux3_b = tmp2.*(Material.k0*C0.*d_phi.*Aux.Hu).*Invariant.gweight*dt;
aux3_c = tmp2.*(Material.k0*C0.*phi.*Aux.Hu).*Invariant.gweight*dt;

%aux3 = - tmp2.* (Material.k0*Constraints.C0*phi.*Aux.Hpu).*Invariant.gweight*dt;

aux4 =  tmp2.* (RT_a0 * Material.phi_max./(Material.phi_max - phi) ...
    + (Material.k0*C0.^2 + Chi+ tmp4).*phi).* d_phi .*Invariant.gweight*dt;


if (nargin == 8)

    L_phi_m = (tmp3.* d_phi + (C(:,1)./ Aux.au).*dd_phi)./tmp1;
    
    aux5 = tmp3.*phi.*L_phi_m.*Lmd .*Invariant.gweight*dt;
    aux6 = tmp2.*d_phi.*L_phi_m.*Lmd .*Invariant.gweight*dt;
    aux7 = tmp2.*phi.*L_phi_m.*Lmd .*Invariant.gweight*dt;
end


for iel=1:size(connect,1)
    ig=(1:ini_Val.gaussno);
    igglob=(iel-1)*ini_Val.gaussno+ig;
    R1_loc = (aux1(igglob)')*Invariant.Nshape(igglob,:);
    R2_loc = (aux2(igglob)')*Invariant.Nshape(igglob,:);
    R3_loc_a = (aux3_a(igglob)')*Invariant.dNsh(igglob,:);
    R3_loc_b = (aux3_b(igglob)')*Invariant.dNsh(igglob,:);
    R3_loc_c = (aux3_c(igglob)')*Invariant.ddNsh(igglob,:);
    R4_loc = (aux4(igglob)')*Invariant.dNsh(igglob,:);
    R_loc = R1_loc + R2_loc + R3_loc_a +R3_loc_b+ R3_loc_c + R4_loc;
    if (nargin == 8)
        R5_loc = (aux5(igglob)')*Invariant.dNsh(igglob,:);
        R6_loc = (aux6(igglob)')*Invariant.dNsh(igglob,:);
        R7_loc = (aux7(igglob)')*Invariant.ddNsh(igglob,:);
        R_loc = R_loc + R5_loc + R6_loc + R7_loc;
    end                
    R(connect(iel,:)) = ...
        R(connect(iel,:)) + R_loc(:);
end

if (F_Chem ==1)
    R(end-2) = 0;
    R(end-1) = 0;
    R(end) = 0;
end

