function HYPER2D(PROP,NE,NDOF,XYZ,LE,DYNF)
%******************************************************************************
%  MAIN PROGRAM COMPUTING GLOBAL STIFFNESS MATRIX RESIDUAL FORCE FOR          %
%  HYPERELASTIC MATERIAL MODELS                                               %
%******************************************************************************
%%                                                                            %
 global LTAN UPDATE DISPTD FORCE GKF SIGMA                   % Global variables
 XG=[-1/sqrt(3), 1/sqrt(3)]; WGT=[1.0, 1.0];     % Integration points & weights
 INTN=0;                    % Index for history variables (each integration pt)
 for IE=1:NE                             % Loop over elements. Assemble K and F
  ELXY=XYZ(LE(IE,:),:);                             % Element nodal coordinates
  IDOF=zeros(1,NDOF*size(LE,2));            % Local to global DOF mapping index
  for I=1:size(LE,2)                                   % Loop for element nodes
   IDOF(I*NDOF+[-1 0])=LE(IE,I)*NDOF+[-1 0];                % DOF mapping index
  end                                               % End of element nodes loop
  DSP=reshape(DISPTD(IDOF),NDOF,4);    % NDOFxNENODE element total displacement
  for LX=1:2, for LY=1:2                          %LOOP OVER INTEGRATION POINTS
   E1=XG(LX); E2=XG(LY);                                    % Integration point
   INTN = INTN + 1;                        % Integration point sequential index
   [~, SHPD, DET] = SHAPEL([E1 E2], ELXY);   % Det & shape function derivatives
   FAC=WGT(LX)*WGT(LY)*DET;                                % Integraton weights
   F=eye(3); F(1:2,1:2)=F(1:2,1:2)+DSP*SHPD';            % Deformation gradient
   [STRESS, DTAN] = MOONEY(F,PROP,LTAN); % Compute stress and tangent stiffness
   if UPDATE                          % Store int. pt. stresses to global array
    STRESS=CAUCHY(F, STRESS);        % Convert 2nd PK stress into Cauchy stress
    SIGMA(:,INTN)=STRESS;                    % Store stress to the global array
    continue;                                        % Do not calculate F and K
   end                                                   % End of update branch
   BN=zeros(4,8); BG=zeros(6,8);                   % displacement-strain matrix
   for I=1:4                                               % Element nodes loop
    COL=(I-1)*NDOF+1:(I-1)*NDOF+NDOF;               % DOF locations of the node
    BN(:,COL)=[SHPD(1,I)*F(1,1), SHPD(1,I)*F(2,1); % Displacement-strain matrix
               SHPD(2,I)*F(1,2), SHPD(2,I)*F(2,2);                            %
                              0,                0;                            %
               SHPD(1,I)*F(1,2)+SHPD(2,I)*F(1,1), ...                         %
               SHPD(1,I)*F(2,2)+SHPD(2,I)*F(2,1)];                            %
    BG(:,COL)=[SHPD(1,I) 0; SHPD(2,I) 0; 0 0;      % Displacement-strain matrix
               0 SHPD(1,I); 0 SHPD(2,I); 0 0];                                %
   end                                              % End of element nodes loop
   FORCE(IDOF) = FORCE(IDOF) - FAC*BN'*STRESS(1:4);  % Residual internal forces
   if LTAN                                        % Calculate tangent stiffness
    SIG=[STRESS(1) STRESS(4) STRESS(6);                         % Stress tensor
         STRESS(4) STRESS(2) STRESS(5);                                       %
         STRESS(6) STRESS(5) STRESS(3)];                                      %
    SHEAD=kron(eye(2),SIG);                          % Initial stiffness matrix
    EKF = BN'*DTAN(1:4,1:4)*BN + BG'*SHEAD*BG;    % Calculate tangent stiffness
    GKF(IDOF,IDOF)=GKF(IDOF,IDOF)+DYNF*FAC*EKF;          % Assemble to global K
   end                                                     % End of LTAN branch
  end, end                                     % End of integration points loop
 end                                                      % End of element loop
end                                                                           %
%
function STRESS=CAUCHY(F, S)
%******************************************************************************
%  CONVERT 2ND PK STRESS INTO CAUCHY STRESS                                   %
%******************************************************************************
%%                                                                            %
 PK=[S(1) S(4) S(6);S(4) S(2) S(5);S(6) S(5) S(3)];     % 2nd PK stress (Voigt)
 DETF = det(F);                           % Determinant of deformation gradient
 ST = F*PK*F'/DETF;                                      % Cauchy stress tensor
 STRESS=[ST(1,1) ST(2,2) ST(3,3) ST(1,2) ST(2,3) ST(1,3)]';    % Voigt notation
end                                                                           %