function [SF,GDSF,DET]=SHAPEL(XI,ELXY)
%******************************************************************************
% COMPUTING SHAPE FUNCTION, DERIVATIVES, AND DETERMINANT AT INTEGRATION POINTS
%******************************************************************************
%%
 XNODE=[-1  1  1 -1 -1  1  1 -1;          % Nodal coordinates of parent element
        -1 -1  1  1 -1 -1  1  1;
        -1 -1 -1 -1  1  1  1  1];
 [NNODE, NDOF]=size(ELXY); A=1/NNODE; ID=1:NDOF;    % Define # of Nodes and DOF
 SF=zeros(NNODE*NDOF,NDOF); DSF=zeros(NDOF,NNODE);   % Shape fns and derivative
 for I=1:NNODE                                         % Loop for element nodes
  XP=XNODE(1:NDOF,I)';                            % I-th column of XNODE matrix
  XI0=ones(1,NDOF)+XI.*XP;                 % Terms in each coordinate direction
  SF(NDOF*(I-1)+1:NDOF*I,1:NDOF)=A*prod(XI0)*eye(NDOF);      % Shape fn. matrix
  for J=1:NDOF                                              % Loop for each DOF
   DSF(J,I)=A*XP(J)*prod(XI0(setdiff(ID,J)));   % Derivative of shape functions
  end                                                         % End of DOF loop
 end                                               % End of elemenet nodes loop
 GJ=DSF*ELXY;                                                 % Jacobian matrix
 DET=det(GJ);                                  % Determinant of Jacobian matrix
 GDSF=GJ\DSF;       % Derivatives of shape functions w.r.t. spatial coordinates
end