function  contour_el_tau ()
Pts_wide=2;
% Copyright 2000, J.E. Akin.  All rights reserved.
% convert any type of mesh to a structured square mesh

% 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_tau_value.tmp; % 
ng = size (el_tau_value, 1) ;     % count points
if ( ng == 0 )
  error ('Error missing file el_tau_value.tmp')
end % if error

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

% Put value at element center
x(1:nt) = 0 ; y(1:nt) = 0 ; f(1:nt) = 0 ;
for it = 1:nt ;
  f (it) = el_tau_value (it) ;
  t_nodes (1:nod_per_el) = 0 ;
  t_nodes = msh_typ_nodes (it, (pre_e+2):(nod_per_el+pre_e+1));
  count = sum ((t_nodes > 0)==1) ;
  if ( count == fix(nod_per_el/4)*4 ) % then a quad
    x(it) = sum ( msh_bc_xyz (t_nodes(1:4), (pre_p+1)) )/4 ;
    y(it) = sum ( msh_bc_xyz (t_nodes(1:4), (pre_p+2)) )/4 ;
  else % a triangle
    x(it) = sum ( msh_bc_xyz (t_nodes(1:3), (pre_p+1)) )/3 ;
    y(it) = sum ( msh_bc_xyz (t_nodes(1:3), (pre_p+2)) )/3 ;
  end % if element type
end % for each element
AveTau = mean (f);
[V_X, L_X] = max (f);
[V_N, L_N] = min (f);

%   get mesh x,y 
x_m = msh_bc_xyz (:, 2); y_m = msh_bc_xyz (:, 3);
xmax = max(x_m); xmin = min (x_m);
ymax = max(y_m); ymin = min (y_m);

xlin = linspace (xmin, xmax, NX);
ylin = linspace (ymin, ymax, NY);
[X, Y] = meshgrid (xlin, ylin);
% [X, Y] = meshgrid (xlin, ylin);
F = griddata (x, y, f, X, Y, 'cubic');

clf
axis ([xmin, xmax, ymin, ymax])   % set axes
axis('equal')
c = contour (X,Y,F,'LineWidth',Pts_wide);
% V=[0.00050, 0.00075, 0.0010, 0.0125, 0.0150, 0.0175]
clabel (c);
% clabel (c,V);
hold on
grid

xlabel (['X for ', int2str(nt),' Elements with ', ...
        int2str(nod_per_el), ' nodes'])
ylabel (['Y for ', int2str(np),' Nodes'])
% xlabel ('X'); ylabel ('Y');
  title (['Smoothed FEA Stabilization Value: min ', ...
       num2str(V_N), ', max ', num2str(V_X)]) % , ...
       % ', ave ', num2str(AveTau)])
%       -depsc -tiff            % for an eps version
  % print ('-dpsc', ['contour_el_tau_'])
  hold off
  % v_text = ['Created contour_el_tau_.ps'];
  % fprintf (1,'%s', v_text) ; fprintf (1, ' \n' )

% end of contour_el_tau