function ELAST2D(ETAN,NE,NDOF,XYZ,LE,DYNF)
%******************************************************************************
% COMPUTE STIFFNESS & RESIDUAL FOR 2D PLANE-STRAIN LINEAR ELASTIC MATERIAL
%******************************************************************************
%%
 global LTAN UPDATE DISPTD FORCE GKF SIGMA
 %
 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 for 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
   DEPS=DSP*SHPD';                                           % Calculate strain
   DDEPS=[DEPS(1,1) DEPS(2,2) 0 DEPS(1,2)+DEPS(2,1) 0 0]';      % Strain vector
   STRESS = ETAN*DDEPS;                                      % Calculate stress
   if UPDATE, SIGMA(:,INTN)=STRESS; continue; end              % Store stresses
   BM=zeros(4,8);                                  % displacement-strain matrix
   for I=1:4                                                % Loop for B matrix
    COL=I*NDOF+[-1 0];                             % Colum indices for B matrix 
    BM(:,COL)=[SHPD(1,I) 0; 0 SHPD(2,I); 0 0; SHPD(2,I) SHPD(1,I)];
   end                                                   % End of B matrix loop
   FORCE(IDOF) = FORCE(IDOF)-FAC*BM'*STRESS(1:4);    % Residual internal forces
   if LTAN                                        % Calculate tangent stiffness
    EKF = BM'*ETAN(1:4,1:4)*BM;                     % Element tangent stiffness
    GKF(IDOF,IDOF)=GKF(IDOF,IDOF)+DYNF*FAC*EKF;          % Assemble to global K
   end                                        % End of tangent stiffness branch
  end, end                                     % End of integration points loop
 end                                                      % End of element loop
end