function [Dtan,stress,stress_vmisses,strain,strain_vmisses,u,stress_trial] =...
    RVE_response(vol,numelem,element,node,u,TDOF,act,material_type,Indexer,...
    stress,macro_F,strain,stress_vmisses,strain_vmisses,W,Q,elemType,...
    StVenant_Kirchhoff,NeoHookean_Compressible,NeoHookean_Incompressible,Ferrara_Pandolfi,...
    I,L,Lz,L1,IdyadI,max_iter,global_toll)

%--------------------------------------------------------------------------
%initialization begins
%--------------------------------------------------------------------------
Du=zeros(TDOF,1);
K=zeros((3*8)^2,numelem);%Exclusive to H4 element

sctB    = zeros(3*8,numelem) ;
sctB_r  = zeros((3*8)^2,numelem);
sctB_c  = zeros((3*8)^2,numelem);

for iel=1:numelem

    sctr=element(iel,:);
    n=length(sctr);
    sctrB=zeros(1,3*n); %Element degrees of freedom
    sctrB(1:3:3*n)=3*sctr-2;
    sctrB(2:3:3*n)=3*sctr-1;
    sctrB(3:3:3*n)=3*sctr;
    elementu=u(sctrB);
    elementu=reshape(elementu,3,n);

    K_e=0;

    material_type_e=material_type(iel,:);

    for igp=1:size(W,1)

        strain_gp=strain(iel,((igp-1)*6)+1:igp*6);
        strain_vmisses_gp=zeros(1,1);
        stress_gp=stress(iel,((igp-1)*6)+1:igp*6);
        stress_vmisses_gp=zeros(1,1);

        [N,dNdxi] = lagrange_basis(elemType,Q(igp,:));
        J=node(sctr,:)'*dNdxi;
        dNdx=dNdxi/J;
        F=macro_F+elementu*dNdx;

        %%-----------------------------------------------------------------------------
        %%calcMaterialResponse begins.
        %%-----------------------------------------------------------------------------
        Dtan=0;

        switch material_type_e(1,1)
            %%*************************************************************************
            case 1 %StVenant_Kirchhoff
                [Dtan,stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,stress_trial]=...
                    StVenant_Kirchhoff_response(StVenant_Kirchhoff,material_type_e,F,...
                    stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,...
                    I,L,L1,IdyadI);
                %%*************************************************************************
            case 2 %Neo_Hookean_Compressible_response
                [Dtan,stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,stress_trial]=...
                    NeoHookean_Compressible_response(NeoHookean_Compressible,material_type_e,F,...
                    stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp);
                %%*************************************************************************
            case 3 %Neo_Hookean_Incompressible_response
                [Dtan,stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,stress_trial] =...
                    NeoHookean_Incompressible_response(NeoHookean_Incompressible,material_type_e,F,...
                    stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp);
                %%*************************************************************************
            case 5 %Ferrara-Pandolfi_Incompressible_response
                [Dtan,stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,J3,stress_trial] =...
                    FerraraPandolfi_Incomp_response(Ferrara_Pandolfi,material_type_e,F,...
                    stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,[0 1 0],[0 0 1]);
                %%*************************************************************************

        end
        %%-------------------------------------------------------------------------
        %%calcMaterialResponse ends.
        %%-------------------------------------------------------------------------
        [Bf] = getLinearBmatrix( dNdx , elemType );
        K_e=K_e+Bf'*Dtan*Bf*W(igp)*det(J);
        if Dtan==zeros(9)
            K_e=zeros(3*n);
            break;
        end

    end

    K(:,iel)=K_e(:);

    sctB(:,iel)       = sctrB;
    sctrB             = ones(length(sctrB),1)*sctrB;
    sctrBT            = sctrB';
    sctB_c(:,iel)     = sctrB(:);
    sctB_r(:,iel)     = sctrBT(:);

end

clear iel igp col sctB_r sctB_c sctrBT sctB
%--------------------------------------------------------------------------
%initialization ends
%--------------------------------------------------------------------------

iter=0;
converged=false;

r=zeros(TDOF,1);

while(~converged)

    iter=iter+1;

    Du=K\r;

    u=u+Du;

    %------------------------------------------------------------------
    %update and assemble begins
    %------------------------------------------------------------------
    fint=zeros(3*8,numelem) ;%Exclusive to H4 element
    K=zeros((3*8)^2,numelem);%Exclusive to H4 element

    Macro_stress=0;
    VoigtMacro_Dtan=0;
    LMacro_Dtan=zeros((3*8)*9,numelem);%Exclusive to B8 element;

    sctB    = zeros(3*8,numelem) ;
    sctB_r  = zeros((3*8)^2,numelem);
    sctB_c  = zeros((3*8)^2,numelem);
    sctB_L_r= zeros((3*8)*9,numelem);

    for iel=1:numelem

        sctr=element(iel,:);
        n=length(sctr);
        sctrB=zeros(1,3*n); %Element degrees of freedom
        sctrB(1:3:3*n)=3*sctr-2;
        sctrB(2:3:3*n)=3*sctr-1;
        sctrB(3:3:3*n)=3*sctr;
        elementu=u(sctrB);
        elementu=reshape(elementu,3,n);

        K_e=0;
        fint_e=0;
        material_type_e=material_type(iel,:);

        LMacro_Dtan_e=0;

        %Cycle for Gauss Points
        for igp=1:size(W,1)

            strain_gp=strain(iel,((igp-1)*6)+1:igp*6);
            strain_vmisses_gp=zeros(1,1);
            stress_gp=stress(iel,((igp-1)*6)+1:igp*6);
            stress_vmisses_gp=zeros(1,1);

            [N,dNdxi] = lagrange_basis(elemType,Q(igp,:));
            J=node(sctr,:)'*dNdxi;
            dNdx=dNdxi/J;
            F=macro_F+elementu*dNdx;

            %%-----------------------------------------------------------------------------
            %%calcMaterialResponse begins.
            %%-----------------------------------------------------------------------------
            Dtan=0;

            switch material_type_e(1,1)
                %%*************************************************************************
                case 1 %StVenant_Kirchhoff
                    [Dtan,stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,stress_trial]=...
                        StVenant_Kirchhoff_response(StVenant_Kirchhoff,material_type_e,F,...
                        stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,...
                        I,L,L1,IdyadI,0,0);
                    %%*************************************************************************
                case 2 %Neo_Hookean_Compressible_response
                    [Dtan,stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,stress_trial]=...
                        NeoHookean_Compressible_response(NeoHookean_Compressible,material_type_e,F,...
                        stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,...
                        I,L,L1,0,0);
                    %%*************************************************************************
                case 3 %Neo_Hookean_Incompressible_response
                    [Dtan,stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,stress_trial] =...
                        NeoHookean_Incompressible_response(NeoHookean_Incompressible,material_type_e,F,...
                        stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,...
                        I,L,L1,IdyadI,0,0);
                    %%*************************************************************************
                case 5 %Ferrara-Pandolfi_Incompressible_response
                    [Dtan,stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,J3,stress_trial] =...
                        FerraraPandolfi_Incomp_response_Macro(Ferrara_Pandolfi_Macro,material_type_e,...
                        reset_load,F,...
                        stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,...
                        I,Lz,L1,elementu,0,0,[0 1 0],[0 0 1]);
                    %%*************************************************************************
            end
            %%-------------------------------------------------------------------------
            %%calcMaterialResponse ends.
            %%-------------------------------------------------------------------------


            [Bf] = getLinearBmatrix( dNdx , elemType );

            fint_e=fint_e+Bf'*stress_trial*W(igp)*det(J);

            K_e=K_e+Bf'*Dtan*Bf*W(igp)*det(J);


            Macro_stress= Macro_stress+stress_trial*W(igp)*det(J);
            VoigtMacro_Dtan= VoigtMacro_Dtan+Dtan*W(igp)*det(J);
            LMacro_Dtan_e=LMacro_Dtan_e+Bf'*Dtan*W(igp)*det(J);

            strain(iel,((igp-1)*6)+1:igp*6)=strain_gp;
            strain_vmisses(iel,igp)=strain_vmisses_gp;
            stress(iel,((igp-1)*6)+1:igp*6)=stress_gp;
            stress_vmisses(iel,igp)=stress_vmisses_gp;

        end

        fint(:,iel)=fint_e;
        K(:,iel)=K_e(:);
        LMacro_Dtan(:,iel)=LMacro_Dtan_e(:);

        sctB(:,iel)       = sctrB;

        sctrB             = ones(length(sctrB),1)*sctrB;
        sctrBT            = sctrB';
        sctB_c(:,iel)     = sctrB(:);
        sctB_r(:,iel)     = sctrBT(:);

        sctrB_L           = ones(9,1)*sctB(:,iel)';
        sctrB_LT          = sctrB_L';
        sctB_L_r(:,iel)   = sctrB_LT(:);

    end

    sctB_L_c=ones(3*8*9,numelem);
    for lcol=1:9
        sctB_L_c(3*8*(lcol-1)+1:3*8*lcol,:)=sctB_L_c(3*8*(lcol-1)+1:3*8*lcol,:)*lcol;
    end
    LMacro_Dtan=sparse(sctB_L_r,sctB_L_c,LMacro_Dtan,TDOF,9);


    clear iel igp col sctB_r sctB_c sctrBT sctB sctB_L_c sctB_L_r sctrB_LT sctrB_L

    %------------------------------------------------------------------
    %update and assemble ends
    %------------------------------------------------------------------

    if iter==1
        scale=norm(fint);
    end

    r=-fint;

    res=norm(fint(act))/scale;


    if iter==1
        res=1;
    end

    if res<=global_toll
        converged=true;
    end

    if iter>max_iter
        load_break=true;
        break;
    end

end


stress_trial=Macro_stress/vol;


fluc_Dtan=LCond(Indexer(:,1)',:)'*K^(-1)*LCond(Indexer(:,1)',:);
Dtan=(VoigtMacro_Dtan-(fluc_Dtan))/vol;

Dtan_4th=zeros(3,3,3,3);
stress_trial_2nd=zeros(3);

idx1=[1 2 3 1 2 2 3 3 1];
idx2=[1 2 3 2 1 3 2 1 3];

for i=1:9      %voigt to 4th
    for j=1:9
        Dtan_4th(idx1(i),idx2(i),idx1(j),idx2(j))=Dtan(i,j);
    end
end
for i=1:9    %voigt to 2nd
    stress_trial_2nd(idx1(i),idx2(i))=stress_trial(i);
end

for i=1:9  % 4th to voigt
    for j=1:9
        Dtan(i,j)=Dtan_4th(idx1(i),idx2(i),idx1(j),idx2(j));
    end
end
for i=1:9 % 2nd to voigt
    stress_trial(i)=stress_trial_2nd(idx1(i),idx2(i));
end