function ELAST3D(ETAN,NE,NDOF,XYZ,LE,DYNF)
%******************************************************************************
%  COMPUTING STIFFNESS & RESIDUAL FOR 3D 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+[-2 -1 0])=LE(IE,I)*NDOF+[-2 -1 0];          % DOF mapping index
  end                                               % End of element nodes loop
  DSP=reshape(DISPTD(IDOF),NDOF,8);    % NDOFxNENODE element total displacement
  for LX=1:2, for LY=1:2, for LZ=1:2              %LOOP OVER INTEGRATION POINTS
   E1=XG(LX); E2=XG(LY); E3=XG(LZ);                         % Integration point
   INTN = INTN + 1;                        % Integration point sequential index
   [~,SHPD,DET]=SHAPEL([E1 E2 E3], ELXY);       % Det, shape func & derivatives
   FAC=WGT(LX)*WGT(LY)*WGT(LZ)*DET;                        % Integraton weights
   DEPS=DSP*SHPD';                                           % Calculate strain
   DDEPS=[DEPS(1,1) DEPS(2,2) DEPS(3,3) ...                     % Strain vector
          DEPS(1,2)+DEPS(2,1) DEPS(2,3)+DEPS(3,2) DEPS(1,3)+DEPS(3,1)]';
   STRESS = ETAN*DDEPS;                                      % Calculate stress
   if UPDATE, SIGMA(:,INTN)=STRESS; continue; end              % Store stresses
   BM=zeros(6,24);                                 % displacement-strain matrix
   for I=1:8                                                % Loop for B matrix
    COL=I*NDOF+[-2 -1 0];                          % Colum indices for B matrix 
    BM(:,COL)=[SHPD(1,I) 0 0; 0 SHPD(2,I) 0; 0 0 SHPD(3,I);
          SHPD(2,I) SHPD(1,I) 0; 0 SHPD(3,I) SHPD(2,I); SHPD(3,I) 0 SHPD(1,I)];
   end                                                   % End of B matrix loop
   FORCE(IDOF) = FORCE(IDOF)-FAC*BM'*STRESS;         % Residual internal forces
   if LTAN                                        % Calculate tangent stiffness
    EKF = BM'*ETAN*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                                % End of integration points loop
 end                                                      % End of element loop
end