function  contour_qp_tau(inc_q)
Pts_wide=2;
% Copyright 2000, J.E. Akin.  All rights reserved.
% plot element Gauss point Tau magnitude as a surface 
% convert any type of mesh to a structured square mesh

if ( nargin == 0 )
  inc_q = 1 ;
end % if 

% x,y   ==  nodal coordinates of the FEA
% f     ==  nodal solution of the FEA
% X,Y   ==  nodal coordinates of the structured square mesh
% F     ==  interpolated solution for the structured square mesh
% NX    ==  number of nodes in x-direction for structured mesh 
% NY    ==  number of nodes in y-direction for structured mesh

load msh_bc_xyz.tmp;
np = size (msh_bc_xyz, 1);
fprintf ('Read %g mesh coordinate pairs \n', np)
pre_p = 1;
ns = size (msh_bc_xyz,2) - pre_p ;   % space dimension
if ( ns < 2 )
  error ('This is not a 2D mesh')
end % if not 2D data

pre_e  = 0 ;
load  msh_typ_nodes.tmp ; % nod_per_el nodes per element
nt         = size (msh_typ_nodes,1);           % number of elements
if ( nt == 0 )
  error ('Error missing file msh_typ_nodes.tmp')
end % if error
nod_per_el = size(msh_typ_nodes,2) - pre_e -1 ; %nodes per element

load el_qp_xyz_tau.tmp; % x,y followed by tau list
ng = size (el_qp_xyz_tau, 1) ;     % count points
nf = size (el_qp_xyz_tau, 2) - 2;  % count components
if ( ng == 0 )
  error ('Error missing file el_qp_xyz_tau.tmp')
end % if error
if ( nf <= 0 )
  error ('Error, no fluxes in file el_qp_xyz_tau.tmp')
end % if error

% NX, NY = 31 are default values, these can be altered by user
NX = 51; NY = 51;

x = el_qp_xyz_tau (1:inc_q:ng, 1);
y = el_qp_xyz_tau (1:inc_q:ng, 2);
f = el_qp_xyz_tau (1:inc_q:ng, 3);
% x(1)
% y(1)
AveTau = mean(f);
V_X=max(f);
fprintf ('Max Tau: %g ', max(f))
fprintf ('Min Tau: %g ', min(f))
fprintf ('Ave Tau: %g \n', AveTau)

xlin = linspace (min(x), max(x), NX);
ylin = linspace (min(y), max(y), NY);
[X, Y] = meshgrid (xlin, ylin);
F = griddata (x, y, f, X, Y, 'cubic');

clf
% V=[0.02,0.04,0.06,0.08,0.10];
% V=[0.04,0.08,0.12,0.16,0.20];
% V=[0.05,0.10,0.15,0.20,0.25];
% V=[0.10,0.20,0.40,0.60,0.80,1.00,1.20];
  V=[0.20,0.40,0.60,0.80,1.00,1.20];
 c = contour (X,Y,F,V,'LineWidth',Pts_wide);
% c = contour (X,Y,F);
axis('square')
%b clabel (c); 
clabel (c,'manual') ; % V);
%b clabel (c,V);
hold on
grid
Show=0;
% show actual qp tau values ?
if ( Show == 1 )
  for kount=1:ng
    t_text = sprintf ('%g', f(kount)) ;
    text (x(kount), y(kount), t_text, 'Rotation', 30)
  end % for qp loop
end % if show

%b xlabel ('X')
%b ylabel ('Y')
%xlabel (['X: ', int2str(nt),' Elements with ', ...
%        int2str(nod_per_el), ' nodes'])
%ylabel (['Y: ', int2str(np),' Nodes'])
 %b title (['Smoothed Tau values at ', ...
          %b int2str(ng), ' Gauss Points . (Max, Ave ', ...
          %b num2str(V_X), ', ', num2str(AveTau), ')']);
%       -depsc -tiff            % for an eps version
  %bprint ('-dpsc', ['contour_qp_tau'])
  hold off
  %bv_text = ['Created contour_qp_tau.ps'];
  %bfprintf (1,'%s', v_text) ; fprintf (1, ' \n' )

% end of contour_qp_tau