% tDL2 is for L2 dissipation

function [tD,dD_dCP] = Incremental_mem_diss(Aux,C,C_1,C1)
global Invariant ini_Val Material phi dt 

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

p = ini_Val.p;
gaussno = ini_Val.gaussno;

dr = C(:,1)-C_1(:,1);
dz = C(:,2)-C_1(:,2);

if (xi ~=0)
    
    tDL2 = sum((-dr .*C1(:,2) + dz.* C1(:,1)).^2./Aux.au.*C(:,1).* Invariant.gweight);
    %tDL2 = tmp* pi* xi/dt;
else
    tDL2=0;
end

tD = sum( dr.^2.*Aux.au./C(:,1).* Invariant.gweight);
tD = tD * 4*pi/dt*Material.mu + tDL2* pi/dt * xi;

if nargout == 2
    %Derivatives of E (local)
    temp1 = (dr.*Aux.au./C(:,1).*(2-dr./C(:,1)))*ones(1,p+1);
    temp2 = (dr.^2./C(:,1))*ones(1,p+1);
    
    dD_dr_l = (temp1.* Invariant.Nshape + temp2.*Aux.da_dr_l   ).*Aux.Mgweight;
    dD_dz_l = (temp2.*Aux.da_dz_l   ).*Aux.Mgweight;    
    ns = length(dD_dr_l)/gaussno;
    vs = [1:1:ns];
    dD_dCPr = reshape(sum(reshape(dD_dr_l,gaussno,ns,p+1),1),ns,p+1);
    dD_dCPz = reshape(sum(reshape(dD_dz_l,gaussno,ns,p+1),1),ns,p+1);
    dD_dCP = zeros(2,ini_Val.ncontrol);

    for i = 0:p
        dD_dCP(1,vs(:)+i) = dD_dCP(1,vs(:)+i) + dD_dCPr(:,i+1)';
        dD_dCP(2,vs(:)+i) = dD_dCP(2,vs(:)+i) + dD_dCPz(:,i+1)';
    end
    
    if (xi ~=0)
        tempA = (- C1(:,2) .* dr + C1(:,1).*dz).^2*ones(1,p+1);
        tempB = 2*(-dr.*C1(:,2) + dz.*C1(:,1)).*(C(:,1)./Aux.au)*ones(1,p+1);
        
        tempC = (1./Aux.au)*ones(1,p+1);
        tempD = (C(:,1)./Aux.au.^2)*ones(1,p+1);
        tempE = C1(:,1)*ones(1,p+1);
        tempF = C1(:,2)*ones(1,p+1);
        
        
        dD_dr_l2 = ((tempC.*Invariant.Nshape - tempD.*Aux.da_dr_l) .* tempA ...
            + tempB.*(-tempF.*Invariant.Nshape + dz*ones(1,p+1).* Invariant.dNsh)).*Aux.Mgweight;
        dD_dz_l2 = (- tempD.*Aux.da_dz_l.* tempA ...
            + tempB.*(-dr*ones(1,p+1).*Invariant.dNsh + Invariant.Nshape.*tempE)).*Aux.Mgweight;
        
        dD_dCPr2 = reshape(sum(reshape(dD_dr_l2,gaussno,ns,p+1),1),ns,p+1);
        dD_dCPz2 = reshape(sum(reshape(dD_dz_l2,gaussno,ns,p+1),1),ns,p+1);
        dD_dCP2 = zeros(2,ini_Val.ncontrol);
        
        for i = 0:p
            dD_dCP2(1,vs(:)+i) = dD_dCP2(1,vs(:)+i) + dD_dCPr2(:,i+1)';
            dD_dCP2(2,vs(:)+i) = dD_dCP2(2,vs(:)+i) + dD_dCPz2(:,i+1)';
        end
    end
        

    dD_dCP= dD_dCP * 4*pi/dt*Material.mu;
    if (xi ~=0)
        dD_dCP = dD_dCP + (dD_dCP2 * pi/dt*xi);
    end
end
