function MASS3D(NE,NDOF,XYZ,LE,DYN)
%******************************************************************************
%  MAIN PROGRAM COMPUTING MASS MATRIX AND INERTIAL VECTOR FOR                 %
%  3D HEXAHEDRAL ELEMENTS                                                     %
%******************************************************************************
%%
 global ACC FORCE GKF                                        % Global variables
 XG=[-1/sqrt(3), 1/sqrt(3)]; WGT=[1.0, 1.0];     % Integration points & weights
 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
  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
   [SHP, ~, DET] = SHAPEL([E1 E2 E3], ELXY); % Det & shape function derivatives
   FAC=DYN.RHO*WGT(LX)*WGT(LY)*WGT(LZ)*DET;                % Integraton weights
   FORCE(IDOF) = FORCE(IDOF)-FAC*SHP*(SHP'*ACC(IDOF));        % Sbtract Ma term
   GKF(IDOF,IDOF)=GKF(IDOF,IDOF)+FAC*(SHP*SHP');                  % Mass Matrix 
  end, end, end                                % End of integration points loop
 end                                                      % End of element loop
end                                                                           %