function R = ResidualP2(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

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

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


% 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);

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

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

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

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

aux4 =  tmp2.* (Material.RT_a0 * Material.phi_max./(Material.phi_max - phi) ...
    + (Material.k0*Constraints.C0^2 + Material.Chi).*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 (iopen == 3 || iopen ==6 || iopen ==7)
    %R(1)=0;
    R(end) = 0;
    R(end-1)=0;
    %R(end-2)=0;
end

    
%     aux5 = tmp3.^2./tmp1.* d_phi *(Lmd) .*Invariant.gweight;
%     aux6 = tmp2.* tmp3./tmp1.* d_phi * (Lmd) .* Invariant.gweight;
%     aux7 = tmp2.* tmp3./tmp1 .* dd_phi * (Lmd) .* Invariant.gweight;
%     aux8 = tmp2.^2./tmp1 .* dd_phi * (Lmd) .* Invariant.gweight;