function JD = matrix_JD(Aux,C,C1,C2,phi,d_phi,dd_phi)
%JD: Jacobian Dirichlet

global Invariant ini_Val Material Constraints dt

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

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


JD=zeros(ini_Val.ncontrol);

KD=zeros(ini_Val.ncontrol);
AD=zeros(ini_Val.ncontrol);
BD=zeros(ini_Val.ncontrol);
CD=zeros(ini_Val.ncontrol);
DD=zeros(ini_Val.ncontrol);
ED=zeros(ini_Val.ncontrol);

connect = Invariant.connect;
b = (C1(:,1).*Aux.au-C(:,1).*Aux.dau)./Aux.au.^2;
L_phi = (b.* d_phi + (C(:,1)./ Aux.au).*dd_phi);

tmp = (1./ Aux.au).^2;


aux_K1 =  tmp .* L_phi .* Lmd.* Invariant.gweight*dt;

aux_K2 =  tmp.* d_phi.* b.* Lmd.* Invariant.gweight*dt;

aux_K3 =  b.^2 .* phi ./ C(:,1) ./ Aux.au .* Lmd.* Invariant.gweight*dt;

aux_A1 =  C(:,1)./Aux.au.*tmp.* d_phi.* Lmd.* Invariant.gweight*dt;

aux_A2 = b .* phi .* tmp .* Lmd.* Invariant.gweight*dt;

aux_B =  b./ C(:,1) ./ Aux.au.* L_phi .* Lmd.* Invariant.gweight*dt;

aux_C =  tmp .* L_phi .* Lmd.* Invariant.gweight*dt;

aux_D =  tmp .* phi .* b .* Lmd.* Invariant.gweight*dt;

aux_E =  C(:,1)./Aux.au.*tmp .* phi .* Lmd.* Invariant.gweight*dt;


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_K1(igglob)+aux_K2(igglob)+aux_K3(igglob))*Invariant.dNsh(igglob,:);
    KD(connect(iel,:),connect(iel,:)) = ...
        KD(connect(iel,:),connect(iel,:)) + K_loc;
    
    A_loc=Invariant.ddNsh(igglob,:)'*diag(aux_A1(igglob)+aux_A2(igglob))*Invariant.dNsh(igglob,:);
    AD(connect(iel,:),connect(iel,:)) = ...
        AD(connect(iel,:),connect(iel,:)) + A_loc;
    
    B_loc=Invariant.Nshape(igglob,:)'*diag(aux_B(igglob))*Invariant.dNsh(igglob,:);
    BD(connect(iel,:),connect(iel,:)) = ...
        BD(connect(iel,:),connect(iel,:)) + B_loc;
    
    C_loc=Invariant.Nshape(igglob,:)'*diag(aux_C(igglob))*Invariant.ddNsh(igglob,:);
    CD(connect(iel,:),connect(iel,:)) = ...
        CD(connect(iel,:),connect(iel,:)) + C_loc;
    
    D_loc=Invariant.dNsh(igglob,:)'*diag(aux_D(igglob))*Invariant.ddNsh(igglob,:);
    DD(connect(iel,:),connect(iel,:)) = ...
        DD(connect(iel,:),connect(iel,:)) + D_loc;
    
    
    E_loc=Invariant.ddNsh(igglob,:)'*diag(aux_E(igglob))*Invariant.ddNsh(igglob,:);
    ED(connect(iel,:),connect(iel,:)) = ...
        ED(connect(iel,:),connect(iel,:)) + E_loc;
end

JD = AD + BD + CD + DD + ED + KD;
% 
% if (iopen == 3)
%     JD(:,end) = 0;
%     JD(end,:) = 0;
% end
