function [SMAX,SMIN]=PLOTELEM(XYZ,LE,S,TIME,OPTION,SCALE)
%******************************************************************************
% Plot initial & deformed geometry, contour plot                              %
%******************************************************************************
%%                                                                            %
 [ND,NDIM]=size(XYZ); [NE,NEND]=size(LE);                   % Size of the model
 fid=fopen(S,'r');                                           % Open output file
 FACE = [1 4 3 2;1 5 8 4;1 2 6 5;2 3 7 6;3 4 8 7;5 6 7 8];   % 6 faces of solid
 if NDIM==2                                                          % 2D model
  NFACE=1;                                                           % One face
  T=[1.866 -0.5 -0.5 0.134;-0.5 0.134 1.866 -0.5;        % Stress extrapolation 
     0.134 -0.5 -0.5 1.866;-0.5 1.866 0.134 -0.5];                   % to nodes
 else                                                                % 3D model
  NFACE=6;                                                          % Six faces
  A=2.5490381056;B=-0.6830127018;C=0.1830127018;D=-0.0490381056;    % Constants
  T=[A B C B B C D C;B C D C A B C B;B C B A C D C B;C D C B B C B A;         %
     B A B C C B C D;C B C D B A B C;C B A B D C B C;D C B C C B A B]';       %
 end                                                  % End of dimension branch
 if OPTION~=-1                        % Plot deformed geometry or contour plots
  DISP=zeros(ND,NDIM);                                    % Nodal displacements
  STRS=zeros(NE,NEND);                                   % Element nodal stress
  S=zeros(NEND,1);                                   % Temporary stress storage
  while 1                                      % Loop for all lines in the file
   tline=fgetl(fid);                                              % Read a line
   if strfind(tline,'TIME =')                                  % Find TIME line
    t=sscanf(tline,'TIME = %f');                               % Read time data
    if t==TIME                                      % Process this time results
     for i=1:4, fgetl(fid); end                                  % Skip 4 lines
     for i=1:ND                                            % Loop for all nodes
      V=sscanf(fgetl(fid),'%d %f %f %f');              % Read displacement data
      DISP(V(1),:)=V(2:NDIM+1);                     % Store nodal displacements
     end                                                % End of all nodes loop
	 if OPTION<1, break; end                   % No need to process stress data
     for i=1:4, fgetl(fid); end                                  % Skip 4 lines
     for e=1:NE                                                  % Element loop
      fgetl(fid);                                             % Read empty line
      for i=1:NEND                                % Loop for integration points
       V=sscanf(fgetl(fid),'%f %f %f %f %f %f');             % Read stress data
       if OPTION==7                                % Calculate von Mises stress
        S(i)=sqrt((V(1)-V(2))^2+(V(2)-V(3))^2+(V(3)-V(1))^2 ...               %
                 +6*(V(4)^2+V(5)^2+V(6)^2))/sqrt(2);                          %
       else                                                  % Stress component
        S(i)=V(OPTION);                                      % Stress component
       end                                          % End of stress type branch
      end                                       % End of integration point loop
      STRS(e,:)=(T*S)';                           % Extrapolate stress to nodes
     end                                                  % End of element loop
     break;                                   % Stop processing the output file
    end                                           % End of matching time branch
   end                                             % End of current time branch
  end                                                % End of reading file loop
 end                                                    % End of results branch
 figure('Color','w'); axis equal; axis off; hold on           % Create a figure
 P=struct('Faces',FACE(1:NFACE,:),'FaceColor','w','EdgeColor','k');   % Setting
 if OPTION>=0, P.FaceColor='none';P.LineStyle=':'; end        % Dashed geometry
 for e=1:NE, P.Vertices=XYZ(LE(e,:)',:); patch(P); end  % Plot initial geometry                                              % Element loop
 SMAX=max(max(XYZ)); SMIN=min(min(XYZ));              % Max & min of dimensions
 if OPTION>=0                                          % Plot deformed geometry
  if isempty(SCALE), SCALE=0.2*max(max(abs(XYZ)))/(max(max(abs(DISP))));end   %
  P.LineWidth=1.5;P.LineStyle='-';P.FaceColor='w';               % Plot setting
  for e=1:NE                                                     % Element loop
   P.Vertices=XYZ(LE(e,:)',:)+SCALE*DISP(LE(e,:)',:);    % Deformed coordinates
   if OPTION>0, P.FaceColor='interp';P.FaceVertexCData=STRS(e,:)'; end % Stress
   patch(P);                                           % Plot deformed geometry
  end                                                     % End of element loop
  SMAX=max(max(DISP)); SMIN=min(min(DISP));        % Max & min of displacements
 end                                          % End of deformed geometry branch
 if OPTION>0                                          % For stress contour plot
   colorbar('horiz'); SMAX=max(max(STRS)); SMIN=min(min(STRS));% Contour legend
 end                                               % End of stress contour plot
 if NDIM==3, view(3); xlabel('X'); ylabel('Y'); zlabel('Z'); end      % 3D view
end